Runge-kutta 4th Order To Solve 2nd Order Ode System Using Python

I'm trying to solve system of two odes numerically by runge-kutta 4th order method. initial system: system to solve: And I have very strange solution graph... I have: Correct

Solution 1:

Your state has 4 components, thus you need 4 slopes in each stage. It should be obvious that the slope/update for z can not come from k1_zdot, it has to be k1_z which is to be computed previously as

k1_z = zdot

and in the next stage

k2_z = zdot + dt/2*k1_zdot

etc.

---

But better is to use a vectorized interface

def derivs(t,u):
    z, theta, dz, dtheta = uddz= -omega_z * z - epsilon / 2 / m * thetaddtheta= -omega_theta * theta - epsilon / 2 / I * z 
    return np.array([dz, dtheta, ddz, ddtheta]);

and then use the standard formulas for RK4

i = 0whileTrue:
    # runge_kutta
    k1 = derivs(t[i], u[i])
    k2 = derivs(t[i] + dt/2, u[i] + dt/2 * k1)
    k3 = derivs(t[i] + dt/2, u[i] + dt/2 * k2)
    k4 = derivs(t[i] + dt, u[i] + dt * k3)

    u.append (u[i] + (k1 + 2*k2 + 2*k3 + k4) * dt / 6)
    i += 1

and later unpack as

z, theta, dz, dtheta = np.asarray(u).T