How do you include a mass matrix in ode45?

조회 수: 6(최근 30일)
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표시 이전 댓글 수: 3숨기기 이전 댓글 수: 3
Paul 2023년 3월 17일
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 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표시 이전 댓글 수: 3숨기기 이전 댓글 수: 3
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 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표시 이전 댓글 수: 1숨기기 이전 댓글 수: 1
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.]

댓글을 달려면 로그인하십시오.

범주

Find more on General Applications in Help Center and File Exchange

R2021a

Community Treasure Hunt

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

Start Hunting!

Translated by