CODEWITHSUNDEEP
matlab
user3277001
user3277001

Reputation: 21

system of 6 differential equations

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

Upvotes: 1

Views: 179

Answers (1)

Rody Oldenhuis
Rody Oldenhuis

Reputation: 38052

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.

Upvotes: 2

Related Questions