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 · Answer 1 · 2020-04-20T05:35:52.103

0

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

edited Apr 20 '20 at 05:35

answered Apr 19 '20 at 08:16

Lutz Lehmann

25,219
2
22
51

so if i change xdot into an array using np.array, how would i be able to call each term of the array into the rk4. ive tried calling k1[0] after changing it into an array, but it says 'list' object is not callable – Ryan Apr 19 '20 at 08:38
You don't. Also, check again what argument is the scalar time and what is the state vector, `EOM` only has the state vector as argument. So you would call `k2 = EOM(xold+0.5*dt*k1)`, provided all tuples are converted to `np.array`. – Lutz Lehmann Apr 19 '20 at 08:43

proper way to code runge kutta 4th order steps

1 Answers1