Solve a system of four differential equations

Hello!

My question might be very simple but I am very new to Matlab. I am trying to solve a system of 4 differential equations and then plot the results. After searching online and have some help I ended up on the following code:

syms S(t) P(t) c(t) z(t) A b g h v k l u d r w t Y 
eqn1 = diff(S,t) == (A*exp(h*t)*b*S)^(1-g) * z^g - v*z - c - A*exp(h*t)*b*S*(k*exp(h*t*(-1))+1);
eqn2 = diff(P,t) == u*z - d*P + (1-b)*S;
eqn3 = diff(c,t) == ( ((A*exp(h*t)*b) * ((1-g)*b*z^g*(A*exp(h*t)*S*b)^(-g)*(A*exp(h*t)+S) - 1 - k*exp(l*t*(-1))))  + (((1-b)/u) * (v - g*z^(g-1)*(A*exp(h*t)*S*b)^(1-g))) - r ) * c ;
eqn4 = diff(z,t) == ( ( (g*(A*exp(h*t)*S*b)^(1-g)*z^(g-1) - v)/(g*(g-1)*(A*exp(h*t)*S*b)^(1-g)*z^(g-2)) ) * ( d + ((c*u*(1-w))/(P*w*(v-g*(A*exp(h*t)*S*b)^(1-g)*z^(g-1)))) + ((A*exp(h*t)*b) * ((1-g)*z^g*(A*exp(h*t)*b*S)^(-g)*(A*exp(h*t)*b + b*S) - 1 - (k*exp(l*t*(-1))))) + (((1-b)*(v-g*(A*exp(h*t)*b*S)^(1-g)*z^(g-1)))/(u)) ) ) + ( ( (b*(A*exp(h*t)+S)*z)/(A*exp(h*t)*S) ) * ( S*A*exp(h*t)*h + (A*exp(h*t)) * ((A*exp(h*t)*b*S)^(1-g)*z^g - v*z - c - (1+(k*exp(l*t*(-1))))*A*exp(h*t)*b*S) ) );
[VF,Subs] = odeToVectorField(eqn1,eqn2,eqn3,eqn4);
odefcn = matlabFunction(VF, 'Vars',{t,A,Y,b,d,g,h,k,l,r,u,v,w});
A=3;
b=0.4;
g=0.3;
h=0.3;
v=2;
k=4;
l=0.7;
u=10;
d=0.2;
r=0.5;
w=0.7;
tspan = [0 10];
y0 = [100 60 30 50];    
[t,Y] = ode45(@(t,Y) odefcn(t,A,Y,b,d,g,h,k,l,r,u,v,w) , tspan, y0);
S_1 = Y(:,1);	
P_1 = Y(:,2);
c_1 = Y(:,3);
z_1 = Y(:,4);
AF_1 = A*exp(h*t);
x_1 = b*S.*AF;
figure(1)
plot(c,x_1)
figure(2)
plot(x,P_1)

But I am unsure yet for what it really does. If I am correct the "odeToVectorField" function takes the 4 equations and set them in one simpler system, the "matlabFunction" converts this system to a matlab function and the "ode45" solves the system. So basically, from the "ode45" I will get back four columns with the solution for each variable. My question is what is the order of the results. I mean, does the first column of "ode45" give the solution for the variable of eqn1 (which is S), the second column the solution for the variable of eqn2 (which is P) and so on? Or in other words, the order of the results of "ode45" is in accordance to the order "odeToVectorField" takes the equations?

Thanks a lot in advance!

Answer 1

[Star Strider]

Youor description is essentially correct. However, odeToVectorField does not ‘simplify’ the system, instead it converts it to a vector of first-order differential equations that — after converting them to an anonymous function — the MATLAB numeric ODE integration functions can use.

‘My question is what is the order of the results.’

That is provided in the ‘Subs’ output of odeToVectorField.

Taking the posted code only as far as the ode45 call, and plotting the results demonstrates this —

syms S(t) P(t) c(t) z(t) A b g h v k l u d r w t Y 
eqn1 = diff(S,t) == (A*exp(h*t)*b*S)^(1-g) * z^g - v*z - c - A*exp(h*t)*b*S*(k*exp(h*t*(-1))+1);
eqn2 = diff(P,t) == u*z - d*P + (1-b)*S;
eqn3 = diff(c,t) == ( ((A*exp(h*t)*b) * ((1-g)*b*z^g*(A*exp(h*t)*S*b)^(-g)*(A*exp(h*t)+S) - 1 - k*exp(l*t*(-1))))  + (((1-b)/u) * (v - g*z^(g-1)*(A*exp(h*t)*S*b)^(1-g))) - r ) * c ;
eqn4 = diff(z,t) == ( ( (g*(A*exp(h*t)*S*b)^(1-g)*z^(g-1) - v)/(g*(g-1)*(A*exp(h*t)*S*b)^(1-g)*z^(g-2)) ) * ( d + ((c*u*(1-w))/(P*w*(v-g*(A*exp(h*t)*S*b)^(1-g)*z^(g-1)))) + ((A*exp(h*t)*b) * ((1-g)*z^g*(A*exp(h*t)*b*S)^(-g)*(A*exp(h*t)*b + b*S) - 1 - (k*exp(l*t*(-1))))) + (((1-b)*(v-g*(A*exp(h*t)*b*S)^(1-g)*z^(g-1)))/(u)) ) ) + ( ( (b*(A*exp(h*t)+S)*z)/(A*exp(h*t)*S) ) * ( S*A*exp(h*t)*h + (A*exp(h*t)) * ((A*exp(h*t)*b*S)^(1-g)*z^g - v*z - c - (1+(k*exp(l*t*(-1))))*A*exp(h*t)*b*S) ) );
[VF,Subs] = odeToVectorField(eqn1,eqn2,eqn3,eqn4);
odefcn = matlabFunction(VF, 'Vars',{t,A,Y,b,d,g,h,k,l,r,u,v,w});
A=3;
b=0.4;
g=0.3;
h=0.3;
v=2;
k=4;
l=0.7;
u=10;
d=0.2;
r=0.5;
w=0.7;
tspan = [0 10];
y0 = [100 60 30 50];    
[t,Y] = ode15s(@(t,Y) odefcn(t,A,Y,b,d,g,h,k,l,r,u,v,w) , tspan, y0);

figure
hp = plot(t,real(Y));
hold on
plot(t, imag(Y),'--')
hold off
grid
legend(hp, string(Subs))

Also, note that ode15s is more appropriate for this, since the system appears to be ‘stiff’.

.

[Myrto Kasioumi]

Thank you very much, I will use ode15s. I just want to be sure about something: When I am running the code up to the "odeToVectorField(eqn1,eqn2,eqn3,eqn4)", it returns a Subs matrix which has 4 rows and they are filled like that: P, S, c, z. Does that mean that when I am solving the system later using ode15s, the first column will be the P variable?

Also, is it wrong to save the ode15s columns into different variables like so:

S_1 = Y(:,1);	
P_1 = Y(:,2);
c_1 = Y(:,3);
z_1 = Y(:,4);

and then plot them using those names instead of plotting them right after the ode?

Thank you very much!

[Star Strider]

My pleasure!

‘Does that mean that when I am solving the system later using ode15s, the first column will be the P variable?’

Yes. The columns are in the order of the ‘Subs’ vector, so the first column of ‘Y’ corresponds to ‘P’, the second column corresponds to ‘S’ and so for the others.

Note — This is also the order of the initial conditions vector.

.

[Myrto Kasioumi]

Thank you so much. Enjoy your weekend!

[Star Strider]

My pleasure!

If my Answer helped you solve your problem, please Accept it!

.