SciPy求微分方程数值解

一阶微分方程(组)

scipy.integrate.odeint()可以数值解一阶微分方程(组)

不同于微分方程的符号解,可以给出未知函数的表达式,数值解只能求得指定位置的函数值

sol=odeint(func, y0, t)

  • func:表示一个微分方程\(\frac{\text dy}{\text dt} = func(y, t)\),其中的y可以是向量,以表示微分方程组
  • y0:未知函数的初值条件,表示\(y(t[0])=y_0\),即第一个自变量取值处的函数值,同样可以是向量
  • t:自变量的取值序列

返回值sol是一个len(t)*len(y0)的数组,表示每个取值位置的函数值

示例: \[ \begin{cases} y'=-2y+x^2+2x \\ y(1)=2 \end{cases} \] 代码在下面一起给出

高阶微分方程(组)

通过换元可以将高阶微分方程转化为一阶微分方程: \[ y_0=y \\ y_1=y_0'=y' \\ y_2=y_1'=y'' \\ \dots \\ y_k=y_{k-1}'=y^{(k)} \\ \] 这样就同样可以使用odeint()来获得数值解了

示例: \[ \begin{cases} y''+2y'+2y=0 \\ y(0)=0, y'(0)=1 \end{cases} \] 先转化成: \[ \begin{cases} y_0'=y_1 \\ y_1'=-2y_1-2y_0 \\ y_0(0)=0, y_1(0)=1 \end{cases} \]

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
import numpy as np
from scipy.integrate import odeint
import pylab as plt


def dy(y, x):
return -2*y+x**2+2*x


x = np.arange(1, 10.5, 0.5)
sol = odeint(dy, 2, x)
print(sol.T)


def dy(y, x):
y0, y1 = y
return [y1, -2*y1 - 2*y0]


x = np.arange(0, 10, 0.1)
sol = odeint(dy, [0, 1], x)
plt.scatter(x, sol.T[0], marker='.') # 只需要y的0阶导数据
plt.savefig('odeint.png', dpi=500)
plt.show()

输出:

1
2
3
4
[[ 2.          2.08484933  2.9191691   4.18723381  5.77289452  7.63342241
9.75309843 12.12613985 14.75041934 17.62515427 20.75005673 24.12502089
27.7500077 31.62500278 35.75000104 40.1250004 44.75000015 49.62500006
54.75000002]]

Lorenz模型的混沌效应

Lorenz方程: \[ \begin{cases} \dot x=\sigma(y-x) \\ \dot y=\rho x-y-xz \\ \dot z=xy-\beta z \end{cases} \] 模拟\(\sigma=10,\rho=28,\beta=\frac 8 3\)\(t\in [0,50]\)的运行轨迹,以及给初值微小扰动后解的偏差演化曲线

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
import numpy as np
from scipy.integrate import odeint
import pylab as plt


def lorenz(w, t):
sigma = 10
rho = 28
beta = 8/3
x, y, z = w
return [sigma*(y-x), rho*x-y-x*z, x*y-beta*z]


t = np.arange(0, 50, 0.01)
sol1 = odeint(lorenz, [0., 1., 0.], t)
sol2 = odeint(lorenz, [0., 1.0001, 0.], t)

ax1 = plt.subplot(121, projection='3d')
ax1.plot(sol1[:, 0], sol1[:, 1], sol1[:, 2], linewidth=1)
ax2 = plt.subplot(122, projection='3d')
ax2.plot(sol1[:, 0]-sol2[:, 0], sol1[:, 1] -
sol2[:, 1], sol1[:, 2]-sol2[:, 2], linewidth=1)
plt.savefig('lorenz.png', dpi=500)
plt.show()