Ashish
Ashish

Reputation: 739

How to solve pair of nonlinear differential equations using Sympy

I have a system of two coupled nonlinear differential equations

for which I wrote the following code to get values of two functions using Sympy:

from sympy import *
t = symbols('t')
c, b, B, alp, mu = symbols('c b B alp mu', integer=True)
f, g = symbols('f g', cls=Function)
print(solve([Eq(f(t).diff(t)+c*f(t)+0.5*B*f(t)**2-b*g(t), 0), Eq(g(t).diff(t)+b*g(t)-c*f(t), (1-alp)*mu)], [f(t),g(t)]))

However, the results what I get is still in the form of derivatives, which cannot be used for my purpose. I also tried Mathematica for the same but it keeps on running indefinitely with no result. enter image description here Can someone please suggest a solution to such a system using Sympy or Mathematica?

Upvotes: 0

Views: 1201

Answers (2)

Nabil
Nabil

Reputation: 26

You can use ParametricNDSolve to allow any number of parameters to be decided later. Then use Manipulate to explore. Here's sample code:

{fx, gx} = Block[ {f, g, c, \[Beta], b, \[Alpha], \[Mu], f0, g0, tmax},
   {f, g} /. ParametricNDSolve[ {
      f'[t] + c f[t] + 0.5 \[Beta] f[t]^2 - b g[t] == 0,
      g'[t] + b g[t] - c f[t] - (1 - \[Alpha]) \[Mu] == 0,
      f[0] == f0, g[0] == g0
      }, {f, g}, {t, tmax}, {c, \[Beta], b, \[Alpha], \[Mu], f0, g0, 
      tmax} ]
   ];
Manipulate[
 ParametricPlot[ {fx[c, \[Beta], b, \[Alpha], \[Mu], f0, g0, 2^
     ln2tmax][t], 
   gx[c, \[Beta], b, \[Alpha], \[Mu], f0, g0, 2^ln2tmax][t]},
  {t, 0, 2^ln2tmax} ],
 {{c, 1}, -1, 1, 0.05},
 {{\[Beta], 1}, -1, 1, 0.05},
 {{b, 1}, -1, 1, 0.05},
 {{\[Alpha], 1}, -1, 1, 0.05},
 {{\[Mu], 1}, -1, 1, 0.05},
 Delimiter,
 {{f0, 1}, -1, 1, 0.05},
 {{g0, 1}, -1, 1, 0.05},
 Delimiter,
 {{ln2tmax, 0}, -3, 10, 0.5}
 ]

Upvotes: 0

zhk
zhk

Reputation: 331

You can solve your system very easily in Mathematica with NDSolve providing you with a numerical solution. But for this you will need to specify the numerical values for the parameters and two initial conditions. Here I choose random ones.

  ClearAll["Gobal"]
Eq1 = f'[t] + c1*f[t] + 0.5*B1*f[t]^2 - b1*g[t] == 0;
Eq2 = g'[t] + b1*g[t] - c1*f[t] == (1 - alp)*mu;
ibcs = {f[0] == 1, g[0] == 0};
sys = Join[{Eq1, Eq2}, ibcs];
B1 = 1; b1 = 1; c1 = 1; alp = 1; mu = 1;
sol1 = NDSolve[sys, {f[t], g[t]}, {t, 0, 10}];
Plot[{f[t], g[t]} /. sol1, {t, 0, 10}]

enter image description here

Upvotes: 0

Related Questions