2

The implementation of scipy.integrate.RK45 does not specify where to mention the times at which the integration should be performed.

The input option "t_bound" seems to be the final epoch.

To integrate, one must use the "step" option, which has no input, so one can't specify the time stamp/s there.

Compared to the other integrator: scipy.integrate.solve_ivp, which has the option "t_eval" where one can mention all the epochs (as an array) where the solution is needed.

The following code works:

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

with "fun_integrator" as my function with input (t,y) t_span = [0 max(t)] and

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

^= epochs at which the answer is needed

However, Matlab ode113 gives better results than scipy solve_ivp, so I wanted to try scipy.integrate.RK45 to check if I get better results with it.

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

BUT how to input step size in the above code?

For more info: Function and 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])
Vishu Singh
  • 21
  • 1
  • 2
  • How do you know that the matlab results (implicit multistep with adaptive order 1-13) are better? Why do you think that explicitly replicating what `solve_ivp` does internally with the `RK45` class will give substantially different results? Did you try to use other solver classes like `LSODA` or `RADAU` in `solve_ivp`? – Lutz Lehmann Jun 20 '19 at 12:54
  • Yes, I did try ```LSODA``` and the other options in ```solve_ivp```, they all have more or less similar results. I know the MATLAB results (both with ode45 and ode113, with a stepsize of 1800) are better because I have values for comparison for different epochs. I'm assuming that one of the integrators (```RK45```, ```RK23```, etc.) could possibly give results as good as MATLAB, or better. Since MATLAB uses on 16 digits for calculation, and my results need to be highly accurate, I'm hoping to implement everything in python to get better accuracy at the end. – Vishu Singh Jun 20 '19 at 22:17
  • You might want to change `t = np.linspace(0,0+0.3*(pts-1)*86400,20)` to `t = np.linspace(0,0+0.3*(pts-1)*86400,pts)` for consistency when changing `pts` from the value `20`. – Lutz Lehmann Jun 20 '19 at 22:38
  • While `solve_ivp` and the classes behind it are pure python, if I remember correctly what I once read, I'm not so certain that it will work flawlessly with multi-precision data types. You might have to use your own implementation for that, if I understand your goal of higher accuracy correctly. – Lutz Lehmann Jun 20 '19 at 22:42

1 Answers1

2

To answer the actual question and not the underlying question:

You can use the dense_output method to calculate the value at a particular desired time point. Otherwise you will lose a lot of accuracy in your solution if you, for example, try to use your own interpolation method.

But, I would note that solve_ivp (using RK45 method) is essentially an easy to use wrapper for this code, so I wouldn't expect anything significantly different this way.

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

Matthew Flamm
  • 448
  • 3
  • 13