Numerically solving the equation of a simple pendulum with Runge-Kutta.

Question

I am trying to solve the equation $\dfrac{\mathrm{d}^2\theta}{\mathrm{d}t^2} + \dfrac{g}{L} \sin{\theta} = 0$ using Runge-Kutta. I have alread split it into the following equations

$\dfrac{\mathrm{d}\omega}{\mathrm{d}t} = - \dfrac{g}{L} \sin{\theta}$

$\dfrac{\mathrm{d}\theta}{\mathrm{d}t} = \omega$

I have implemented it correctly for the most part however It doesn't appear to care what value of theta I put in initially, the period does work correctly (at least it appears this way) and is dependent on the omega I put in initially. The angle also increases first, ideally I want the first frame to be where the pendulum is at its maximum angle and for it to fall down from there. The maximum angle (no matter what I set it to initially) looks to be around 1 RAD.

I'm only a high school student so I'm a bit lost, I tried to solve the equations analytically at first however I haven't read far into elliptic integrals.

Just a quick note: a ShmpIntance object will be initialised in the main procedure and then a frame and canvas will be created which call the getGFX method repeatedly where x is incremented each by 1 each cycle. There some code in action behind the scenes making sure there is a maximum of 60 cycles per seconds (60FPS).

package sim;

import java.awt.Color;
import java.awt.Graphics;
import java.text.DecimalFormat;

public class ShmpInstance extends SimInstance {

private double omega, theta, length;
private double currentT, lastT;
private int xPosPix, yPosPix;
private double a, b, c, d, h;

private int lengthInPix, lengthFInPix;

private DecimalFormat df = new DecimalFormat("#.#");

public ShmpInstance(double length, double theta_0, int width, int height) {

    super(width, height);

    this.theta = theta_0;
    this.omega = Math.sqrt(g / length);
    this.length = length;

    if(8*(height/10) > 4*(width/10)){
    lengthInPix = 4*(width/10);
    } else {
    lengthInPix = 8*(height/10);
    }

    if(height > width){ 
        lengthFInPix = 2*(width/25);
    } else {
        lengthFInPix = 2*(height/25);
    } 

}

public Graphics getGFX (Graphics g, int x){

    lastT = currentT;
    currentT = (double)(x) / 60;
    h = (currentT-lastT);

    RungeKutta(theta, omega);

    xPosPix = (int)(width/2 + lengthInPix*Math.sin(theta));
    yPosPix = (int)(height/10 + lengthInPix*Math.cos(theta));

    g.setColor(Color.RED);
    g.drawLine(xPosPix, yPosPix, xPosPix - (int)(lengthFInPix*Math.sin(theta)*Math.cos(theta)), yPosPix + (int)(lengthFInPix*Math.sin(theta)*Math.sin(theta)));
    g.drawLine(xPosPix, yPosPix, xPosPix + (int)(lengthFInPix*Math.cos(theta)*Math.sin(theta)), yPosPix + (int)(lengthFInPix*Math.cos(theta)*Math.cos(theta)));

    g.setColor(Color.BLACK);
    g.fillOval(width/2 - 3, height/10 - 3, 6, 6);
    g.drawLine(width/2, height/10, xPosPix, yPosPix);
    g.fillOval(xPosPix - 4, yPosPix - 4, 8, 8);

    g.drawString("T = " + df.format(currentT), 25, 25);

    return g;

}

public void RungeKutta(double theta, double omega){

    a = omegaDot(theta);
    b = omegaDot(theta + 0.5 * h * a);
    c = omegaDot(theta + 0.5 * h * b);
    d = omegaDot(theta + h * c);

    omega = omega + (h / 6) * (a + 2 * b + 2 * c + d);

    a = thetaDot(omega);
    b = thetaDot(omega + 0.5 * h * a);
    c = thetaDot(omega + 0.5 * h * b);
    d = thetaDot(omega + h * c);

    theta = theta + (h / 6) * (a + 2 * b + 2 * c + d);

    this.theta = theta;
    this.omega = omega;

}

public double omegaDot(double theta){

    return -(g / length) * Math.sin(theta);

}

public double thetaDot(double omega){

    return omega;

}
}

Answer 1

Your Runge-Kutta implementation is correctly taken from the scalar case. However, coupled systems have to be treated as vector-valued ODE, the state of the next stage vector depends on all the slopes of the previous stage. You can not separate the updates of the components of the state vector as you tried.

Correct is

aomega = omegaDot(theta);
atheta = thetaDot(omega);
bomega = omegaDot(theta + 0.5 * h * atheta);
btheta = thetaDot(omega + 0.5 * h * aomega);
comega = omegaDot(theta + 0.5 * h * btheta);
ctheta = thetaDot(omega + 0.5 * h * bomega);
domega = omegaDot(theta + h * ctheta);
dtheta = thetaDot(omega + h * comega);

omega = omega + (h / 6) * (aomega + 2 * bomega + 2 * comega + domega);
theta = theta + (h / 6) * (atheta + 2 * btheta + 2 * ctheta + dtheta);

Note how intertwined the slope computations are.

---

One energy function for the pendulum is $$ E(θ,ω)=\frac12ω^2 + \frac{g}{L} (1-\cosθ) $$ At the high point the angular velocity is zero, so $E(θ_H,0)=\frac{g}{L} (1-\cosθ_H)$. At the low point the angle is zero, thus the angular velocity there can be determined from $E(0,ω_L)=\frac12ω_L^2$. The energy is constant, thus $$ ω_L=\sqrt{2\frac{g}{L} (1-\cosθ_H)}=2\sin(θ_H/2)\sqrt{\frac{g}{L}}. $$