Problem solving differential equations using odeint and sympy

Question

I am trying to solve and display a graph of the following equation:

f'=af²-bf

Therefore I have tried to use scipy.integrate.odeint library function to solve it without success.

Here is what I have done so far:

from IPython.display import display
import sympy as sy
from sympy.solvers.ode import dsolve
import matplotlib.pyplot as plt
import numpy as np
from sympy import functions
from sympy import Function, Symbol

sy.init_printing()

t = sy.symbols("t", real=True)
f = sy.symbols("f", cls=functions)

eq1 = sy.Eq(f(t).diff(t), (0.1*f(t)**2 - 2*f(t)))
sls = dsolve(eq1)

print("For ode")
display(eq1)
print("the solutions are:")
for s in sls:
    display(s)

# plot solutions:
x = np.linspace(0, 3, 100)
fg, axx = plt.subplots(5000, 3)

Then I would like to display it for different f₀ (condition at t=0) as you would do for a normal equation which would look like this:

Answer 1

It always helps to know the closed-form solution if you can get it. I asked Wolfram Alpha to solve your ODE. Here is the answer they gave me.

I'd recommend plotting that so you'll know the right answer when you see it.

Answer 2

I solved your problem using Sympy with the isympy shell (that ease a little defining names, etc) — you have to be more careful to do all the necessary imports.

1. I solved the differential equation, took the rhs (right hand side) of the resulting equation and assigned said rhs expression to the variable sol.
2. I solved for c1 the equation f(0)-x0=0.
3. I assigned (arbitrarily) values to the different symbols involved.
4. I plotted the function after substituting actual values for all the variables except for t, the free variable of our plot.

18:43 boffi@debian:~ $ isympy 
IPython console for SymPy 1.4 (Python 3.7.4-64-bit) (ground types: gmpy)

These commands were executed:
>>> from __future__ import division
>>> from sympy import *
>>> x, y, z, t = symbols('x y z t')
>>> k, m, n = symbols('k m n', integer=True)
>>> f, g, h = symbols('f g h', cls=Function)
>>> init_printing()

Documentation can be found at https://docs.sympy.org/1.4/


In [1]: a, b = symbols('a b')                                                             

In [2]: sol = dsolve(Derivative(f(t), t) + a*f(t)**2 + b*f(t), f(t)).rhs                  

In [3]: c1, x0 = symbols('C1 x0')                                                         

In [4]: c1_from_x0, = solve(sol.subs(t, 0) - x0, c1)                                      

In [5]: sol, c1_from_x0, sol.subs(c1, c1_from_x0)                                         
Out[5]: 
⎛                       ⎛  a⋅x₀  ⎞                                ⎞
⎜        C₁⋅b        log⎜────────⎟                                ⎟
⎜     b⋅ℯ               ⎝a⋅x₀ + b⎠               b⋅x₀             ⎟
⎜──────────────────, ─────────────,  ──────────────────────────────⎟
⎜  ⎛   C₁⋅b    b⋅t⎞        b                   ⎛    a⋅x₀      b⋅t⎞⎟
⎜a⋅⎝- ℯ     + ℯ   ⎠                 (a⋅x₀ + b)⋅⎜- ──────── + ℯ   ⎟⎟
⎝                                              ⎝  a⋅x₀ + b       ⎠⎠

In [6]: values = {x0:10, a:3, b:4, c1:c1_from_x0}                                         

In [7]: plot(sol.subs(values), (t, 0, 0.5));

In [8]: sol.subs(values).simplify()                                                      
Out[8]: 
     20     
────────────
    4⋅t     
17⋅ℯ    - 15

In [9]:                                                                                  

---

Addendum

For completeness, the numerical solution using scipy.integrate.odeint ¹

from numpy import exp, linspace 
from scipy.integrate import odeint 
import matplotlib.pyplot as plt 
 
t = linspace(0, 0.5, 201) ; f0 = 10; a = 3 ; b = 4                                
f = odeint(lambda f, t: (-a*f -b)*f, f0, t)                                       

fig0, ax0 = plt.subplots(figsize=(8,3))                                           
ax0.plot(t, f) ; ax0.set_title('Numerical Solution')
plt.show()

exact = 20 / (17*exp(4*t)-15)                                                     
fig1, ax1 = plt.subplots(figsize=(8,3))                                           
ax1.plot(t, (f.flat-exact)*1E6) ; ax1.set_title('(Numerical-Analytical)*10**6')   
plt.show()

^{1. For new code, use scipy.integrate.solve_ivp to solve a differential equation.}