How do you include a mass matrix in ode45?

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

[Matthew Hunt]

I have coded up the mass matrix as a function:

function M=mass(y,h,c)
N=length(y);
n=int((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

Would I make this a function handle for this when I use ode45?

[Torsten]

options = odeset('Mass',@(t,y)mass(y,h,c))

[Matthew Hunt]

Cheers

[Paul]

@Matthew Hunt

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.

Answer 1

[Jan]

Simply by a matrix division inside the function to be integrated:

dydt = M(t,y) \ f(t,y)

[Matthew Hunt]

Good grief, the answer can't be that simple and straight forward, can it?

[Matthew Hunt]

Just had a look, and they recommend not to do this for speed of calculation.

[VBBV]

Please look at this solution given in below link

https://in.mathworks.com/help/matlab/ref/ode15s.html#bu7kkx9

Looks like you have a similar question as here

https://in.mathworks.com/matlabcentral/answers/571207-coupled-odes-involving-mass-matrix?s_tid=answers_rc1-3_p3_Topic

[Jan]

@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.

Answer 2

[Steven Lord]

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.

[Matthew Hunt]

This has no information on how to implement the mass matrix.

[Steven Lord]

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.]