高等数学问题的数值解（6）--解微分方程

更新于
求微分方程的符号解，其关键词为odeint。

Python对常微分方程的数值求解是基于一阶方程进行的，高阶微分方程必须化成一阶方程组，通常采用龙格-库塔方法。
scipy.integrate模块的odeint函数可求常微分方程的数值解。

例：求下述微分方程
$ y^{\prime }=-2y+x^{2}+2x，y(1) = 2$
在$1\leq x \leq 10$内步长间隔为0.5点上的数值解

1
2
3
4
5
6
7
8
from scipy.integrate import odeint
from numpy import arange

dy = lambda y, x: -2*y+x**2+2*x
x = arange(1, 10.5, 0.5)
sol = odeint(dy, 2, x)

print('x={}\n对应的数值解y={}'.format(x, sol.T))

例：求下述微分方程的数值解，并在同一个图形界面上画出符号解和数值解的曲线
$\frac{d^{2}y}{dx^{2}}+2\frac{dy}{dx}+2y=sinx,y(0)=0,y^{\prime }(0)=1$
解：引入$y_{1}=y，y_{2}=y^{\prime }$，则可以把原来的二阶微分方程化为如下一阶微分方程组

$$
\begin{cases}
y_{1}^{\prime } = y_{2} & y_{1}(0)=0 \\
y_{2}^{\prime } = -2y_{1}-2y_{2}&y_{2}(0)=0
\end{cases}
$$

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
from scipy.integrate import odeint
from sympy.abc import t
import numpy as np
import matplotlib.pyplot as plt

def Pfun(y,x):
    y1, y2=y
    return np.array([y2, -2*y1-2*y2])

x = np.arange(0, 10, 0.1)  # 创建时间点
sol1 = odeint(Pfun, [0.0, 1.0], x)  # 求数值解

plt.rc('font',size=16); plt.rc('font',family='SimHei')
plt.plot(x, sol1[:,0],'r*',label="数值解")
plt.plot(x, np.exp(-x)*np.sin(x), 'g', label="符号解曲线")
plt.legend()
# plt.savefig("figure.png")
plt.show()
0%