我的RK4轨道繁殖器怎么了？

Question

我正在RK4轨道繁殖器中工作，并且遵循了：

k1 = dt×f（tn，yn）

k2 = dt×f（tn +dt/2，yn +k1/2）

k3 = dt×f（tn +dt/2，yn +k2/2）

k4 = dt×f（tn +dt，yn +k3）

yn+1 = yn+1/6（k1+2k2+2k3+k4）

tn+1 = tn+dt

我已经以更简单的速度来提出此问题，这是我为RK1实现而创建的：

#NOTE: sat is an array with the structure [x,y,z,vx,vy,vz,ax,ay,az]
def pasosimple(sat,dt):
   
    #physics constants  
    G = 6.6742e-20
    m_s = 1.989e30
    mu=G*m_s

    #Position Change from previous values of velocity and acceleration
    dx=sat[3]*dt+sat[6]*(dt**2)/2
    dy=sat[4]*dt+sat[7]*(dt**2)/2
    dz=sat[5]*dt+sat[8]*(dt**2)/2

   #Velocity change due to previous acceleration
    dvx=sat[6]*dt
    dvy=sat[7]*dt
    dvz=sat[8]*dt

   #xyz update
    x=dx+sat[0]
    y=dy+sat[1]
    z=dz+sat[2]
   
    #With the xyz update, we calculate new accelerations
    ax=(-mu*(x)/(np.sqrt((x)**2+(y)**2+(z)**2)**3))
    ay=(-mu*(y)/(np.sqrt((x)**2+(y)**2+(z)**2)**3))
    az=(-mu*(z)/(np.sqrt((x)**2+(y)**2+(z)**2)**3))
    
    #Substraction to obtain the difference acceleration
    dax=ax-sat[6]
    day=ay-sat[7]
    daz=az-sat[8]
   

    dsat=np.array([dx,dy,dz,dvx,dvy,dvz,dax,day,daz])           
    sat=np.array([x,y,z,dvx+sat[3],dvy+sat[4],dvz+sat[5],ax,ay,az])
    return dsat,sat

据我所知，此代码有效，并且我已经对其进行了测试。 现在，为了掩盖RK4，我正在这样做：

def rk4(sat,dt):
    d1,s1=pasosimple(sat,dt)            
    d2,s2=pasosimple(sat+d1/2,dt+0.5*dt)    
    d3,s3=pasosimple(sat+d2/2,dt+0.5*dt)
    d4,s4=pasosimple(sat+d3,dt+dt)     
    sat=sat+(1/6)*(d1+d2*2+d3*2+d4) 

    return sat

它不起作用。 我希望任何人都可以对我做错了什么有所了解。 谢谢大家。

编辑：

我已经尝试了此替代配置的pasosimple，以下（或尝试遵循）以下注释。结果也无效。但是，仅使用新的pasosimple代码，而不是将其嵌套到rk4中效果很好，这使我认为问题现在在rk4中。

def pasosimple(sat,dt):
    G = 6.6742e-20
    m_t=5.972e24
    m_s = 1.989e30
    m_m=6.39e23
    mu=G*m_s
    mu_t=G*m_t
    mu_m=G*m_m
    
    ax=(-mu*(sat[0])/(np.sqrt((sat[0])**2+(sat[1])**2+(sat[2])**2)**3))
    ay=(-mu*(sat[1])/(np.sqrt((sat[0])**2+(sat[1])**2+(sat[2])**2)**3)) 
    az=(-mu*(sat[2])/(np.sqrt((sat[0])**2+(sat[1])**2+(sat[2])**2)**3))
                        
    x_punto=np.array([sat[3],sat[4],sat[5],ax,ay,az]) 
    x=np.array([sat[0]+x_punto[0]*dt,sat[1]+x_punto[1]*dt,sat[2]+x_punto[2]*dt,sat[3]+x_punto[3]*dt,sat[4]+x_punto[4]*dt,sat[5]+x_punto[5]*dt])
    
    
    dsat=np.array([x_punto[0]*dt,x_punto[1]*dt,x_punto[2]*dt,x_punto[3]*dt,x_punto[4]*dt,x_punto[5]*dt,ax-sat[6],ay-sat[7],ay-sat[8]])
    sat=np.array([x[0],x[1],x[2],x[3],x[4],x[5],ax,ay,az])
    return dsat,sat