# How do you include a mass matrix in ode45?

16 Ansichten (letzte 30 Tage)
Matthew Hunt am 14 Mär. 2023
Bearbeitet: Matthew Hunt am 21 Mär. 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 Kommentare2 ältere Kommentare anzeigen2 ältere Kommentare ausblenden
Matthew Hunt am 17 Mär. 2023
Cheers
Paul am 17 Mär. 2023
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.

Melden Sie sich an, um zu kommentieren.

### Antworten (2)

Jan am 14 Mär. 2023
Bearbeitet: Jan am 15 Mär. 2023
Simply by a matrix division inside the function to be integrated:
dydt = M(t,y) \ f(t,y)
##### 4 Kommentare2 ältere Kommentare anzeigen2 ältere Kommentare ausblenden
VBBV am 15 Mär. 2023
Jan am 15 Mär. 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.

Melden Sie sich an, um zu kommentieren.

Steven Lord am 14 Mär. 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 KommentareKeine anzeigenKeine ausblenden
Matthew Hunt am 14 Mär. 2023
This has no information on how to implement the mass matrix.
Steven Lord am 14 Mär. 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.]

Melden Sie sich an, um zu kommentieren.

### Kategorien

Mehr zu Ordinary Differential Equations finden Sie in Help Center und File Exchange

R2021a

### Community Treasure Hunt

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

Start Hunting!

Translated by