MATLAB Answers

How to solve and plot system of nonlinear differential equations?

243 views (last 30 days)
Daniel Jovin
Daniel Jovin on 22 Feb 2018
Commented: thanos tria on 26 Mar 2019
I'm trying to recreate graphs from a modeling paper by plotting a system of differential equations in MatLab. Unfortunately, I don't have much MatLab experience if any. I've found other questions on systems of nonlinear equations asked in MatLab answers and have managed to produce a plot for my own system, but this plot is not the same as the one in the paper I'm using.
https://imgur.com/a/hBQ3z This should show the differential equations I'm using, the parameter values, the graph I'm trying to recreate, and the plot I ended up with.
Here is the code I used to set up my system and plot the equations.
syms p(t) m(t) l(t) T Y
Eqns = [diff(p(t),t) == (3*p(t)*(1-p(t)))-(30*(p(t)*m(t))); diff(m(t),t) == ((25*p(t)+l(t))*m(t)*(1-m(t)))-m(t); diff(l(t),t) == (15*(1+tanh(2*m(t)-2)))-l(t)]
[DEsys,Subs] = odeToVectorField(Eqns);
DEFcn = matlabFunction(DEsys, 'Vars',{T,Y});
tspan = [0,25];
y0 = [0.01 0.05 0.539];
[T,Y] = ode45(DEFcn, tspan, y0);
figure(1)
plot(T,Y)
legend('p(t)','m(t)','l(t)')
grid
These are the same equations shown in the imgur album which I tried to recreate through diff(p(t),t), diff(m(t),t), and diff(l(t),t)
Here are the parameter values and graph I'm trying to recreate
And this is the graph I ended up with
I based my code off of a response provided to a question asked here (https://www.mathworks.com/matlabcentral/answers/365203-how-to-solve-a-system-of-nonlinear-differential-equations), but I would really appreciate anyone that could help me solve my plotting problems and show me what I did incorrectly and how to get the right plot.
Thanks.

  0 Comments

Sign in to comment.

Accepted Answer

Daniel Jovin
Daniel Jovin on 26 Feb 2018
Edited: Daniel Jovin on 26 Feb 2018
% code
kp = 3;
kpm = 30;
kmp = 25;
klm = 15;
kl = 1;
f = @(t,y) [(kp*y(1)*(1-y(1)))-(kpm*(y(1)*y(2))); ((kmp*y(1)+y(3))*y(2)*(1-y(2)))-y(2);(klm*(1+tanh(2*y(2)-2)))-(kl*y(3))];
tspan = [0, 25];
xinit = [0.01, 0.05, 0.539];
ode45(f, tspan, xinit)
legend('p(t)', 'm(t)', 'l(t)')

  2 Comments

Torsten
Torsten on 26 Feb 2018
And what was the problem with your first version from above ?
Best wishes
Torsten.
Daniel Jovin
Daniel Jovin on 26 Feb 2018
I'm new to Matlab, so I don't really understand what I did incorrectly and what differentiates my failed solution from the correct solution. But, the problem was that the plot I was generating, Figure 1, was incorrect- the values from the graph were not in the correct range and lacked the periodic nature of the graph from the modeling paper, Fig. 1.

Sign in to comment.

More Answers (2)

Rahiti Arapari Hargous
Rahiti Arapari Hargous on 31 Jul 2018
Hi everyone,
I have a similar problem with the next governing equation, with y(t) and x(t). The system of equation is the next one:
ay''' - b*x*y' + w*cos(y) = 0
x' - w*L*sin(y)=0
If anyone can help thanks

  0 Comments

Sign in to comment.


thanos tria
thanos tria on 26 Mar 2019
Edited: Walter Roberson on 26 Mar 2019
Hello, i am trying to solve a non linear system but when i plot the solutions only a straight line in one of the variables appears on the graph, and all the variables just can't get away from the initial conditions. That's my code:
syms c1(t) c2(t) c3(t) c4(t) T Y
Eqns = [diff(c1(t),t) == -0.015*c1+1.8*10^(-14)*c2*c4;
diff(c2(t),t) == 0.015*c1-1.8*10^(-14)*c2*c4;
diff(c3(t),t) == -6*10^(-34)*c3*0.051765*10^(20)*0.2465*10^(20)+0.015*c1;
diff(c4(t),t) == 6*10^(-34)*c3*0.051765*10^(20)*0.2465*10^(20)-1.8*10^(-14)*c2*c4];
[DEsys,Subs] = odeToVectorField(Eqns);
DEFcn = matlabFunction(DEsys, 'Vars',{T,Y});
tspan = [0,50000];
y0 = [0.040923*10^(12) 0 0 0];
[T,Y] = ode45(DEFcn, tspan, y0);
plot(T,Y)
And that's my graph:
It sould have had some curves because it describes how ozon is formulated in the atmosphere. Any suggestions??

  2 Comments

Walter Roberson
Walter Roberson on 26 Mar 2019
I agree that when I use those equations with Maple's numeric ODE solver, that I get distinct curves for each of the functions, for all of Maple's numeric methods that I have tried.
I also find in Maple that using stiff methods improves performance a lot without changing the solution. That would correspond to using routines such as ode23s() in MATLAB --- which produces the same (not so good) result as ode45 except faster.

Sign in to comment.