Problems implementing Runge Kutta to solve a Damped Pendulum

Question

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.

(The equations can be seen here: http://www.maths.tcd.ie/~smurray/Pendulumwriteup.pdf)

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:

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

I am aware that I am multiplying each step by too many h's - I am just lost on what is happening.

My full 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

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!

Lutz Lehmann Lutz Lehmann · Accepted Answer · 2015-02-26T21:51:57

It helps if you separate the differential equations from the RK4 implementation. Then you can implement RK4 as in the documentation

k1 = f1( y1, y2, x)
l1 = f2( y1, y2, x)
k2 = f1( y1 + 0.5*h*k1, y2 + 0.5*h*l1, x + 0.5*h)
l2 = f2( y1 + 0.5*h*k1, y2 + 0.5*h*l1, x + 0.5*h)
k3 = f1( y1 + 0.5*h*k2, y2 + 0.5*h*l2, x + 0.5*h)
l3 = f2( y1 + 0.5*h*k2, y2 + 0.5*h*l2, x + 0.5*h)
k4 = f1( y1 + h*k3, y2 + h*l3, x + h)
l4 = f2( y1 + h*k3, y2 + h*l3, x + h)

y1 = y1 + h/6*(k1+2*(k2+k3)+k4)
y2 = y2 + h/6*(l1+2*(l2+l3)+l4)
x = x + h

Not only does it help to avoid the duplication of the multiplication by h, but also the duplication of the addition of the base point values Ome and The.

Problems implementing Runge Kutta to solve a Damped Pendulum

1 Answers