Double Pendulum RK4

Question

I am trying to implement the classical RK4 algorithm to solve the system of differential equations which govern the equations of motion. However, I am having several problems for a while. The equations of motion can be written as below:

Since this is a system of first order differential equations, we are ready to use RK4 to solve it. However, I am really lost in regards to writing the system above in python code.

There is a related question here: System of second order ODEs Runge Kutta 4th order which writes the system of ODES for a different physical system into Python code and also uses RK4 to solve it. However, I have not been able to use that to make my own code work. There is another question I want to have the pendulum blob in my double pendulum which is closer to what I want, but it does not take care of the process by using arrays, but rather by writing several equations instead.

import matplotlib.pyplot as plt
import numpy as np


def RK4(t, y, h, f):  
    
    #Runge Kutta standard calculations
    k1 = f(t, y)
    k2 = f(t + h/2, y + h/2 * k1)
    k3 = f(t + h/2, y + h/2 * k2)
    k4 = f(t + h, y + h * k3)    

    return 1/6*(k1 + 2 * k2 + 2 * k3 + k4) 

    
def RHS(t, y):
    
    ################################
    #Critical physical Parameters
    g = 9.8
    l1 = 1
    l2 = 1
    m1 = 1
    m2 = 1
    w1 = w2 = 0
    theta_1 = theta_2 = np.pi/4
    delta_theta = theta_1 - theta_2
    ################################
    
    #Writing the system of ODES.
    f0 = w1
    f1 = w2

    f2 = (m2*l1*w1**2*np.sin(2*delta_theta) + 2*m2*l2*w2**2*np.sin(delta_theta) + 2*g*m2*np.cos(theta_2)*np.sin(delta_theta)
    + 2*g*m1*np.sin(theta_1))/(-2*l1*(m1 + m2*np.sin(delta_theta)**2))
    
    f3 = (m2*l2*w2**2*np.sin(2*delta_theta) + 2*(m1 + m2)*l1*w1**2*np.sin(delta_theta) 
    + 2*g*(m1 + m2)*np.cos(theta_1)*np.sin(delta_theta))/(2*l2*(m1 + m2*np.sin(delta_theta)**2))
    
    return np.array([f0, f1, f2, f3])


def main():
    #Specify time interval and number of domain subintervals
    t0 = 0 
    T = 100
    n = 2048
    # initial conditions
    y0 = [np.pi/4, np.pi/4, 0, 0]                  
    
    #Domain discretization
    t_n = np.linspace(t0, T) 
    y_n = [np.array(y0)]  
    
    #Step size
    h  = (T - t0)/n        
    
    while t_n[-1] < T:
        #Keep going until T is reached.
        #Store numerical solutions into y_n
        y_n.append(y_n[-1] + h * RK4(t_n[-1], y_n[-1], h, RHS))
        t_n.append(t_n[-1] + h)
    
    
    print(y_n)
 

main()

However, the terminal is giving me the following output:

[array([0.78539816, 0.78539816, 0.        , 0.        ])]

which probably indicates the system of equations is not being solved. Why is that? I guess that I am not appropriately passing the initial conditions also, and I just dont know how to do it properly.

I am really really struggling to make this work. Can someone help me fix the code as to appropriately integrate the system of ODES?

Thanks in advance, Lucas

Answer 1

You are not using the passed state variable in the derivatives function, you just set the state to a constant. You should have a line like

theta1, theta2, w1, w2 = y

instead.