Thread: System non-linear ODE's

Hi all,
I have a dynamical system of two non-linear autonomous ODE's.
The general form is

q'[t] = F_1(q,z)
z'[t] = F_2(q,z)

Actually, I'm interested in the function q[t], given some initial conditions q[0] and z[0]. See attachment for more details and the specific form of F_1 and F_2.

Is there some methodology that this can be answered? Also any literature or advice would be highly appreciated.

Thanks in advance,
idaios

Ok I looked at them. Why not solve them numerically in Mathematica? Just do an NDSolve on them. Tell you what, you clearly supply the values of the parameters and initial conditions for a particular case, and I'll run it and post the Mathematica code so that you can run it by changing the setup.

Originally Posted by shawsend

Ok I looked at them. Why not solve them numerically in Mathematica? Just do an NDSolve on them. Tell you what, you clearly supply the values of the parameters and initial conditions for a particular case, and I'll run it and post the Mathematica code so that you can run it by changing the setup.

Hi,
thanks! yes I can supply initial conditions:
q[0] = 0.001;
z[0] = -10
a1=0.5
a2=1
s1=0.0005
s2=0.001
omega=1
sigma^2 = 1

please, if you do it with mathematica could you upload the file as well?
Do you think that with mathematica is impossible to find a more general solution that will be a function of (arbitrary) parameter values and initial conditions?

best
idaios

Here's the deal with Mathematica: in order to cut and paste Mathematica code, I need to first convert it to Raw-Input form which I did and pasted it below but it's hard to read. So cut and paste it back into your Mathematica session, select it, then do a Cell/Convert To/Standard Form so that it's easier to read. I first defined the two equations as eqn1 and eqn2 so make sure I input it correctly. Then defined the parameters and supplied the initial conditions to NDSolve and ran it from t=0 to t=10, obtained the solutions and then did a plot of them which I posted below. Red is q and Blue is z. Also, I don't think you can obtain an explicit solution although Mathematica has DSolve for that. You can try running DSolve[{eqn1,eqn2},{q,z},t]. I tried but stopped it after about one minute as it did not come back with an answer.

Code:

eqn1 = Derivative[1][q][t] == 
    (1/(Exp[-(z[t]^2/(2*(\[Sigma]^2 + \[Omega]^2)))]*(-1 + q[t])^2 + 
       2*Exp[-((a1 + z[t])^2/(2*(\[Sigma]^2 + \[Omega]^2)))]*(-1 + s1)*(-1 + q[t])*
        q[t] - Exp[-((a2 + z[t])^2/(2*(\[Sigma]^2 + \[Omega]^2)))]*(-1 + s2)*
        q[t]^2))*(1 - q[t])*q[t]*(Exp[-(z[t]^2/(2*(\[Sigma]^2 + \[Omega]^2)))]*
       (-1 + q[t]) + Exp[-((a1 + z[t])^2/(2*(\[Sigma]^2 + \[Omega]^2)))]*(-1 + s1)*
       (-1 + q[t]) + Exp[-((a1 + z[t])^2/(2*(\[Sigma]^2 + \[Omega]^2)))]*(-1 + s1)*
       q[t] - Exp[-((a2 + z[t])^2/(2*(\[Sigma]^2 + \[Omega]^2)))]*(-1 + s2)*q[t]); 
eqn2 = Derivative[1][z][t] == (1/(\[Sigma]^2 + \[Omega]^2)^(3/2))*\[Omega]*
     ((-Exp[-(z[t]^2/(2*(\[Sigma]^2 + \[Omega]^2)))])*(-1 + q[t])^2*z[t] - 
      2*Exp[-((a1 + z[t])^2/(2*(\[Sigma]^2 + \[Omega]^2)))]*(-1 + s1)*(-1 + q[t])*
       q[t]*(a1 + z[t]) + Exp[-((a2 + z[t])^2/(2*(\[Sigma]^2 + \[Omega]^2)))]*
       (-1 + s2)*q[t]^2*(a2 + z[t])); 
a1 = 0.5; 
a2 = 1; 
s1 = 0.0005; 
s2 = 0.001; 
\[Omega] = 1; 
\[Sigma] = 1; 
mysol = NDSolve[{eqn1, eqn2, q[0] == 0.001, z[0] == -10}, {q, z}, 
   {t, 0, 10}]
Plot[Evaluate[{q[t], z[t]} /. mysol], {t, 0, 10}, 
  PlotStyle -> {Red, Blue}]

thanks a lot,
is it possible to solve it q[t] against t and not q against z?

thanks a lot for the mathematica code and advices

idaios

Yes:

myq[t_] := Evaluate[q[t] /. mysol];

Now you can treat myq[t] as the solution q(t) although may be a better way to do this. Also can take derivatives of the numerical results. For example, I can do the above code, take the derivative of it via:

mydq[t_]=D[myq[t],t]

Solve it for t=2, then check the right and left side of the first equation:

Code:

In[172]:=
myq[t_] := Evaluate[q[t] /. mysol]; 
myz[t_] := Evaluate[z[t] /. mysol]; 
mydq[t_] = D[q[t], t]; 
N[(1/(Exp[-(z[t]^2/(2*(\[Sigma]^2 + \[Omega]^2)))]*(-1 + q[t])^2 + 
      2*Exp[-((a1 + z[t])^2/(2*(\[Sigma]^2 + \[Omega]^2)))]*(-1 + s1)*(-1 + q[t])*
       q[t] - Exp[-((a2 + z[t])^2/(2*(\[Sigma]^2 + \[Omega]^2)))]*(-1 + s2)*
       q[t]^2))*(1 - q[t])*q[t]*(Exp[-(z[t]^2/(2*(\[Sigma]^2 + \[Omega]^2)))]*
      (-1 + q[t]) + Exp[-((a1 + z[t])^2/(2*(\[Sigma]^2 + \[Omega]^2)))]*(-1 + s1)*
      (-1 + q[t]) + Exp[-((a1 + z[t])^2/(2*(\[Sigma]^2 + \[Omega]^2)))]*(-1 + s1)*
      q[t] - Exp[-((a2 + z[t])^2/(2*(\[Sigma]^2 + \[Omega]^2)))]*(-1 + s2)*
      q[t])] //. {z -> myz, q -> myq, t -> 2}
N[mydq[2]]

Out[175]=
{0.25751817581793557}

Out[176]=
{0.2575182526588393}

Pretty close.