贡献者: shizy0808
这里我们再次以 “天体运动的简单数值计算” 中的问题为例,利用 Python 语言实现微分方程组
在对该问题求解之前,我们先简单介绍 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
的微分方程解的数据。
基于对 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
行均为作图部分。