Solving a system of ODE

Question

Esraa Abdelkhaleq 2017-2-14

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/324976-solving-a-system-of-ode

评论： Star Strider 2017-2-15

采纳的回答： Star Strider

在 MATLAB Online 中打开

Hello,

How can I solve a system of two ODE's as follows,

dNcdt= Nc(t)*Kgr- dc*Nc(t)*Ni(t)
dNidt= ai*Nc(t) - di*Ni(t)

To obtain Nc(t) and Ni(t).

Thanks in advance.

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

请先登录，再进行评论。

请先登录，再回答此问题。

Answer 1

Star Strider 2017-2-14

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/324976-solving-a-system-of-ode#answer_254745

在 MATLAB Online 中打开

It is easier to let the Symbolic Math Toolbox do the algebra:

syms ai dc di Kgr Nc(t) Ni(t) t Y
Eq1 = diff(Nc) == Nc(t)*Kgr- dc*Nc(t)*Ni(t);
Eq2 = diff(Ni) == ai*Nc(t) - di*Ni(t);
[odesys, vrs] = odeToVectorField(Eq1, Eq2);
odefcn = matlabFunction(odesys, 'Vars',[t Y ai dc di Kgr])
odefcn =
    function_handle with value:
      @(t,Y,ai,dc,di,Kgr)[ai.*Y(2)-di.*Y(1);Kgr.*Y(2)-dc.*Y(1).*Y(2)]

The rest should be relatively straightforward for you to complete. To solve it numerically, begin with ode45, and if your equation turns out to be ‘stiff’ because of a wide variation in parameter magnitudes, and ode45 has problems, use ode15s or one of the other stiff solvers appropriate to your system.

8 个评论
显示 6更早的评论隐藏 6更早的评论

Star Strider 2017-2-14

在 MATLAB Online 中打开

My pleasure.

You have defined more parameters than in your original system of differential equations, so I will let you decide how best to pass them to odefcn in my initial Answer. If you have a new equation or set of equations, you will have to use the Symbolic Math Toolbox to create the vector field and MATLAB anonymous function to use in your numerical integration code.

You can easily solve it numerically with one of the ODE solvers. Define your parameters with their values, then the ODE function call would be:

odefcn = @(t,Y,ai,dc,di,Kgr) [ai.*Y(2)-di.*Y(1); Kgr.*Y(2)-dc.*Y(1).*Y(2)];
ai  = ...;
dc  = ...;
di  = ...;
Kgr = ...;
tspan = [0 10];
initcond = [Nc0  Ni0];
[t,y] = ode45(@(t,Y) odefcn(t,Y,ai,dc,di,Kgr), tspan, initcond);
figure(1)
plot(t, y)
grid
xlabel('Time')
ylabel('N')
legend('Ni(t)', 'Nc(t)')

Choose the appropriate values for ‘tspan’ (time span or time vector for the integration).

Try ode45 first. As I mentioned, you may need a stiff solver such as ode15s rather than ode45 if ode45 seems to ‘get stuck’ and take more than a few minutes to integrate your equations.

Esraa Abdelkhaleq 2017-2-15

OK. Thanks a lot for your help.

Star Strider 2017-2-15

在 MATLAB Online 中打开

My pleasure.

Out doing stuff for 3 hours (life intrudes) so just got back to this topic.

This code:

syms ai dc di fc fch fi fih Ke Kgr Nc(t) Ni(t) Rc Rch Ri Rih t Y Nc0 Nch0 Nih0 Ni0 qpl(t)
Eq1 = diff(Nc) == Nc(t)*Kgr- dc*Nc(t)*Ni(t);
Eq2 = diff(Ni) == ai*Nc(t) - di*Ni(t);
Eq3 = diff(qpl) == fc*Rc*Nc(t)+ fi*Ri*Ni(t)+fch*Rch*Nch0+ fih*Rih*Nih0-Ke*qpl(t) ;
[odesys, vrs] = odeToVectorField(Eq1, Eq2, Eq3)
odefcn = matlabFunction(odesys, 'Vars',[t Y ai dc di fc fch fi fih Ke Kgr Nc0 Nch0 Nih0 Ni0 Rc Rch Ri Rih])

produces:

odesys =
                                             ai*Y[2] - di*Y[1]
                                       Kgr*Y[2] - dc*Y[1]*Y[2]
Nch0*Rch*fch - Ke*Y[3] + Nih0*Rih*fih + Rc*fc*Y[2] + Ri*fi*Y[1]
vrs =
    Ni
    Nc
   qpl
odefcn = @(t,Y,ai,dc,di,fc,fch,fi,fih,Ke,Kgr,Nc0,Nch0,Nih0,Ni0,Rc,Rch,Ri,Rih) [ai.*Y(2)-di.*Y(1);Kgr.*Y(2)-dc.*Y(1).*Y(2);-Ke.*Y(3)+Nch0.*Rch.*fch+Nih0.*Rih.*fih+Rc.*fc.*Y(2)+Ri.*fi.*Y(1)];

The ‘vrs’ output you can interpret as:

Y(1) =  Ni
Y(2) =  Nc
Y(3) = qpl

You have to supply all the values for the parameters, so you can then integrate it with ode45 or ode15s (or other solver, depending on your requirements).

请先登录，再进行评论。