Clear Filters
Clear Filters

How do you include a mass matrix in ode45?

21 views (last 30 days)
Matthew Hunt
Matthew Hunt on 14 Mar 2023
Edited: Matthew Hunt on 21 Mar 2023
ode45 solves a set of equations
The length of is 3N+1, I have a non identity mass matrix M, but I've not seen any real instructions on how to implement it. Something about odeset, but nothing that really clear on how to a) define the mass matrix, and b) once defined how it is implemented. I've tried to read the instructions but there is little to no instructions on how to do this.
So my question is: How do you go about including the mass matrix into the calculation?
I can defined a function that gives the mass matrix, but how do I include it into the calculation?
The code I wrote is:
clear;clc;
N=10;
c=1;
h=1/N;
tspan=[0 20];
y_0=zeros(1,3*N+1);
y_0(1:N)=1;
y_0(N+1:2*N)=0;
y_0(2*N+1:3*N)=1;
y_0(3*N+1)=1;
opts = odeset('Mass',@(y,h,c) mass(y,h,c));
[t,y] = ode45(@(t,y) sint(N,y),tspan,y_0,opts);
function dydt = sint(N,y)
%This is the RHS of the sintering equations
g=9.81;
alpha=0.1;
h=1/N;
kappa=1;
x=linspace(0,1,N);
dydt=zeros(3*N+1,1);
dydt(3*N+1)=y(2*N);
%These are the bulk equations
for i=2:N-1
dydt(i)=y(2*N)*(y(i+1)-y(i-1))/(2*h*y(3*N+1))-(y(3*N+1)/(2*h))*(y(i+1)*y(N+1+1)-y(i-1)*y(N+i-1));
A=(y(2*N)*y(i)/(2*h*y(3*N+1)))*(y(i+1)-y(i-1))-(y(i)*y(N+i)/(2*h*y(3*N+1)))*(y(N+i+1)-y(N+i-1))+1/(2*h*y(3*N+1))*(P_L(y(i+1))-P_L(y(i-1)));
C=(1/(h^2*y(3*N+1)^2))*0.5*(zeta(y(i+1))+zeta(y(i)))*y(N+i+1)-0.5*(zeta(y(i+1))+zeta(y(i))+zeta(y(i))+zeta(y(i-1)))*y(N+i)+0.5*(zeta(y(i+1))+zeta(y(i-1)))*y(N+i-1);
D=-g+alpha/(2*h*y(3*N+1))*(y(2*N+i+1)-y(2*N+i-1));
dydt(N+i)=A+C+D;
E=-((y(N+i)/y(3*N+1))-y(2*N)/y(3*N+1).^2)*((y(N+i+1)-y(N+i-1))/(2*h));
F=(kappa*y(i)/(h^2*y(3*N+1)))*(y(2*N+i+1)-2*y(2*N+i)+y(2*N+i-1))-(P_L(y(i))/y(3*N+1))*((y(N+i+1)-y(N+i-1))/(2*h));
G=(alpha*y(i)*y(2*N+i)/h^2)*(y(N+i)/y(3*N+1)^2-y(2*N)/y(3*N+1))*(y(N+i+1-2*y(N+i)+y(N+i-1)));
dydt(2*N+i)=E+F+G;
end
%This is the boundary condition for the density;
dydt(1)=0;
dydt(N)=y(2*N)*(y(N)-y(N-1))/(h*y(3*N+1));
%This is the boundary condition for the velocity
dydt(N+1)=(P_L(x(2))-P_L(x(1)))/(h*y(3*N+1))+(y(N+2)*(zeta(x(2))-zeta(x(1))))/(h^2*y(3*N+1)^2)-g+alpha*(y(2*N+2)-y(2*N+1))/(h*y(3*N+1));
A=y(2*N)*y(N)*(y(N)-y(N-1))/(h*y(3*N+1))+2*y(N)*y(2*N)*(P_L(y(N))+alpha*T_a)/(y(3*N+1)*zeta(y(N)))+(P_L(y(N))-P_L(y(N-1)))/(h*y(3*N+1));
B=(0.5*(3*zeta(y(N))-zeta(y(N-1)))*(y(2*N-1)-2*h*(P_L(y(N)))+alpha*T_a)/zeta(y(N))+0.25*y(2*N)*(7*zeta(y(N))-zeta(y(N-1)))+0.5*y(2*N-1)*(zeta(y(N))+zeta(y(N-1))))/(h^2*y(3*N+1)^2);
C=-g+alpha*(y(3*N)-y(3*N-1))/(h*y(3*N+1));
dydt(2*N)=A+B+C;
%This is the boundary condition for Temperature
dydt(2*N+1)=y(1)*(kappa*y(3*N+1)-2*kappa*T_a+kappa*(2*T_a-y(3*N+1)))/(h^(2)*y(3*N+1)^2)-P_L(y(1))*y(N+2)/h;
dydt(3*N)=y(N)*(kappa*(2*T_a)-y(3*N-1)-2*kappa*T_a+kappa*y(3*N-1))-P_L(y(N))*(y(2*N)-y(2*N-1))/(h*y(3*N+1));
end
function M=mass(y,h,c)
N=length(y);
n=floor((N-1)/3);
l=ones(1,n);
M=zeros(N,N);
%Insert the conservation of mass terms
M(1:n,1:n)=diag(l,0);
%Insert the coefficients for the conservation of momentum
M(n+1:2*n,n+1:2*n)=diag(y(1:n),0);
%Insert the coefficients for the conservation of energy
M(2*n+1:3*n,2*n+1:3*n)=diag(c*y(1:n),0); %Diagonal elements
%Compute the off diagonal elements
l_sub=ones(1,n-1)/(2*h);
M_sub=diag(l_sub,-1)-diag(i_sub,1);
M_sub(1,1)=1/h;
M_sub(1,2)=-1/h;
M_sub(n,n-1)=1/h;
M_sub(n,n)=-1/h;
M(2*n,n+1:2*n)=M_sub;
M(N,N)=1;
end
function y=P_L(x)
gamma=1;
r_0=1;
y=3*gamma*(1-x).^2/r_0;
end
function y=zeta(x)
eta_0=1;
y=2*(1-x).^3*eta_0/(3*x);
end
function y=T_a(t)
y=1;
end
  4 Comments
Paul
Paul on 17 Mar 2023
Can you show how mass(h,y,c) works for example input data? I'm particularly curious about this line
N = 5
N = 5
n=int((N-1)/3);
Incorrect number or types of inputs or outputs for function 'int'.
Unless you've defined your own function called int, I don't understand how this will run because neither function int that comes with Matlab will apply.

Sign in to comment.

Answers (2)

Jan
Jan on 14 Mar 2023
Edited: Jan on 15 Mar 2023
Simply by a matrix division inside the function to be integrated:
dydt = M(t,y) \ f(t,y)
  4 Comments
Jan
Jan on 15 Mar 2023
@Matthew Hunt: There is no other option for ODE45. The DAE integrators allows to define a mass matrix. Then e.g. an approximated mass matrix can be applied to multiple evaluatione of a predictor-corrector method. But in ODE45 the only difference would be, that this matrix division is performed internally instead of manually.

Sign in to comment.


Steven Lord
Steven Lord on 14 Mar 2023
The "Summary of ODE Examples and Files" section on this documentation page lists a number of examples and indicates which of the ODE options each demonstrates. If the only option you need to set is the Mass matrix I recommend reading through the published batonode example (this documentation page) for an example of how to translate the mathematical form of your ODEs with a mass matrix into the code needed to solve it.
  2 Comments
Matthew Hunt
Matthew Hunt on 14 Mar 2023
This has no information on how to implement the mass matrix.
Steven Lord
Steven Lord on 14 Mar 2023
The "Code Equations" section on the batonode example page (the second link I posted) most certainly does show the mathematical form of the mass matrix, the code written to evaluate that mass matrix, and the code to specify that mass matrix in the ode45 solver call. [Well, okay, the code that solves the system with the mass matrix is actually in the "Solve System of Equations" section on that page not the "Code Equations" section.]

Sign in to comment.

Products


Release

R2021a

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!