proper way to code runge kutta 4th order steps

Question

xdot returns [dp_dt, dphi_dt] but the second step of RK4 method doesn't allow k1 to be called because k1 is a list...

def EOM(s)

    "Equations of Motion"
    p = s[0]
    phi = s[1]
    p_hat = p*b/(2*V)
    C = -0.06*phi+0.033*p_hat+0.073*p_hat**3-0.36*p_hat*phi**2+1.47*p_hat**2*phi
    dp_dt = 1/(2*Ixx)*rho*V**2*S*b*C
    dphi_dt = p 
    sdot = []
    sdot.append(dp_dt)
    sdot.append(dphi_dt)
    return sdot



def rk4(xold):


    xnew = []
    xdot = EOM(xold)
    for i, state in enumerate(xold):
        k1=xdot[i]
        #this is a list of [p,phi] ... how 
        k2=EOM(xold(1)+dt/2,xold(2)+dt/2*k1)
        k3=EOM(xold(1)+dt/2,xold(2)+dt/2*k2)
        k4=EOM(xold(1)+dt,xold(2)+dt*k3)
        xnew.append(state+dt/6*(k1+2*k2+2*k3+k4))
    return xnew

Lutz Lehmann Lutz Lehmann · Accepted Answer · 2020-04-19T08:16:14

Because it is a list. python does not do arithmetic vector operations on lists, for that use numpy to get arrays that act as vectors.

For an implementation of RK4 without numpy see Solving the Lorentz model using Runge Kutta 4th Order in Python without a package. It is not exactly the same as it uses a list of scalar ODE functions, but the principles should be visible.

Your code can be corrected to work without numpy, but that is only a solution for 2 or 3 variables, with more it becomes error-prone

def rk4(xold):
    k1 = EOM(xold)
    #this is a list of [p_dot,phi_dot]
    k2=EOM(xold[0]+dt/2*k1[0],xold[1]+dt/2*k1[1],) # trailing comma generates list
    k3=EOM([xold[0]+dt/2*k2[0],xold[1]+dt/2*k2[1]]) # or do it directly as list
    k4=EOM([xold[i]+dt*k3[i] for i in [0,1]]) # or use list generation
    return [ xold[i]+dt/6*(k1[i]+2*k2[i]+2*k3[i]+k4[i] for i in [0,1] ]))

proper way to code runge kutta 4th order steps

1 Answers