system of 6 differential equations

Question

I am trying to solve a system of 6 differential equations using matlab. I created a set of 6 differential equations as follows in a function m file named as Untitled.m

function ydot=Untitled(t,y)
ydot = zeros(6,1);
%y(1)=A
%y(2)=B
%y(3)=C
%y(4)=D
%y(5)=P
%y(6)=T;

A=0.50265;
k11=(333/106.7)*1.15*1000*exp(-59660/(8.314*960));
k31=(333/40)*73.6*exp(-47820/(8.314*960));
k32=(333/14.4)*1.79*exp(-30950/(8.314*960));
k21=(106.7/40)*426*exp(-68830/(8.314*960));
k22=(106.7/14.4)*0.000599*exp(-57740/(8.314*960));
Pcat=1450;
g=9.81;
%phi=exp(-(59100*exp(-67210/(8.314*y(6))))*((0.50265*1450*33)/143.64));
H11=393000;
H31=795000;
H32=1200000;
H21=1150000;
H22=151000;
E=1-((285.765*17.56)/((6.1*1450)+(17.56*285.765)));
Fcat=143.64;
Cpcat=1087;
%Cp=1000*(y(1)*3.3+y(2)*3.3+y(3)*3.3+y(4)*1.087);
F=19.95;




ydot(1)= -(A*(1-E)*Pcat*(k11+k31+k32)*(exp(-(59100*exp(-67210/(8.314*y(6))))*((0.50265*1450*33)/143.64)))*y(1)*y(1))/F;
ydot(2)=  (A*(1-E)*Pcat*(k11*y(1)*y(1)-(k21+k22)*y(2))*(exp(-(59100*exp(-67210/(8.314*y(6))))*((0.50265*1450*33)/143.64))))/F;
ydot(3)=  (A*(1-E)*Pcat*(k31*y(1)*y(1)+ k21*y(2))*(exp(-(59100*exp(-67210/(8.314*y(6))))*((0.50265*1450*33)/143.64))))/F;
ydot(4)=  (A*(1-E)*Pcat*(k32*y(1)*y(1)+k22*y(2))*(exp(-(59100*exp(-67210/(8.314*y(6))))*((0.50265*1450*33)/143.64))))/F;
ydot(5)= -(Pcat*g*(1-E));
%dydz(6)= ((1-E)*Pcat*A*(((k11+k31+k32)*phi*y(1)*y(1)*-H1)+ ((k11*y(1)*y(1)-(k21+k22)*y(2))*phi*-H2)+((k31*y(1)*y(1)+ k21*y(2))*phi*H3)+((k32*y(1)*y(1)+k22*y(2))*phi*H4)))/(F*Cp+Fcat*Cpcat)
ydot(6) = (((exp(-(59100*exp(-67210/(8.314*y(6))))*((0.50265*1450*33)/143.64)))*(1-E)*Pcat*A)*(y(1)*y(1)*((k11*H11)+(k31*H31)+(k32*H32))+y(2)*((k21*H21)+(k22*H22)))/((F*(1000*(y(1)*3.3+y(2)*3.3+y(3)*3.3+y(4)*1.087)))+(Fcat*Cpcat)));

then i have created another file for using the ODE45 solver to solve the equations

function ode()

options = odeset('RelTol',1e-1,'AbsTol',[1e-1 1e-1 1e-1 1e-1 1e-1 1e-1 ]);
Y0=[1.0;0.0;0.0;0.0;180000.0;959.81];
zspan=0:0.1:33;

[z,y]= ode45(@Untitled,zspan,Y0,options);

figure
hold on
plot(z,y(:,5));

the code is compiling but not solving the system of differential equations and only giving straight line graph at initial value of variables. eg initial value of y(1) is 1, so i get a graph at y=1.

can anyone tell me what is going wrong

Untitled.m gives a warning that "input argument z might be unused, but a later input argument is used "; so i guess it is not getting the diff equations properly...but i have no idea what to do about it...please help — user3277001
That just means that t isn't being used in calculating the derivatives. I tested Untitled(0,Y0) and got ydot=[0;0;0;0;-5148;0] so it looks like they should be straight lines apart from y(5). Check your derivatives! — David

Rody Oldenhuis Rody Oldenhuis · Accepted Answer · 2014-02-12T09:15:15

Seems your derivative is badly scaled, or just plain incorrect.

The 5th element of ydot is

-(Pcat*g*(1-E));

which contains only constants and is therefore constant. Therefore, that will result in a straight line with negative slope.

The other 5 elements of ydot contain terms like

exp(-59100*exp(-67210/(...))*(...))

That will almost certainly underflow, meaning, the value is too small to be expressed in double precision numbers, and will therefore be equal to 0. Those will also result in straight lines, this time horizontal and at a value equal to their initial value.

The value exp(-745.1332...) (equal to 4.940656458412465e-324) is the lowest you'll be able to go. And this value is already tricky, since it's degenerate (exp(-724) < realmin).

So either you messed up your equations, or you need to re-scale the equations so that the probability of over/underflow inside the solution space is reduced to zero.

system of 6 differential equations

1 Answers