Simulating a Simple Pendulum - Increasing amplitude on each swing?

Question

I am trying to animate a dampened pendulum using RK4 for a highschool project.

The equations that describe the dampened system are as follows: (from http://www.maths.tcd.ie/~smurray/Pendulumwriteup.pdf) $$\begin{aligned} \frac{d\theta}{dt}&=\omega,\\ \frac{d\omega}{dt}&=-\beta^2\sin\theta-k\omega+A\cos\Omega. \end{aligned}$$

Where $B$ is a constant, $k$ is the coefficient of damping, $A$ is the driving amplitude and $\Omega$ is the driving frequency.

My problem is that I am getting an increasing amplitude on each swing, varying the coefficient of damping does not change this.

I believe my problem is either my programming is wrong, or I have the wrong idea of what the driving amplitude and the driving frequency are.
Ideally, I want the program to be able to showcase simple harmonic motion (when $k = 0$) and also allow the user to investigate dampening.

This is my first time attempting to use the Runge Kutta method and I have most likely incorrectly implemented it.

My code is written in Visual Basic.net 2010 and is as follows:

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

Public Function y(ByRef the As Decimal, ByRef omega As Decimal, ByRef time As Decimal)
    Return omega
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 = y(The, Ome, t)
    k1 = f(The, Ome, t)
    l2 = y(The + (0.5 * h * l1), Ome + (0.5 * h * k1), t + (0.5 * h))
    k2 = f(The + (0.5 * h * l1), Ome + (0.5 * h * k1), t + (0.5 * h))
    l3 = y(The + (0.5 * h * l2), Ome + (0.5 * h * k2), t + (0.5 * h))
    k3 = f(The + (0.5 * h * l2), Ome + (0.5 * h * k2), t + (0.5 * h))
    l4 = y(The + (h * l3), Ome + (h * k3), t + h)
    k4 = f(The + (h * l3), Ome + (h * k3), 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

Here is a picture of the Dis, Vel and Acc graphs (in that order)
As you can see, the acceleration falls apart during the simulation (forgive my lack of scientific terms) Why does this happen? (I will also update my code with my new RK4 implementation)

Answer 1

It looks like your revised code still has a problem. What is this term (initheta * Math.Cos(omega * time)) in your function f() supposed to be? It's like a forcing term that has amplitude equal to the initial angle and a time-dependent frequency equal to the instantaneous angular frequency. Shouldn't it be (A * Math.Cos(Omega * time)) and you have to choose A and Omega?

With the (incorrect) term you have now and your initial conditions this forcing will initially act in opposition to the motion, slowing the pendulum. Weird things might happen eventually since the forcing frequency is constantly changing.

Answer 2

You are multiplying $l$'s and $k$'s by too many $h$'s. In particular look out for expressions like $l_2 = h l_1/2+\Omega$, where $\Omega$ is $O(1)$ but $l_1$ is $O(h)$ and $l_2$ is $O(h^2)+O(1)$. This is not a correct implementation of RK4.

The idea of RK methods is that you have terms like $$ k_2 = f(t_0 + c_2h, y_0 + h a_{21}k_1), $$ but $k_2$ has the magnitude of the rhs of the equation, it is an approximation to the rhs at some point, not like in your code.

Also, consider asking on http://codereview.stackexchange.com.