Sympy入门之求解微分方程-优快云博客

本文链接：https://blog.youkuaiyun.com/weixin_73953650/article/details/146436140

Sympy是一个专注于符号数学计算的数学工具，使得用户可以轻松地进行复杂的符号运算，如求解方程、求导数、积分、级数展开、矩阵运算等。本文，我们将详细讲解Sympy在求解微分方程中的应用。

获取方式

pip install -i https://mirrors.tuna.tsinghua.edu.cn/pypi/web/simple sympy

注:本文所使用python环境为jupyter notebook,所有代码只有在jupyter ipynb环境中才可以正确输出

微分方程定义

$f(x,y,y',y''..y^{(n)})=0$

微分方程，是指含有未知函数及其各阶导数的关系式。简而言之,解微分方程就是找出未知函数。未知函数是一元函数的，叫做常微分方程；未知函数是多元函数的叫做偏微分方程。

实际上,只有少数简单的微分方程才可以求得解析解。不过即使没有找到其解析解，仍然可以确认其解的部分性质。在无法求得解析解时，可以利用数值分析的方式，来找到其数值解。

微分方程求解

微分方程在求解时可以分为ODE(常微分方程)与PDE(偏微分方程)两类,但无论是哪一类,只要其有显著的解析解,我们都可以使用Sympy进行求解。

Sympy中微分方程的表示

Sympy中微分方程的表示需要按照以下的格式

import sympy as sp
# 定义符号变量
'''
sp.Function()用来声明未知函数
比如常微分方程中带求解函数y

sp.Eq(表达式1,表达式2)
用来表示表达式1=表达式2
我们经常将微分方程表示为f(x,y,y')=0的形式
因此便可以所使用sp.Eq(f'(x,y,y'),0)来将其使用
sympy表示

Function对象下的diff(x,n)方法表示该未知函数对x的n阶导数

'''
x=sp.symbols('x')#自变量x
y=sp.Function('y')#待求解函数y
#y''(x)-y=0
diff_eq=sp.Eq(y(x).diff(x,2)-y(x),0)

其中，sp.Eq()函数用来表示一个等式。已经声明过的Function对象下的diff(x,n)方法表示该函数对x的n阶导数。

微分方程的求解

求解微分方程时，我们需要用到Sympy下的dsolve()函数，他可以用来求解一些简单，存在显著解析解的常微分与偏微分方程。

sp.dsolve()函数参数详解:

参数名	描述
`eq`	需要求解的微分方程。可以是 `Eq` 对象或表达式。
`func`	需要求解的未知函数。
`hint`	指定求解方法（可选）。默认情况下，`dsolve` 会自动选择合适的方法。
`ics`	初始条件（可选）。以字典形式提供，键为函数及其导数的条件，值为对应的值。
`simplify`	是否对结果进行简化（可选）。默认值为 `True`。
`**kwargs`	其他传递给求解方法的参数（可选）。

示例1：

$y''-2y'+y=0, y(0)=1,y'(0)=0$

import sympy as sp
# 定义符号变量
'''
sp.Function()用来声明未知函数
比如常微分方程中带求解函数y

sp.Eq(表达式1,表达式2)
用来表示表达式1=表达式2
我们经常将微分方程表示为f(x,y,y')=0的形式
因此便可以所使用sp.Eq(f'(x,y,y'),0)来将其使用
sympy表示

Function对象下的diff(x,n)方法表示该未知函数对x的n阶导数

'''
x=sp.symbols('x')#自变量x
y=sp.Function('y')#待求解函数y
#y''(x)-y=0
diff_eq=sp.Eq(y(x).diff(x,2)-2*y(x).diff(x,1)+y(x), 0)
# 求解微分方程
# 求解带初始条件的微分方程
solution=sp.dsolve(diff_eq,y(x),ics={y(0):1,y(x).diff(x,1).subs(x,0):0})
display(solution)

结果:

示例2：稍微复杂一点点的一阶微分方程 $(x^2y+xy)dy+(x^3-y^2)dx=0$

这是一个典型的代换方法求解的微分方程，这里我先给出我的计算步骤：

解: $\frac{dy}{dx}=\frac{y^2-x^3}{x^2y+xy}$ （1）

(1)式右侧分子分母同除 $xy$ 有 $\frac{dy}{dx}=\frac{\frac{y}{x}-\frac{x^2}{y}}{x+1}$ （2）

令 $\frac{y}{x}=u$ ，则 $y=xu$ , $\frac{dy}{dx}=x\frac{du}{dx}+u$ （3）

(3)式替换(2)式有 $\frac{du}{dx}=\frac{-\frac{u^2+1}{u}}{x+1}$

整理一下为: $\frac{-2}{x+1}dx=\frac{2u}{1+u^2}du$

两边各自积分有:

$C_{1}-ln(1+x)^2=ln(1+u^2)$

这里我们不妨将C1看做 $ln(e^{C_{2}})$ ,那么原等式为

$ln(\frac{e^{C_{2}}}{(1+x)^2})=ln(u^2+1)$

这里我们再令 $C=e^{C_{2}}$ ,那么便有

$ln(\frac{C}{(1+x)^2})=ln(u^2+1)$

这么做的目的就是为了把函数法则外的常数放到函数里边去这样在我们解方程时,便于化简。

考虑到 $ln(x)$ 为单调递增函数，且无对称轴,那么 $f(x_1)=f(x_2)$

必然有 $x_{1}=x_{2}$ ，即 $\frac{C}{(x+1)^2}=1+u^2$

即 $u=\frac{\sqrt{C-1-x^2-2x}}{x+1}$ ,这里我们不妨将 $C-1$ 看做整体并将 $u=\frac{y}{x}$ 带入化简，即可得到 $y=\pm\frac{x\sqrt{C-x^2-2x}}{x+1}$

Sympy验证

#sp.dsolve()求解复杂一些的常微分方程
import sympy as sp
from IPython.display import display,Math
# 定义符号变量
'''
sp.Function()用来声明未知函数
比如常微分方程中带求解函数y

sp.Eq(表达式1,表达式2)
用来表示表达式1=表达式2
我们经常将微分方程表示为f(x,y,y')=0的形式
因此便可以所使用sp.Eq(f'(x,y,y'),0)来将其使用
sympy表示

Function对象下的diff(x,n)方法表示该未知函数对x的n阶导数

'''
x=sp.symbols('x')#自变量x
y=sp.Function('y')#待求解函数y
diff_eq=sp.Eq(y(x).diff(x,1), (y(x)**2-x**3)/(x**2*y(x)+x*y(x)))
solution=sp.dsolve(diff_eq,y(x))
Latex_result=sp.latex(solution)
display(Math(Latex_result))

结果：