使用 odeint 和 sympy 求解微分方程的问题

Question

我正在尝试求解并显示以下等式的图形：

f'=af²-bf

因此我尝试用scipy.integrate.odeint库函数解决，没有成功。

这是我到目前为止所做的：

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)

然后我想针对不同的 f₀（条件为 t=0）显示它，就像您对正规方程所做的那样看起来像这样：

Answer 1

我使用带有 isympy shell 的 Sympy 解决了你的问题（简化了一些定义名称等）——你必须更加小心地完成所有必要的 import s.

我解了微分方程，取了所得方程的 rhs（右侧）并将 rhs 表达式赋值给变量 sol.
我求解了 c1 方程 f(0)-x0=0。
我为涉及的不同符号分配了（任意）值。
我在用实际值代替所有变量后绘制了函数，t 是我们绘图的自由变量。

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]:

附录

为了完整性，数值解使用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。对于新代码，使用 scipy.integrate.solve_ivp 求解微分方程。}

Answer 2

了解 closed-form 解决方案（如果可以）总是有帮助的。我请 Wolfram Alpha 求解您的 ODE。 Here是他们给我的答案

我建议你把它画出来，这样你看到它就会知道正确的答案。

使用 odeint 和 sympy 求解微分方程的问题

Problem solving differential equations using odeint and sympy

python

math

sympy

odeint