SciPy 求解常微分方程组的初值问题

1. 微分方程组

　　 这里我们再次以 “天体运动的简单数值计算” 中的问题为例，利用 Python 语言实现微分方程组

2. solve_ivp 求解器介绍

　　 在对该问题求解之前，我们先简单介绍 Python 中关于微分方程数值求解库极相关函数的用法。一般，最常用的数值计算库为 scipy,而微分方程求解对应的模块为 scipy.integrate.solve_ivp.

　　 solve_ivp 的一般格式为

solve_ivp(fun, t_span, y0, method='RK45', t_eval=None,max_step=np.inf,
        dense_output=False, events=None, 
        vectorized=False, args=None, **options)

• fun 是微分方程（组）的右端；
• t_span=(t0,tf) 代表积分区间为 t0 到 tf;
• y0 为初始条件；
• method 为数值积分方法，这里可以使用的积分方法有 RK45、RK23、DOP853(8 阶显式龙格库塔法)、Radau(5 阶 Radau IIA 族的隐式 Runge-Kutta 方法)、BDF(隐式多步变阶)、LSODA(具有自动刚度检测和切换的 Adams/BDF 方法)等。需要主要的是显式 Runge-Kutta 方法（RK23、RK45、DOP853）应用于非刚性问题，隐式方法（Radau、BDF）应用于刚性问题。
• t_eval 是可选参数，需要数组类型数据。如果给定就在这些时间点进行求解，否则求解器自己选择时间点进行求解。
• max_step 允许的最大步长。默认为 np。inf，即步长不受约束，仅由解算器确定。如果求解不能满足需求，可以手动改变最大步长已达到精度要求。相应的还有最小步长参数 min_step.
• dense_output 为布尔类型，默认为假，是否计算连续解。
• args 为元组（tuple）类型的参数，用于向微分方程传递必要的参数
• 后面还有许多可选参数，很少会使用。

　　 函数返回值分别为

• t 返回计算的时间点数据，是一个 ndarray 数据类型，长度为 (n_points,).
• y 大小为 (n, n_points) 的 ndarray 的微分方程解的数据。
• 还有许多返回参数，感兴趣读者可以去官网查看。

3. 代码实现

　　 基于对 solve_ivp 的使用说明，我们接下来对微分方程组\ref{eq_PyIVP_10} 进行数值计算。具体代码如下：

# 导入必要的库
from scipy.integrate import solve_ivp
import  numpy as np
import matplotlib.pyplot as plt

# 定义一些参数
GM = 1 # 万有引力常数乘以中心天体质量
x0,y0 = 1,0 # 初始位置
vx0,vy0 = 0,0.7 # 初始速度
tspan = (0, 4) # 总时间和步数
init0=(x0,y0,vx0,vy0)

# 定义微分方程
def odefun(t,z,GM):
    x,y,vx,vy  = z
    temp = -GM /(x**2 + y**2)**(3/2)
    return [vx,vy, temp * x, temp * y]
# 调用求解器求解微分方程
sol = solve_ivp(odefun,tspan,init0,dense_output=True,
          max_step=0.001,args=(GM,))
t = np.linspace(0, 4, 3000)
z = sol.sol(t)


# 微分方程解的可视化

# xy的时序图
plt.plot(t, z.T)
plt.xlabel('t')
plt.legend(['x', 'y','vx','vy'], shadow=True)
plt.show()

# xy的相图
plt.plot(z[0],z[1])
plt.xlabel('t')
plt.xlabel('x')
plt.ylabel('y')
plt.title('Lotka-Volterra System')
plt.show()

　　 对应的输出如图所示。

　　 我们现在来看代码的第 2-4 行是一些基本必要库的导入。 7-9 行为方程的一些参数与初始条件。 10 行为积分区间。注意它是一个元包数组，也可以写成 [t0,tf] 形式。 14-17 行为微分方程的定义，这里需要传递一个参数 GM. 19 行为调用求解器求解微分方程组，我们使用了 dense_output 参数，这样我们就可以在连续的对时间取值计算对应的 x,y,vx,vy. 28-39 行均为作图部分。