How do you include a mass matrix in ode45?

22 次查看（过去 30 天）
Matthew Hunt
Matthew Hunt 2023-3-14
编辑： Matthew Hunt 2023-3-21
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 个评论
显示 2更早的评论
Paul
Paul 2023-3-17
Ran in:
@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.

请先登录，再进行评论。

请先登录，再回答此问题。

回答（2 个）

Jan
Jan 2023-3-14
编辑：Jan 2023-3-15
Simply by a matrix division inside the function to be integrated:
dydt = M(t,y) \ f(t,y)
  4 个评论
显示 2更早的评论
Jan
Jan 2023-3-15
@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.

请先登录，再进行评论。

Steven Lord
Steven Lord 2023-3-14
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 个评论
Matthew Hunt
Matthew Hunt 2023-3-14
This has no information on how to implement the mass matrix.
Steven Lord
Steven Lord 2023-3-14
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.]

请先登录，再进行评论。

类别

在 Help Center 和 File Exchange 中查找有关 Ordinary Differential Equations 的更多信息

标签

产品


版本

R2021a

Community Treasure Hunt

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

Start Hunting!

Translated by