title: Using Runge-Kutta method to solved differential equation of motion
body: <p>In our physics class we have to model a damping torsional pendulum.</p>
  <p>Ilustration of torsional pendulum:
  <img src="https://i.stack.imgur.com/QiKkK.png" alt="Torsional Pendulum" /></p>
  <p>We came up with this equation of motion:
  <img src="https://i.stack.imgur.com/Mpfgr.png" alt="Equation of motion" /></p>
  <p>Where <em><strong>θ</strong></em> is the angle, <strong>A</strong> is <em>torsion parameter</em>, <strong>B</strong> is <em>Newton's parameter</em>, <strong>C</strong> is <em>Stokes' parameter</em>, and <strong>D</strong> is <em>friction parameter</em>. We also use the sign function <em>sgn</em> that determines the direction of the acting force upon the pendulum, depending on the current angle from the reference point.</p>
  <p><strong>The problem is</strong>, that I'm unable to solve it using <em><strong>Runge-Kutta</strong></em> method in <em>Python</em>.</p>
  <p>I got a working solution in <em>MATLAB</em> by using Euler's method, which has some flaws, but it is something.</p>
  <p><em><strong>MATLAB Code:</strong></em></p>
  <pre><code>function [theta, dtheta, epsilon] = drt(t, theta0, dtheta0, A, B, C, D)
      epsilon = zeros(1, length(t));
      theta = epsilon;
      dtheta = epsilon;
  
      theta(1) = theta0;
      dtheta(1) = dtheta0;
  
  
      epsilon(1) = A * alpha0 - B * dtheta^2 - C * dtheta - D;
  
      dt = t(2) - t(1);
  
      for i = 1 : (length(t) - 1)
         epsilon(i + 1)= - A * theta(i) - B * dtheta(i)^2 * sign(dtheta(i)) - C * dtheta(i) - D * sign(dtheta(i));
         dtheta(i + 1)= dtheta(i) + dt * epsilon(i); 
         theta(i + 1) = theta(i) + dt * dtheta(i);
      
      end
  end
  </code></pre>
  <p>We call this <em>MATLAB</em> function like this for example:</p>
  <pre><code>t = linspace(0, 10, 100);
  theta0 = 90;
  dtheta0 = 0;
  A = 1;
  B = 0.1;
  C = 0.1;
  D = 0.1;
  [theta, dtheta, epsilon] = drt(t, theta0, dtheta0, A, B, C, D);
  </code></pre>
  <p>We can then plot the <em>theta</em> and other values in a graph, which shows us, how the torsional pendulum is being damped by the external forces acting on it.</p>
  <p><em><strong>Python Code:</strong></em></p>
  <pre><code>import numpy as np
  import matplotlib.pyplot as plt
  
  
  # Damping torsional pendulum
  def drp(drp_alpha, drp_d_alpha, drp_params):
      a = drp_params[&quot;tors&quot;]
      b = drp_params[&quot;newt&quot;]
      c = drp_params[&quot;stok&quot;]
      d = drp_params[&quot;fric&quot;]
  
      result = a * drp_alpha - b * np.power(drp_d_alpha, 2) * np.sign(drp_d_alpha) - c * drp_d_alpha - d * np.sign(drp_d_alpha)
  
      return result
  
  
  # Runge-Kutta 4th Order
  # f - function DamRotPen
  # x0 - initial condition
  # t0 - initial time
  # tmax - maximum time
  # dt - sample time
  def RG4(rg4_f, rg4_x0, rg4_t0, rg4_tmax, rg4_dt):
      # Time vector
      rg4_t = np.arange(rg4_t0, rg4_tmax, rg4_dt)
      # Time vector size
      rg4_t_sz = rg4_t.size
  
      # Initialize the array
      rg4_alpha = np.zeros(rg4_t_sz)
      # Initial value of the system
      rg4_alpha[0] = rg4_x0
  
      for k in range(rg4_t_sz - 1):
          k1 = dt * f(rg4_t[k], rg4_alpha[k])
          k2 = dt * rg4_f(rg4_t[k] + dt / 2, rg4_alpha[k] + k1 / 2)
          k3 = dt * rg4_f(rg4_t[k] + dt / 2, rg4_alpha[k] + k2 / 2)
          k4 = dt * rg4_f(rg4_t[k] + dt, rg4_alpha[k] + k3)
  
          rg4_d_alpha = (k1 + 2 * k2 + 2 * k3 + k4) / 6
  
          rg4_alpha[k + 1] = rg4_alpha[k] + rg4_d_alpha
  
      return rg4_alpha, rg4_t
  
  
  # Parameters of the forces acting on the system
  # tors - torsion parameter
  # newt - Newton's parameter
  # stok - Stokes' parameter
  # fric - friction parameter
  params = {&quot;tors&quot;: 1, &quot;newt&quot;: 0.1, &quot;stok&quot;: 0.1, &quot;fric&quot;: 0.1}
  
  # Start parameters
  alpha = 90
  d_alpha = 0
  
  # Initial time
  t0 = 0
  # Maximum time
  tmax = 120
  # Sample time
  dt = 0.01
  
  # Define DamRotPen function as 'f' using lambda
  f = lambda t, alpha : drp(alpha, d_alpha, params)
  
  # Try to solve this shit
  alpha, t = RG4(f, alpha, t0, tmax, dt)
  
  # Plot this shit
  plt.plot(t, alpha, &quot;r&quot;, &quot;r&quot;, label=&quot;Position I guess&quot;)
  plt.xlabel(&quot;Time t / s&quot;)
  plt.grid()
  plt.show()
  </code></pre>
  <p>After plotting the values, we can see on the graph, that the <em><strong>θ</strong></em> sky rockets after certain amount of time. I don't know what I'm doing wrong, I tried practically everything, so that's why I'm asking you for help. (Though I think part of the problem might be my misunderstanding on how to implement the Runge-Kutta method, maybe I got the math wrong, etc...)</p>