0

I am working in a RK4 orbital propagator, and I have followed the basic structure of:

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

I have bassed this procces in a simpler one, created by me for RK1 implementation:

#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

This code works, as far as I know, and I have already tested it. Now, for imlementing the RK4, I'm doing this:

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

Which does not work. I was hoping that anyone could give some insight about what i am doing wrong. Thank you all.

EDIT:

I have tried this alternate configuration for pasosimple, following (or trying to follow) the comments below. The result did not work either. However, using only the new pasosimple code and not nesting it into rk4works nicely, which makes me think the problem is now in 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
dfernr
  • 1
  • 1
  • Nothing is wrong with the RK4 steps. But the `pasosimple` is too much of a mixture of velocity Verlet and derivatives function. The derivatives are evaluated at the wrong point for RK4 so that its order is reduced, possibly to error order 1. It would be better to separate the acceleration computation out and have separate procedures for the Verlet loop and the RK4 loop. – Lutz Lehmann Apr 29 '22 at 07:51
  • Thank you for your comments, I'm sure there will be very useful. Unfortunately I was not able to fully understand it porperly. Could yo elaborate your answer a bit more? Thank you. – dfernr Apr 30 '22 at 08:41
  • For the RK4 computation you need the derivatives at the current state. However, in your implementation you also do parts of a Verlet step, invalidating the internals that give the order of the RK4 step. Also, in the RK4 step there is no 3rd derivative, the provisions you made for the velocity Verlet step where you indeed need the old acceleration are not needed in RK4. – Lutz Lehmann Apr 30 '22 at 11:36

0 Answers0