How to solve two differential equations using ode45.

87 次查看(过去 30 天)
My system is this
x"+ x'+ x + y'=0;
y"+ y'+ y + x'=0;
i need to solve these differential equations using ode's.
thanks in advance.
  2 个评论
Torsten
Torsten 2018-2-22
You will need four initial conditions (for x,x',y,y') if you want to use ODE45 as numerical solver.
Ebraheem Menda
Ebraheem Menda 2018-2-22
Thank you Torsten. i have the initial conditions. but my question is how to convey these equations to ode45 or any other solver. Because they are coupled equations. thanks for your help.

请先登录,再进行评论。

采纳的回答

Torsten
Torsten 2018-2-22
fun=@(t,y)[y(2);-y(2)-y(1)-y(4);y(4);-y(4)-y(3)-y(2)];
x0=...; %supply x(t0)
x0prime=...; %supply x'(t0)
y0=...; %supply y(t0)
y0prime=...; %supply y'(t0)
Y0=[x0; x0prime; y0; y0prime];
tspan=[0 5];
[T,Y]=ode45(fun,tspan,Y0);
plot(T,Y(:,1)) %plot x
plot(T,Y(:,3)) %plot y
Best wishes
Torsten.
  9 个评论
Andreas Zimmermann
Andreas Zimmermann 2020-4-16
编辑:Andreas Zimmermann 2020-4-16
Hi Joe,
Using Torsten's relations
x = y(1)
x' = y(2)
y = y(3)
y' = y(4)
the first line
fun=@(t,y)[y(2); -y(2)-y(1)-y(4); y(4); -y(4)-y(3)-y(2)];
can be translated to the following, using Torsten's relations:
[x'; -x'-x-y'; y'; -y'-y-x']
which is essentially
[x'; x"; y'; y"].
He got these relations by solving for x" and y":
x"+ x'+ x + y'=0; => x" = -x' -x -y' = -y(2) -y(1) -y(4)
y"+ y'+ y + x'=0; => y" = -y' -y +x' = -y(4) -y(3) -y(2)
I have tried to use ode45 to solve the SIR model from this fantastic Geogebra Tutorial: (https://youtu.be/k6nLfCbAzgo)
% Just some start parameters and coefficients
Istart = 0.01;
Sstart = 0.99;
Rstart = 0;
transm = 3.2;
recov = 0.23;
maxT = 20;
% define S',I',R'
% S' = - transm * S * I
% I' = transm * S I - recov * I
% R' = recov * I
% S is y(1)
% I is y(2)
% R is y(3)
%so here fun = @(t,y)[S'; I'; R'];
fun = @(t,y)[-transm*y(1)*y(2); (transm*y(1)*y(2))-(recov*y(2)); recov*y(2)];
% Provide the starting parameters
Y0 = [Sstart; Istart; Rstart;];
% Define the range of t
tspan = [0 maxT];
% Magic happens and matrix Y contains S,I,R
[T,Y] = ode45(fun,tspan,Y0);
% Plot plot
figure
plot(T,abs(Y(:,1)),'b-') % S
hold on
plot(T,abs(Y(:,2)),'r-') % I
plot(T,abs(Y(:,3)),'g-') % R
xlim([0 23])
%Hope this helps :)
% Many thanks to Torsten and Ebraheem for your very very useful discussion!

请先登录,再进行评论。

更多回答(3 个)

Abraham Boayue
Abraham Boayue 2018-2-24
Hey Ebraheem There are many excellent methods that you can use to solve your problem, for instance, the finite difference method is a very powerful method to use. I can try with that.The ode45 function is a matlab built in function and was designed to solve certain ode problems, it may not be suitable for a number of problems. What are the initial values of your equations? Do you have any plot of the solution that one can use as a guide?
  1 个评论
Ebraheem Menda
Ebraheem Menda 2018-2-24
thank you Abraham for your response. i am yet to solve those equations and obtain initial conditions. once i got them i will post the equations. Meanwhile what is this finite difference method ? is it available in matlab 2009b ?

请先登录,再进行评论。


Abraham Boayue
Abraham Boayue 2018-2-24
The finite difference method is used to solve differential and partial equations. It is easier to implement in matlab. You can do the coding in any version of matlab, I have taken a course in numerical mathematics before and have a fairly good knowledge of how to solve such problems.
  5 个评论

请先登录,再进行评论。


Abraham Boayue
Abraham Boayue 2020-7-19
Hi Ebraheem Is there anything specific that you want me to do for you?
  2 个评论
Ebraheem Menda
Ebraheem Menda 2021-7-1
fun=@(t,y)[y(2);-y(2)-y(1)-y(4);y(4);-y(4)-y(3)-y(2)];
x0=...; %supply x(t0)
x0prime=...; %supply x'(t0)
y0=...; %supply y(t0)
y0prime=...; %supply y'(t0)
Y0=[x0; x0prime; y0; y0prime];
tspan=[0 5];
[T,Y]=ode45(fun,tspan,Y0);
plot(T,Y(:,1)) %plot x
plot(T,Y(:,3)) %plot y
hi Abraham i want to put this code in for loop and run it for 2 or more tspan. Could you please help out on this ?

请先登录,再进行评论。

类别

Help CenterFile Exchange 中查找有关 Programming 的更多信息

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by