【Python】求解常微分方程的符号解+数值解例题

🌾常微分方程的符号解

Python提供了功能强大的求解常微分方程符号解的函数dsolve。
下面我们来看几道例题。

🌼例1

求解微分方程$y^{{}^{\prime }}=-2y+2x^2+2x,y(0)=1$
代码如下：

import sympy as sp
x=sp.var('x')
y=sp.Function('y')
eq=y(x).diff(x)+2*y(x)-2*x*x-2*x
s=sp.dsolve(eq,ics={y(0):1})
s=sp.simplify(s)
print(s)

运行结果如下：

在交互式窗口中运行得到：

🌼例2

求解二阶微分方程
$y^{{}^{\prime \prime }}-2y^{{}^{\prime }}+y=e^x,y(0)=1,y^{\prime }(0)=-1$
代码如下：

import sympy as sp
x=sp.var('x')
y=sp.Function('y')
eq=y(x).diff(x,2)-2*y(x).diff(x)+y(x)-sp.exp(x)
con={y(0):1,y(x).diff(x).subs(x,0):-1}
s=sp.dsolve(eq,ics=con)
print(s)

解得方程为：

🌼例3

已知输入信号为$u(t)=e^{-t}cos(t)$,
试求下面微分方程的解。
$y^{(4)}(t)+10y^{(3)}(t)+35y^{{}^{\prime \prime }}(t)+50y^{{}^{\prime }}(t)+24y(t)=u^{{}^{\prime \prime }}(t)$
$y(0)=0,y^{{}^{\prime }}(0)=-1,y^{{}^{\prime \prime }}(0)=1,y^{{}^{\prime \prime \prime }}(0)=1$
代码如下：

import sympy as sp
t=sp.var('t')
y=sp.Function('y')
u=sp.exp(-t)*sp.cos(t)
eq=y(t).diff(t,4)+10*y(t).diff(t,3)+35*y(t).diff(t,2)+50*y(t).diff(t)+24*y(t)-u.diff(t,2)
con={y(0):0,y(t).diff(t).subs(t,0):-1,y(t).diff(t,2).subs(t,0):1,y(t).diff(t,3).subs(t,0):1}
s=sp.dsolve(eq,ics=con)
s=sp.expand(s)
print(s)

求得解如下：

🌾常微分方程的数值解

🌼例4

求解例1 微分方程$y^{{}^{\prime }}=-2y+2x^2+2x,y(0)=1$的数值解。
代码：

from scipy.integrate import odeint
import numpy as np
import pylab as plt
import sympy as sp

dy=lambda y,x:-2*y+2*x*x +2*x
xx=np.linspace(0,3,31)#自变量的取值
s=odeint(dy,1,xx)#y的取值
print('x={}\n对应的数值解y={}'.format(xx,s.flatten()))
plt.plot(xx,s,'*')

x=sp.var('x')
y=sp.Function('y')
eq=y(x).diff(x)+2*y(x)-2*x*x-2*x
s2=sp.dsolve(eq,ics={y(0):1})
sx=sp.lambdify(x,s2.args[1],'numpy')

plt.plot(xx,sx(xx)) #观察吻合度
plt.show()

运行结果：

🌼例5

求解微分方程$(1-x)\frac{d^{2}y}{d{x}^2}=\frac{1}{5}\sqrt{1+(\frac{dy}{dx}{)}^2}$ 0＜x≤1
的数值解。
求数值解时，需要把二阶微分方程转化为一阶微分方程，引进$y_1=y,y_2=y^{{}^{\prime }}$,则可以将方程转化为如下的一阶微分方程组：
$$y_1^{{}^{\prime }}=y_2,y_1(0)=0$$
$$y_2^{{}^{\prime }}=\frac{1}{5(1-x)}\sqrt{1+y_2^2},y_2(0)=0$$

代码：

from scipy.integrate import odeint
import numpy as np
import pylab as plt
yx=lambda y,x:[y[1],np.sqrt(1+y[1]**2)/5/(1-x)]
x0=np.arange(0,1,0.00001)
y0=odeint(yx,[0,0],x0)

plt.rc('font',size=16)
plt.plot(x0,y0[:,0])
plt.show()

运行结果：