Solving a system of ODE

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.

Answer 1

[Star Strider]

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.

[Esraa Abdelkhaleq]

Thanks for your answer but I can not understand what Y means?

I really can not understand the code or how to complete the solution.

My problem is to solve the following ODE:

 dqpldt=fc*Rc*Nc(t)+ fi*Ri*Ni(t)+fch*Rch*Nch0+ fih*Rih*Nih0-Ke*qpl(t) ;

So I need to get Nc(t) and Ni(t) to substitute in the model equation in order to solve it analytically and get qpl(t).

[Star Strider]

My pleasure.

I probably should have been a bit more explicit. The ‘vrs’ output of the odeToVectorField call, indicates that ‘Y(1)=vrs(1)’ and ‘Y(2)=vrs(2)’, so ‘Y(1)=Ni’ and ‘Y(2)=Nc’.

An analytic solution may not be possible, since your equation is nonlinear, and only a few nonlinear equations have analytic solutions. When I use dsolve:

NiNc = dsolve(Eq1, Eq2, Nc(0) == Nc0, Ni(0) == Ni0)

the result is:

Warning: Explicit solution could not be found. 
> In dsolve (line 201)

NiNc =

[ empty sym ]

You could use other solvers, perhaps Wolfram Alpha (free online) or Maple, but I would be surprised if you could get an analytic solution. I could not get one using Wolfram Alpha with your system (link).

[Esraa Abdelkhaleq]

Thanks a lot.

What if I try to solve the model equation numerically or try to graph the solution?

I mean if I have the values of the parameters in the following equation:

dqpldt=fc*Rc*Nc(t)+ fi*Ri*Ni(t)+fch*Rch*Nch0+ fih*Rih*Nih0-Ke*qpl(t) ;

The parameters are:

fc; Rc; Nc0; fi*Ri; Ni0; fch*Rch*Nch0; fih*Rih*Nih0; Ke; ai; di; dc; Kgr

Can I solve the equation numerically/ or graph the solution while I do not have Nc(t) or Ni(t) but I only have dNc/dt and dNi/dt?

[Star Strider]

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]

The parameters are for the full model equation (dqpl/dt) in addition to the ODE system (dNc/dt, dNi/dt).

When I try to run the code you wrote in your comment, an error appears:

Error using odearguments (line 110)
Inputs must be floats, namely single or double.

Error in ode45 (line 113)
[neq, tspan, ntspan, next, t0, tfinal, tdir, y0, f0, odeArgs, odeFcn, ...

Error in Untitled5 (line 11)
[t,y] = ode45(@(t,Y) odefcn(t,Y,ai,dc,di,Kgr), tspan, initcond);

In my previous comment, I did not want to solve the ODE system numerically but the full model equation (dqpl/dt):

dqpldt=fc*Rc*Nc(t)+ fi*Ri*Ni(t)+fch*Rch*Nch0+ fih*Rih*Nih0-Ke*qpl(t) ;

But I do not have Nc(t) or Ni(t), I have instead dNc/dt and dNi/dt. So I wondered If I can write a code including the ODE system (dNc/dt and dNi/dt) together with the full model equation (dqpl/dt) and try to at least graph the solution for qpl(t).

I wish I could explain my problem.

Thanks for your time.

[Star Strider]

My pleasure.

I have no idea what could be throwing that error. Be sure that the parameters you are passing to ‘odefcn’ are numeric single- or double-precision values, and not symbolic.

Analytic solutions do not exist for ‘Nc(t)’ and ‘Ni(t)’, so you cannot get an analytic solution for ‘qpl(t)’ in your differential equation. You have to include it as the third differential equation in your system of differential equations, just as I did for the first two, then use the Symbolic Math Toolbox to create a function from it that the numeric ODE solvers can use. That should work.

[Esraa Abdelkhaleq]

OK. Thanks a lot for your help.

[Star Strider]

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).