title: Double Pendulum RK4
body: <p>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:</p>
  <p><a href="https://i.stack.imgur.com/qlIDX.png" rel="nofollow noreferrer"><img src="https://i.stack.imgur.com/qlIDX.png" alt="enter image description here" /></a></p>
  <p>Since this is a system of first order differential equations, we are ready to use RK4 to solve it.
  However, <strong>I am really lost in regards to writing the system above in python code.</strong></p>
  <p>There is a related question here: <a href="https://stackoverflow.com/questions/64962719/system-of-second-order-odes-runge-kutta-4th-order">System of second order ODEs Runge Kutta 4th order</a> 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 <a href="https://stackoverflow.com/questions/65224923/i-want-to-have-the-pendulum-blob-in-my-double-pendulum">I want to have the pendulum blob in my double pendulum</a> 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.</p>
  <pre><code>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] &lt; 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()
  </code></pre>
  <p>However, the terminal is giving me the following output:</p>
  <pre><code>[array([0.78539816, 0.78539816, 0.        , 0.        ])]
  </code></pre>
  <p>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.</p>
  <p>I am really really struggling to make this work. Can someone help me fix the code as to appropriately integrate the system of ODES?</p>
  <p>Thanks in advance, Lucas</p>