Simulating a Simple Pendulum - Increasing amplitude on each swing?

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
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)

• I edited the question to use LaTeX for equations instead of images. Feb 26 '15 at 13:21
• @Fidem, I had mentioned in a comment that there are cleaner ways to write RK4 code. Take a look at this answer for the way I usually do things. Ideally you would not have separate k1 and l1 variables. Mar 1 '15 at 17:25
• @DougLipinski Thank you very much. I initially looked to use vectors, but I couldn't find a data type/structure that could represent one... I guess I could use an array the size of the number of equations I am solving, and loop the RK4 method through that array Mar 4 '15 at 13:25
• I'm not familiar with Visual Basic, but I'd use an array type that supports scalar multiplication and array addition (like std::valarray in C++). Then you don't need to loop over the array. Your current method will work and it's clear enough, but doesn't generalize well to more equations. Mar 4 '15 at 16:08

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.

• I believe the problem lies here too. A is the "Driving Amplitude" and Omega is the "Driving Frequency". I am not 100% clear on what those terms are referring to so it is very possible that I have the wrong idea. Feb 27 '15 at 13:06
• After testing, I can confirm that it was that part of the function that was incorrect - So this is where my problem lies. Maybe Acos Ω is a "driving force" and is not needed when I am only using a damping effect? Feb 27 '15 at 13:14
• That's right. In general, the forcing term can actually be any function you want. A common choice would be simple periodic forcing which would take the form A*cos(Omega*t) (one would generally want time dependence in there). In this case A is the forcing amplitude and Omega is the angular frequency of the forcing. Note that resonance occurs when Omega is equal to the natural frequency of the system. Feb 27 '15 at 14:41

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.