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 行均为作图部分。


致读者: 小时百科一直以来坚持所有内容免费无广告,这导致我们处于严重的亏损状态。 长此以往很可能会最终导致我们不得不选择大量广告以及内容付费等。 因此,我们请求广大读者热心打赏 ,使网站得以健康发展。 如果看到这条信息的每位读者能慷慨打赏 20 元,我们一周就能脱离亏损, 并在接下来的一年里向所有读者继续免费提供优质内容。 但遗憾的是只有不到 1% 的读者愿意捐款, 他们的付出帮助了 99% 的读者免费获取知识, 我们在此表示感谢。

                     

友情链接: 超理论坛 | ©小时科技 保留一切权利