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;
end

[Aquatris]

Something does not feel right.

• you call ode15s with 'ff' argument, I do not think this is right way to call it
• The Methanol_Water function returns xdot as 13 element vector if 'flag', and if not it is an 8x8 matrix, which is not a regular ODE problem.

Answer 1

[Torsten]

tf = 200; 
x0 = [0 0 0 0 2500 0 0.8 0 50 50 50 50 5]; 
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);  % changed and divided xdot(7) by x(5) instead.
M(7,7) = 1;
M(8,8) = 1;
opts = odeset('Mass', M, 'MassSingular', 'yes');
[t, x] = ode15s(@Methanol_Water, [0 tf], x0,opts);
plot(t,x)

function xdot=Methanol_Water(t,x)
    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));   % Divided by M(7,7) = x(5) instead  
    xdot(7)=(L*(x4-xB)-V*(yB-xB))/x(5);  
    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(:); 
end

[Deyun]

Your answer is very helpful !