在rk45的true_divide中遇到无效值

Question

我正在尝试针对RK45来解决地球与太阳的两体问题，但不断除以我不理解的零。加速功能似乎在规范中发生了划分，但我不知道这种情况如何发生或如何解决。这是代码：

from scipy import optimize
from numpy import linalg as LA
import matplotlib.pyplot as plt
from scipy.optimize import fsolve
import numpy as np

AU=1.5e11
a=AU
e=0.5
mss=2E30
ms = 2E30
me = 5.98E24
mv=4.867E24
yr=3.15e7
h=100
mu1=ms*me/(ms+me)
mu2=ms*me/(ms+me)
G=6.67E11
step=24

vi=np.sqrt(G*ms*(2/(a*(1-e))-1/a))
#sun=sphere(pos=vec(0,0,0),radius=0.1*AU,color=color.yellow)
#earth=sphere(pos=vec(1*AU,0,0),radius=0.1*AU)

sunpos=np.array([-903482.12391302, -6896293.6960525, 0.  ])
earthpos=np.array([a*(1-e),0,0])

earthv=np.array([0,vi,0])
sunv=np.array([0,0,0])





def accelerations(t,earthposs, sunposs):
    norme=sum( (earthposs-sunposs)**2 )**0.5
    gravit = G*(earthposs-sunposs)/norme**3
    sunaa = me*gravit
    earthaa = -ms*gravit
    return earthaa, sunaa

def ode45(f,t,y,h):
        """Calculate next step of an initial value problem (IVP) of an ODE with a RHS described
        by the RHS function with an order 4 approx. and an order 5 approx.
        Parameters:
        t: float. Current time.
        y: float. Current step (position).
        h: float. Step-length.
        Returns:
        q: float. Order 2 approx.
        w: float. Order 3 approx.
        """

        s1 = f(t, y[0],y[1])
        s2 = f(t + h/4.0, y[0] + h*s1[0]/4.0,y[1] + h*s1[1]/4.0)
        s3 = f(t + 3.0*h/8.0, y[0] + 3.0*h*s1[0]/32.0 + 9.0*h*s2[0]/32.0,y[1] + 3.0*h*s1[1]/32.0 + 9.0*h*s2[1]/32.0)
        s4 = f(t + 12.0*h/13.0, y[0] + 1932.0*h*s1[0]/2197.0 - 7200.0*h*s2[0]/2197.0 + 7296.0*h*s3[0]/2197.0,y[1] + 1932.0*h*s1[1]/2197.0 - 7200.0*h*s2[1]/2197.0 + 7296.0*h*s3[1]/2197.0)
        s5 = f(t + h, y[0] + 439.0*h*s1[0]/216.0 - 8.0*h*s2[0] + 3680.0*h*s3[0]/513.0 - 845.0*h*s4[0]/4104.0,y[1] + 439.0*h*s1[1]/216.0 - 8.0*h*s2[1] + 3680.0*h*s3[1]/513.0 - 845.0*h*s4[1]/4104.0)
        s6 = f(t + h/2.0, y[0] - 8.0*h*s1[0]/27.0 + 2*h*s2[0] - 3544.0*h*s3[0]/2565 + 1859.0*h*s4[0]/4104.0 - 11.0*h*s5[0]/40.0,y[1] - 8.0*h*s1[1]/27.0 + 2*h*s2[1] - 3544.0*h*s3[1]/2565 + 1859.0*h*s4[1]/4104.0 - 11.0*h*s5[1]/40.0)
        w1 = y[0] + h*(25.0*s1[0]/216.0 + 1408.0*s3[0]/2565.0 + 2197.0*s4[0]/4104.0 - s5[0]/5.0)
        w2 = y[1] + h*(25.0*s1[1]/216.0 + 1408.0*s3[1]/2565.0 + 2197.0*s4[1]/4104.0 - s5[1]/5.0)
        q1 = y[0] + h*(16.0*s1[0]/135.0 + 6656.0*s3[0]/12825.0 + 28561.0*s4[0]/56430.0 - 9.0*s5[0]/50.0 + 2.0*s6[0]/55.0)
        q2 = y[1] + h*(16.0*s1[1]/135.0 + 6656.0*s3[1]/12825.0 + 28561.0*s4[1]/56430.0 - 9.0*s5[1]/50.0 + 2.0*s6[1]/55.0)

        return w1,w2, q1,q2
t=0
T=10**5
xarray=[]
yarray=[]
while t<T:
    ode45(accelerations,t,[earthpos,sunpos],h)
    earthpos=ode45(accelerations,t,[earthpos,sunpos],h)[1]
    sunpos=ode45(accelerations,t,[earthpos,sunpos],h)[3]
    xarray.append(ode45(accelerations,t,[earthpos,sunpos],h)[0][0])
    yarray.append(ode45(accelerations,t,[earthpos,sunpos],h)[0][1])
    print(ode45(accelerations,t,[earthpos,sunpos],h)[0][0],ode45(accelerations,t,[earthpos,sunpos],h)[0][1])
    t=t+h

plt.plot(xarray,yarray)
plt.savefig('orbit.png')
plt.show()

第二次迭代后，代码只返回了Earthpos的nan值。

Answer 1

数值积分方法通常积分一阶系统，y'=f(t,y)。您想要集成一个二阶ODE系统y''=f(t,y)，您首先需要将其转换为一阶系统。

为什么不使用numpy的向量类？

为什么您多次使用相同的参数执行相同的计算，而不是一次捕获所有返回值，然后将它们分配到列表中？

您也可以将scipy.integrate.solve_ivp与"RK45"方法一起使用，而不用自己编程。