数学建模：微分方程模型— Python 求解

使用 Python 求常微分方程的数值求解通常是基于一阶方程进行的，高阶微分方程要化成一阶方程组。

例：使用显式欧拉法和四阶龙格库塔法计算Lorenz模型

Lorenz模型是由美国气象学家Lorenz在研究大气运动时，通过简化对流模型，只保留3个变量提出的一个完全确定性的一阶自治常微分方程组（不显含时间变量），其方程为

$$\{\begin{array}{l}\dot{x}=\sigma (y-x),\\ \dot{y}=\rho x-y-xz,\\ \dot{z}=xy-\beta z.\end{array}$$

其中，参数 $\sigma$ 为 Prandtl 数，$\rho$ 为 Rayleigh 数，$\beta$ 为方向比。

这里分别用显式欧拉法和四阶龙格库塔法给出 Lorenz 模型在 $\sigma =10,\rho =28,\beta =\frac{8}{3}$ 时系统的三维演化轨迹，和系统从两个靠的很近的初值出发（相差仅 $0.0001$）后，解的偏差演化曲线。显示欧拉法及四阶龙格库塔法详见此处。

#显示欧拉法解lorenz方程
import numpy as np
import matplotlib.pyplot as plt

#显式欧拉法
def Euler(f,x0,y0,h,l):
    #f函数, x0,y0初值, h步长, l数量
    xn, yn = x0, y0
    n = 0
    ys = [y0]
    while n <= l:
        n += 1
        xn += h
        yn = yn + f(xn,yn) * h
        ys.append(yn)
    return np.array(ys)

#定义函数
def lorenz(t,w):
    sigma=10; rho=28; beta=8/3
    x, y, z=w
    return np.array([sigma*(y-x), rho*x-y-x*z, x*y-beta*z])
#t=np.arange(0, 50, 0.01)  #创建时间点
sol1=Euler(lorenz, 0, [0.0, 1.0, 0.0], 0.01, 5000)  #第一个初值问题求解
sol2=Euler(lorenz, 0, [0.0, 1.0001, 0.0], 0.01, 5000)  #第二个初值问题求解

#绘图
plt.rc('font',size=16); plt.rc('text',usetex=True)
#图A
ax1=plt.subplot(121,projection='3d')
ax1.plot(sol1[:,0], sol1[:,1], sol1[:,2],'r')
ax1.set_xlabel('$x$'); ax1.set_ylabel('$y$'); ax1.set_zlabel('$z$')
#图B
ax2=plt.subplot(122,projection='3d')
ax2.plot(sol1[:,0]-sol2[:,0], sol1[:,1]-sol2[:,1], sol1[:,2]-sol2[:,2],'g')
ax2.set_xlabel('$x$'); ax2.set_ylabel('$y$'); ax2.set_zlabel('$z$')
plt.show()
#print("sol1=",sol1, '\n\n', "sol1-sol2=", sol1-sol2)

import numpy as np
import matplotlib.pyplot as plt

#四阶龙格库塔
def RK4(f,x0,y0,h,maxx):
    #f函数, x0,y0初值, h步长, maxx: x最大值
    xn, yn = x0, y0
    n = 0
    ys = [yn]
    while xn < maxx:
        n += 1
        xn += h
        k1 = f(xn,yn)
        k2 = f(xn + (h / 2),yn + (h * k1) / 2)
        k3 = f(xn + (h / 2),yn + (h * k2) / 2)
        k4 = f(xn + h,yn + h * k3)
        yn = yn + (h / 6) * (k1 + 2 * k2 + 2 * k3 + k4)
        ys.append(yn)
    return np.array(ys)

#定义函数
def lorenz(t,w):
    sigma=10; rho=28; beta=8/3
    x, y, z=w
    return np.array([sigma*(y-x), rho*x-y-x*z, x*y-beta*z])
sol1=RK4(lorenz, 0, [0.0, 1.0, 0.0], 0.01, 50)  #第一个初值问题求解
sol2=RK4(lorenz, 0, [0.0, 1.0001, 0.0], 0.01, 50)  #第二个初值问题求解

#绘图
plt.rc('font',size=16); plt.rc('text',usetex=True)
#图A
ax1=plt.subplot(121,projection='3d')
ax1.plot(sol1[:,0], sol1[:,1], sol1[:,2],'r')
ax1.set_xlabel('$x$'); ax1.set_ylabel('$y$'); ax1.set_zlabel('$z$')
#图B
ax2=plt.subplot(122,projection='3d')
ax2.plot(sol1[:,0]-sol2[:,0], sol1[:,1]-sol2[:,1], sol1[:,2]-sol2[:,2],'g')
ax2.set_xlabel('$x$'); ax2.set_ylabel('$y$'); ax2.set_zlabel('$z$')
plt.show()
#print("sol1=",sol1, '\n\n', "sol1-sol2=", sol1-sol2)

两个程序的输出图像中，左图显示出我们经常听到的“蝴蝶效应”。右图给出了系统从两个靠的很近的初值出发（相差仅0.0001）后，解的偏差演化曲线，从图中看出随着时间的增大，两个解的差异越来越大，这正是动力系统对初值敏感性的直观表现。

scipy.integrate.odeint 求解微分方程模型

scipy.integrate 模块的 odeint 函数求常微分方程的数值解的基本调用格式为：

from scipy.integrate import odeint
sol=odeint(func, y0, t)

其中func是定义微分方程的函数或匿名函数，y0是初始条件的序列，t是一个自变量取值的序列（t的第一个元素一定为初始时刻），返回值sol是对应于序列t中元素的数值解，如果微分方程组中有$n$个函数，返回值sol是$n$列的矩阵，第$i(i=1,2,\cdots ,n)$列对应于第$i$个函数的数值解。

示例程序: 用 scipy.integrate.odeint 求解 Lorenz 方程

from scipy.integrate import odeint
import numpy as np
import matplotlib.pyplot as plt
#定义函数
def lorenz(w,t):
    sigma=10; rho=28; beta=8/3
    x, y, z=w
    return np.array([sigma*(y-x), rho*x-y-x*z, x*y-beta*z])
t=np.arange(0, 50, 0.01)  #创建时间点
sol1=odeint(lorenz, [0.0, 1.0, 0.0], t)  #第一个初值问题求解
sol2=odeint(lorenz, [0.0, 1.0001, 0.0], t)  #第二个初值问题求解
#绘图
plt.rc('font',size=16); plt.rc('text',usetex=True)
#图A
ax1=plt.subplot(121,projection='3d')
ax1.plot(sol1[:,0], sol1[:,1], sol1[:,2],'r')
ax1.set_xlabel('$x$'); ax1.set_ylabel('$y$'); ax1.set_zlabel('$z$')
#图B
ax2=plt.subplot(122,projection='3d')
ax2.plot(sol1[:,0]-sol2[:,0], sol1[:,1]-sol2[:,1], sol1[:,2]-sol2[:,2],'g')
ax2.set_xlabel('$x$'); ax2.set_ylabel('$y$'); ax2.set_zlabel('$z$')
plt.show()
#print("sol1=",sol1, '\n\n', "sol1-sol2=", sol1-sol2)

scipy.integrate.solve_ivp 求解微分方程模型

求解 $$\{\begin{array}{l}\frac{\mathrm{d}y}{\mathrm{d}t}=f(t,y),\\ y(t0)=y0.\end{array}$$

调用方式:

from scipy.integrate import solve_ivp
slove=scipy.integrate.solve_ivp(fun, t_span, y0, method='RK45', t_eval=None, 
		dense_output=False, events=None, vectorized=False, args=None, **options)

参数:

1. fun: 定义微分方程的函数, fun(t,y)
2. t_span: 2 元组的浮点数, 积分区间 (t0, tf)
3. y0: 数组, 形状 (n,), 初始状态
4. method: 使用的积分方法, 如Runge-Kutta方法 (‘RK23’、‘RK45’、‘DOP853’)
5. t_eval： 数组 或 None, 可选, 求解的时间, 必须在 t_span 内, 如果无(默认)，则使用求解器选择的点

返回值:

print(slove.t) #时间
print{slove.y) #解

示例程序: 用 scipy.integrate.odeint 求解 Lorenz 方程

from scipy.integrate import solve_ivp
import numpy as np
import matplotlib.pyplot as plt
#定义函数
def lorenz(t,w):
    sigma=10; rho=28; beta=8/3
    x, y, z=w
    return np.array([sigma*(y-x), rho*x-y-x*z, x*y-beta*z])
t=np.arange(0, 50, 0.01)  #创建时间点
sol1=solve_ivp(lorenz, (0,50), np.array([0.0, 1.0, 0.0]), 'RK45', t)  #第一个初值问题求解
sol2=solve_ivp(lorenz, (0,50), np.array([0.0, 1.0001, 0.0]), 'RK45', t)  #第二个初值问题求解
#绘图
plt.rc('font',size=16); plt.rc('text',usetex=True)
#图A
ax1=plt.subplot(121,projection='3d')
ax1.plot(sol1.y[0,:], sol1.y[1,:], sol1.y[2,:],'r')
ax1.set_xlabel('$x$'); ax1.set_ylabel('$y$'); ax1.set_zlabel('$z$')
#图B
ax2=plt.subplot(122,projection='3d')
ax2.plot(sol1.y[0,:]-sol2.y[0,:], sol1.y[1,:]-sol2.y[1,:], sol1.y[2,:]-sol2.y[2,:],'g')
ax2.set_xlabel('$x$'); ax2.set_ylabel('$y$'); ax2.set_zlabel('$z$')
plt.show()
#print("sol1=",sol1.y, '\n\n', "sol1-sol2=", sol1.y-sol2.y)