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

                     

贡献者: shizy0808

1. 微分方程组

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

\begin{equation} \begin{cases} x' = v_x\\ y' = v_y\\ v'_x = -GMx/(x^2 + y^2)^{3/2}\\ v'_y = -GMy/(x^2 + y^2)^{3/2} \end{cases}~ \end{equation}
的数值解,其中参数 $GM=1$。

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)
其中,输入参数分别为

   函数返回值分别为

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()

   对应的输出如图所示。

图
图 1:位置与速度关于时间的图形。
图
图 2:$x$ 与 $y$ 的相图。

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

                     

© 小时科技 保留一切权利