How to use ode15s to solve a stiff ode with mass?
Now I get a matlab code using ode15s to solve a function with singular Mass. My question is:
what does it mean (there is no function called ‘M’ in my dir)
opts = odeset(‘Mass’, ‘M’, ‘MassSingular’, ‘yes’);
and why in the ode function, there is an extra input called ‘flag’ ?
The main code is as following:
tf = 200;
x0 = [0 0 0 0 2500 0 0.8 0 50 50 50 50 5];
opts = odeset(‘Mass’, ‘M’, ‘MassSingular’, ‘yes’);
[t, x] = ode15s(‘ff’, [0 tf], x0,opts);
function xdot=Methanol_Water(t,x)
if isempty(flag),
P=1;
x1=x(1);
x2=x(2);
x3=x(3);
x4=x(4);
x0=x(6);
xB=x(7);
T1=x(9);
T2=x(10);
T3=x(11);
T4=x(12);
TB=x(13);
D1=-4617.8;
C1=13.676;
D2=-5042.6;
C2=13.519;
A=0.85;
B=0.48;
R=10;
V=10;
Di=V/(R+1);
L=V-Di;
y1=exp(A*(1-x1^2)/(1-x1+A*x1/B)^2)*exp(C1+D1/(273.15+T1))*x1/P;
y2=exp(A*(1-x2^2)/(1-x2+A*x2/B)^2)*exp(C1+D1/(273.15+T2))*x2/P;
y3=exp(A*(1-x3^2)/(1-x3+A*x3/B)^2)*exp(C1+D1/(273.15+T3))*x3/P;
y4=exp(A*(1-x4^2)/(1-x4+A*x4/B)^2)*exp(C1+D1/(273.15+T4))*x4/P;
yB=exp(A*(1-xB^2)/(1-xB+A*xB/B)^2)*exp(C1+D1/(273.15+TB))*xB/P;
MH0=100;
MH=20;
xdot(1)=1/MH*(L*(x0-x1)+V*(y2-y1));
xdot(2)=1/MH*(L*(x1-x2)+V*(y3-y2));
xdot(3)=1/MH*(L*(x2-x3)+V*(y4-y3));
xdot(4)=1/MH*(L*(x3-x4)+V*(yB-y4));
xdot(5)=L-V;
xdot(6)=(V/MH0*y1-(L+Di)/MH0*x0);
xdot(7)=(L*(x4-xB)-V*(yB-xB));
xdot(8)=Di;
xdot(9)=P-exp(A*(1-x1^2)/(1-x1+A*x1/B)^2)*exp(C1+D1/(273.15+T1))*x1-exp(A*x1^2/(x1+A*(1-x1)/B)^2)*exp(C2+D2/(273.15+T1))*(1-x1);
xdot(10)=P-exp(A*(1-x2^2)/(1-x2+A*x2/B)^2)*exp(C1+D1/(273.15+T2))*x2- exp(A*x2^2/(x2+A*(1-x2)/B)^2)*exp(C2+D2/(273.15+T2))*(1-x2);
xdot(11)=P-exp(A*(1-x3^2)/(1-x3+A*x3/B)^2)*exp(C1+D1/(273.15+T3))*x3-exp(A*x3^2/(x3+A*(1-x3)/B)^2)*exp(C2+D2/(273.15+T3))*(1-x3);
xdot(12)=P-exp(A*(1-x4^2)/(1-x4+A*x4/B)^2)*exp(C1+D1/(273.15+T4))*x4- exp(A*x4^2/(x4+A*(1-x4)/B)^2)*exp(C2+D2/(273.15+T4))*(1-x4);
xdot(13)=P-exp(A*(1-xB^2)/(1-xB+A*xB/B)^2)*exp(C1+D1/(273.15+TB))*xB- exp(A*xB^2/(xB+A*(1-xB)/B)^2)*exp(C2+D2/(273.15+TB))*(1-xB);
xdot = xdot(:);
else
M = zeros(13,13);
M(1,1) = 1;
M(2,2) = 1;
M(3,3) = 1;
M(4,4) = 1;
M(5,5) = 1;
M(6,6) = 1;
M(7,7) = x(5);
M(8,8) = 1;
xdot = M;
endNow I get a matlab code using ode15s to solve a function with singular Mass. My question is:
what does it mean (there is no function called ‘M’ in my dir)
opts = odeset(‘Mass’, ‘M’, ‘MassSingular’, ‘yes’);
and why in the ode function, there is an extra input called ‘flag’ ?
The main code is as following:
tf = 200;
x0 = [0 0 0 0 2500 0 0.8 0 50 50 50 50 5];
opts = odeset(‘Mass’, ‘M’, ‘MassSingular’, ‘yes’);
[t, x] = ode15s(‘ff’, [0 tf], x0,opts);
function xdot=Methanol_Water(t,x)
if isempty(flag),
P=1;
x1=x(1);
x2=x(2);
x3=x(3);
x4=x(4);
x0=x(6);
xB=x(7);
T1=x(9);
T2=x(10);
T3=x(11);
T4=x(12);
TB=x(13);
D1=-4617.8;
C1=13.676;
D2=-5042.6;
C2=13.519;
A=0.85;
B=0.48;
R=10;
V=10;
Di=V/(R+1);
L=V-Di;
y1=exp(A*(1-x1^2)/(1-x1+A*x1/B)^2)*exp(C1+D1/(273.15+T1))*x1/P;
y2=exp(A*(1-x2^2)/(1-x2+A*x2/B)^2)*exp(C1+D1/(273.15+T2))*x2/P;
y3=exp(A*(1-x3^2)/(1-x3+A*x3/B)^2)*exp(C1+D1/(273.15+T3))*x3/P;
y4=exp(A*(1-x4^2)/(1-x4+A*x4/B)^2)*exp(C1+D1/(273.15+T4))*x4/P;
yB=exp(A*(1-xB^2)/(1-xB+A*xB/B)^2)*exp(C1+D1/(273.15+TB))*xB/P;
MH0=100;
MH=20;
xdot(1)=1/MH*(L*(x0-x1)+V*(y2-y1));
xdot(2)=1/MH*(L*(x1-x2)+V*(y3-y2));
xdot(3)=1/MH*(L*(x2-x3)+V*(y4-y3));
xdot(4)=1/MH*(L*(x3-x4)+V*(yB-y4));
xdot(5)=L-V;
xdot(6)=(V/MH0*y1-(L+Di)/MH0*x0);
xdot(7)=(L*(x4-xB)-V*(yB-xB));
xdot(8)=Di;
xdot(9)=P-exp(A*(1-x1^2)/(1-x1+A*x1/B)^2)*exp(C1+D1/(273.15+T1))*x1-exp(A*x1^2/(x1+A*(1-x1)/B)^2)*exp(C2+D2/(273.15+T1))*(1-x1);
xdot(10)=P-exp(A*(1-x2^2)/(1-x2+A*x2/B)^2)*exp(C1+D1/(273.15+T2))*x2- exp(A*x2^2/(x2+A*(1-x2)/B)^2)*exp(C2+D2/(273.15+T2))*(1-x2);
xdot(11)=P-exp(A*(1-x3^2)/(1-x3+A*x3/B)^2)*exp(C1+D1/(273.15+T3))*x3-exp(A*x3^2/(x3+A*(1-x3)/B)^2)*exp(C2+D2/(273.15+T3))*(1-x3);
xdot(12)=P-exp(A*(1-x4^2)/(1-x4+A*x4/B)^2)*exp(C1+D1/(273.15+T4))*x4- exp(A*x4^2/(x4+A*(1-x4)/B)^2)*exp(C2+D2/(273.15+T4))*(1-x4);
xdot(13)=P-exp(A*(1-xB^2)/(1-xB+A*xB/B)^2)*exp(C1+D1/(273.15+TB))*xB- exp(A*xB^2/(xB+A*(1-xB)/B)^2)*exp(C2+D2/(273.15+TB))*(1-xB);
xdot = xdot(:);
else
M = zeros(13,13);
M(1,1) = 1;
M(2,2) = 1;
M(3,3) = 1;
M(4,4) = 1;
M(5,5) = 1;
M(6,6) = 1;
M(7,7) = x(5);
M(8,8) = 1;
xdot = M;
end Now I get a matlab code using ode15s to solve a function with singular Mass. My question is:
what does it mean (there is no function called ‘M’ in my dir)
opts = odeset(‘Mass’, ‘M’, ‘MassSingular’, ‘yes’);
and why in the ode function, there is an extra input called ‘flag’ ?
The main code is as following:
tf = 200;
x0 = [0 0 0 0 2500 0 0.8 0 50 50 50 50 5];
opts = odeset(‘Mass’, ‘M’, ‘MassSingular’, ‘yes’);
[t, x] = ode15s(‘ff’, [0 tf], x0,opts);
function xdot=Methanol_Water(t,x)
if isempty(flag),
P=1;
x1=x(1);
x2=x(2);
x3=x(3);
x4=x(4);
x0=x(6);
xB=x(7);
T1=x(9);
T2=x(10);
T3=x(11);
T4=x(12);
TB=x(13);
D1=-4617.8;
C1=13.676;
D2=-5042.6;
C2=13.519;
A=0.85;
B=0.48;
R=10;
V=10;
Di=V/(R+1);
L=V-Di;
y1=exp(A*(1-x1^2)/(1-x1+A*x1/B)^2)*exp(C1+D1/(273.15+T1))*x1/P;
y2=exp(A*(1-x2^2)/(1-x2+A*x2/B)^2)*exp(C1+D1/(273.15+T2))*x2/P;
y3=exp(A*(1-x3^2)/(1-x3+A*x3/B)^2)*exp(C1+D1/(273.15+T3))*x3/P;
y4=exp(A*(1-x4^2)/(1-x4+A*x4/B)^2)*exp(C1+D1/(273.15+T4))*x4/P;
yB=exp(A*(1-xB^2)/(1-xB+A*xB/B)^2)*exp(C1+D1/(273.15+TB))*xB/P;
MH0=100;
MH=20;
xdot(1)=1/MH*(L*(x0-x1)+V*(y2-y1));
xdot(2)=1/MH*(L*(x1-x2)+V*(y3-y2));
xdot(3)=1/MH*(L*(x2-x3)+V*(y4-y3));
xdot(4)=1/MH*(L*(x3-x4)+V*(yB-y4));
xdot(5)=L-V;
xdot(6)=(V/MH0*y1-(L+Di)/MH0*x0);
xdot(7)=(L*(x4-xB)-V*(yB-xB));
xdot(8)=Di;
xdot(9)=P-exp(A*(1-x1^2)/(1-x1+A*x1/B)^2)*exp(C1+D1/(273.15+T1))*x1-exp(A*x1^2/(x1+A*(1-x1)/B)^2)*exp(C2+D2/(273.15+T1))*(1-x1);
xdot(10)=P-exp(A*(1-x2^2)/(1-x2+A*x2/B)^2)*exp(C1+D1/(273.15+T2))*x2- exp(A*x2^2/(x2+A*(1-x2)/B)^2)*exp(C2+D2/(273.15+T2))*(1-x2);
xdot(11)=P-exp(A*(1-x3^2)/(1-x3+A*x3/B)^2)*exp(C1+D1/(273.15+T3))*x3-exp(A*x3^2/(x3+A*(1-x3)/B)^2)*exp(C2+D2/(273.15+T3))*(1-x3);
xdot(12)=P-exp(A*(1-x4^2)/(1-x4+A*x4/B)^2)*exp(C1+D1/(273.15+T4))*x4- exp(A*x4^2/(x4+A*(1-x4)/B)^2)*exp(C2+D2/(273.15+T4))*(1-x4);
xdot(13)=P-exp(A*(1-xB^2)/(1-xB+A*xB/B)^2)*exp(C1+D1/(273.15+TB))*xB- exp(A*xB^2/(xB+A*(1-xB)/B)^2)*exp(C2+D2/(273.15+TB))*(1-xB);
xdot = xdot(:);
else
M = zeros(13,13);
M(1,1) = 1;
M(2,2) = 1;
M(3,3) = 1;
M(4,4) = 1;
M(5,5) = 1;
M(6,6) = 1;
M(7,7) = x(5);
M(8,8) = 1;
xdot = M;
end ode15s, mass MATLAB Answers — New Questions