如何在 scipy.integrate.RK45 中输入时间步长答案

【问题标题】：How to input time steps in scipy.integrate.RK45如何在 scipy.integrate.RK45 中输入时间步长
【发布时间】：2019-11-03 03:40:15
【问题描述】：

scipy.integrate.RK45 的实现没有指定在哪里提到应该执行集成的时间。

输入选项“t_bound”似乎是最后一个纪元。

要集成，必须使用“step”选项，该选项没有输入，因此无法指定时间戳。

与其他积分器相比：scipy.integrate.solve_ivp，它具有选项“t_eval”，可以在其中提及需要解决方案的所有时期（作为数组）。

以下代码有效：

check=solve_ivp(fun_integrator, t_span, y0, \
             method='RK45', t_eval=t, max_step=1800, rtol=10**(-11), atol=10**(-12))

以“fun_integrator”作为我的函数，输入 (t,y) t_span = [0 max(t)] 和

t = np.linspace(0,0+0.3*(pts-1)*86400,20)

^= 需要答案的时期

但是，Matlab ode113 比 scipy solve_ivp 提供更好的结果，所以我想尝试 scipy.integrate.RK45 来检查我是否得到更好的结果。

check45=RK45(fun_integrator, 0, y0, max(t), max_step=1800, \
             rtol=1e-11, atol=1e-12)
for i in range(len(t)):
    dt=0+0.3*i*86400
    check45.step() #step has no input
    x_E45=check45.y
    print(x_E45[0:3])

但是如何在上面的代码中输入步长？

更多信息：函数和y0：

def fun_integrator(t,y):
    G=6.67428*10**(-11)/(10**3)**3 #Km**3/(Kg s**2) 
    m_E=5.97218639014246e+24 #Kg
    m_M=7.34581141668673e+22 #Kg
    m_S=1.98841586057223e+30 #Kg
#    c=2.997924580000000e+05 # [Km/s], speed of light, IERS TN36 p. 18
    x_E=y[0:3]
    x_M=y[6:9]
    x_S=y[12:15]
    ##
#    v_E=y[0,3:6]
#    v_M=y[0,9:12]
#    v_S=y[0,15:18]
    ##
    diff_EM=x_E-x_M
    diff_SM=x_S-x_M
    diff_ES=x_E-x_S
    d_EM=math.sqrt(diff_EM[0]**2 + diff_EM[1]**2 + diff_EM[2]**2)
    d_SM=math.sqrt(diff_SM[0]**2 + diff_SM[1]**2 + diff_SM[2]**2)
    d_ES=math.sqrt(diff_ES[0]**2 + diff_ES[1]**2 + diff_ES[2]**2)
    ##
    dydt= y[3:6] #E vel
    two=G*((m_M*(-diff_EM)/d_EM**3) + (m_S*(-diff_ES)/d_ES**3)) # E acc
    three=y[9:12] #M vel
    four=G*((m_S*(diff_SM)/d_SM**3) + (m_E*(diff_EM)/d_EM**3)) # M acc
    five=y[15:18] #S vel
    six=G*((m_E*(diff_ES)/d_ES**3) + (m_M*(-diff_SM)/d_SM**3)) # S acc
    dydt=np.append(dydt,two)
    dydt=np.append(dydt,three)
    dydt=np.append(dydt,four)
    dydt=np.append(dydt,five)
    dydt=np.append(dydt,six)
    return dydt

y0=np.array([0.18030617557041088626562258966E+08, #Epos
-0.13849983879561028969707483658E+09,
-0.60067586558743159549398380427E+08,
0.29095339700433727181771510027E+02, #Evel
0.30306447236947439878791802685E+01,
0.13145956648460951506322034682E+01,
0.17909715946785370737211415301E+08, #Mpos
-0.13879823119464902933064808964E+09,
-0.60230238740754862546525906740E+08,
0.30136092114725188798503502634E+02, #Mvel
0.27407201232607066962011074680E+01,
0.11664485102480685879418310910E+01,
0.67356572699027185845872823900E+06, #Spos
0.11475300015697844979868843500E+06,
0.39801697980815553718511148000E+05,
-0.60903913908114150920444000000E-03, #Svel
0.89648366457310610847232000000E-02,
0.38595942076153950302606500000E-02])

【问题讨论】：

你怎么知道matlab结果（隐式多步，自适应阶数1-13）更好？为什么你认为在 RK45 类内部显式复制 solve_ivp 所做的事情会产生截然不同的结果？您是否尝试在 solve_ivp 中使用其他求解器类，例如 LSODA 或 RADAU？
是的，我确实尝试了LSODA 和solve_ivp 中的其他选项，它们都有或多或少相似的结果。我知道 MATLAB 结果（使用 ode45 和 ode113，步长为 1800）更好，因为我有不同时期的比较值。我假设其中一个积分器（RK45、RK23 等）可能会提供与 MATLAB 一样好的结果，或者更好。由于 MATLAB 使用 16 位进行计算，并且我的结果需要非常准确，我希望在 python 中实现所有内容，以便最终获得更好的准确性。
在将 pts 从值 20 更改时，您可能希望将 t = np.linspace(0,0+0.3*(pts-1)*86400,20) 更改为 t = np.linspace(0,0+0.3*(pts-1)*86400,pts) 以保持一致性。
虽然solve_ivp 及其背后的类是纯python，但如果我没记错我曾经读过的内容，我不确定它是否能完美地处理多精度数据类型。如果我正确理解您的更高准确度目标，您可能必须使用自己的实现。

标签： python scipy ode orbital-mechanics integrator

【解决方案1】：

回答实际问题而不是根本问题：

您可以使用 dense_output 方法计算特定所需时间点的值。否则，如果您尝试使用自己的插值方法，您的解决方案将失去很多准确性。

但是，我要注意的是，solve_ivp（使用 RK45 方法）本质上是一个易于使用的代码包装器，所以我不希望这种方式有什么显着不同。

https://docs.scipy.org/doc/scipy/reference/generated/scipy.integrate.RK45.dense_output.html#scipy.integrate.RK45.dense_output

【讨论】：