用Python求常微分方程的解析解

简介

sympy提供了常微分方程求解器，以下面的微分方程为例

$$9f\left(x\right)+\frac{d^2}{d{x}^2}f\left(x\right)=0$$

可调用dsolve来进行求解

from sympy import Function, dsolve, Derivative
from sympy.abc import x
from sympy import print_latex
f = Function('f')
ff = dsolve(Derivative(f(x), x, x) + 9*f(x), f(x))
print_latex(ff)

结果如下

$$f\left(x\right)=C_1\sin \left(3x\right)+C_2\cos \left(3x\right)$$

dsolve

dsolve的完整参数如下

sympy.solvers.ode.dsolve(eq, func=None, hint='default', simplify=True, ics=None, xi=None, eta=None, x0=0, n=6, **kwargs)

各参数含义为

• eq 即常微分方程表达式，是唯一必须输入从桉树
• func 一个单变量函数，即dsolve的求解对象，一般无需指明，dsolve会自动从eq中检测待求解函数。
• hint 准备采用的求解方法，默认采用classify_ode返回的方法。
• simplify 为True时，通过odesimp()函数进行简化
• ics 微分方程的边界条件
• xi, eta 是两个极小函数
• x0 幂级数展开点
• n 幂级数所对应的因变量的指数

hint

sympy.solvers.ode中提供了allhints元组，给出了所有hint的可选方法，共计45个，下面随机挑选两个进行测试。

# 可查看所有的hint
from sympy.solvers.ode import allhints

以下面的函数为例

$$\sin \left(x\right)\cos \left(f\left(x\right)\right)+\sin \left(f\left(x\right)\right)\cos \left(x\right)\frac{d}{dx}f\left(x\right)$$

from sympy import sin, cos
eq = sin(x)*cos(f(x)) + cos(x)*sin(f(x))*f(x).diff(x)
f1 = dsolve(eq, hint='1st_exact')
f2 = dsolve(eq, hint='almost_linear')

二者的求解结果完全一致：

$$\left[f\left(x\right)=?\mathrm{acos}\left(\frac{C_1}{\cos \left(x\right)}\right)+2\pi ,f\left(x\right)=\mathrm{acos}\left(\frac{C_1}{\cos \left(x\right)}\right)\right]$$

常微分方程组

dsolve也可以求解常微分方程组，以下面的方程组为例

$$\begin{gathered} x^{\prime }\left(t\right)=12tx\left(t\right)+8y\left(t\right) \\ {y}^{\prime }\left(t\right)=7ty\left(t\right)+21x\left(t\right) \end{gathered}$$

求解方法如下

from sympy import Eq, symbols
from sympy.abc import t

x, y = symbols('x, y', cls=Function)
eq = (Eq(Derivative(x(t),t), 12*t*x(t) + 8*y(t)), Eq(Derivative(y(t),t), 21*x(t) + 7*t*y(t)))
xys = dsolve(eq)
print_latex(xys)

$$\begin{gathered} x\left(t\right)=C_1{x}_0\left(t\right)+C_2{x}_0\left(t\right)\int \frac{8\left(e^{\int 7tdt}\right)e^{\int 12tdt}}{x_0^2\left(t\right)}dt, \\ y\left(t\right)=C_1{y}_0\left(t\right)+C_2\left(y_0\left(t\right)\int \frac{8\left(e^{\int 7tdt}\right)e^{\int 12tdt}}{x_0^2\left(t\right)}dt+\frac{\left(e^{\int 7tdt}\right)e^{\int 12tdt}}{x_0\left(t\right)}\right) \end{gathered}$$