title: Problems implementing Runge Kutta to solve a Damped Pendulum
body: <p>I am a high school student working on a "home project" to animate a damped pendulum by solving differential equations using the Runge Kutta method. <br><br>
   (The equations can be seen here: <a href="http://www.maths.tcd.ie/~smurray/Pendulumwriteup.pdf" rel="nofollow">http://www.maths.tcd.ie/~smurray/Pendulumwriteup.pdf</a>) <bR> <br>
  I have been informed that in my code, my implementation of RK4 is not correct, and to be honest I have been struggling to understand it.
  The program is written in VB.net 2010.
  My code for solving the equations are as follows:</p>
  
  <pre><code>Public Sub RK4Solve(ByRef The As Decimal, ByRef Ome As Decimal, ByRef h As Decimal)
  l1 = h * Ome
  k1 = h * f(The, Ome, h)
  l2 = h * (0.5 * l1) + Ome
  k2 = f((The + (0.5 * h * k1)), (Ome + (0.5 * h * l1)), (t + (0.5 * h)))
  l3 = h * (0.5 * l2) + Ome
  k3 = f((The + (0.5 * h * k2)), (Ome + (0.5 * h * l2)), (t + (0.5 * h)))
  l4 = h * l3 + Ome
  k4 = f((The + (h * k3)), (Ome + (h * l3)), (t + h))
  'Setting next step of variables
  The = The + (h / 6 * (l1 + (2 * l2) + (2 * l3) + l4))
  Ome = Ome + (h / 6 * (k1 + (2 * k2) + (2 * k3) + k4))
  t += h
  End Sub
  </code></pre>
  
  <p>I am aware that I am multiplying each step by too many h's - I am just lost on what is happening.</p>
  
  <p>My full code:</p>
  
  <pre><code>Public Class Form1
  
  Dim l As Decimal = 1 'Length of rod (1m)
  Dim g As Decimal = 9.81 'Gravity
  Dim w As Decimal = 0 ' Angular Velocity
  Dim initheta As Decimal = -Math.PI / 2 'Initial Theta
  Dim theta As Decimal = -Math.PI / 2 'Theta (This one changes for the simulation)
  Dim t As Decimal = 0 'Current time of the simulation
  Dim h As Decimal = 0.01 'Time step
  Dim b As Decimal = Math.Sqrt(g / l) 'Constant used in the function for dw/dt
  Dim k As Decimal = 0 'Coefficient of Damping
  Dim initialx = l * Math.Sin(initheta) 'Initial Amplitude of the pendulum
  
  Private Sub Form1_Load(ByVal sender As System.Object, ByVal e As System.EventArgs) Handles MyBase.Load
  End Sub
  
  'Function for dw/dt
  Public Function f(ByRef the As Decimal, ByRef omega As Decimal, ByRef time As Decimal)
  
      Return ((-b ^ 2) * Math.Sin(the)) - (k * omega) + (initheta * Math.Cos(omega * time))
  End Function
  
  Dim k1, k2, k3, k4, l1, l2, l3, l4 As Decimal 'Initialising RK4 variables
  
  
  Public Sub RK4Solve(ByRef The As Decimal, ByRef Ome As Decimal, ByRef h As Decimal)
  
      l1 = h * Ome
      k1 = h * f(The, Ome, h)
      l2 = h * (0.5 * l1) + Ome
      k2 = f((The + (0.5 * h * k1)), (Ome + (0.5 * h * l1)), (t + (0.5 * h)))
      l3 = h * (0.5 * l2) + Ome
      k3 = f((The + (0.5 * h * k2)), (Ome + (0.5 * h * l2)), (t + (0.5 * h)))
      l4 = h * l3 + Ome
      k4 = f((The + (h * k3)), (Ome + (h * l3)), (t + h))
  
      'Setting next step of variables
      The = The + (h / 6 * (l1 + (2 * l2) + (2 * l3) + l4))
      Ome = Ome + (h / 6 * (k1 + (2 * k2) + (2 * k3) + k4))
      t += h
  End Sub
  
  
  'Timer ticking every 0.1s
  'Time step is 0.01s to increase accuracy of results for testing
  Private Sub Timer1_Tick(ByVal sender As System.Object, ByVal e As System.EventArgs) Handles Timer1.Tick
      ComboBox1.Items.Add(theta) 'Adding theta to a drop down box to test data
      RK4Solve(theta, w, h)
  End Sub
  
  Private Sub Button1_Click(ByVal sender As System.Object, ByVal e As System.EventArgs) Handles Button1.Click
      Timer1.Enabled = False
  End Sub
  End Class
  </code></pre>
  
  <p>I have been trying to solve this for a while now, and I'm on my last legs so I am resorting to asking for help.  Thanks to anyone that can!</p>