How do I simulate a simple pendulum?

Question

I have the equation of motion of a simple pendulum as $$\frac{d^2\theta}{dt^2} + \frac{g}{l}\sin \theta = 0$$ It's a second order equation. I am trying to simulate it using a SDL library in C++. I know how to solve first order differential equation using Runge-Kutta method. But I can't combine all these. Can anybody help me to solve the differential equation to get the correct simulation?

Update
BTW, I am not sure whether it should be posted here or in any stack-exchange forums. But If you could help me or redirect me to any other forums, that would be helpful.

I have simulated it using Runge-Kutta method and I am adding my code below.

//simplified equations 1
double thetadot(double u)
{
    return u;
}

//simplified equations 2
double udot(double theta)
{
    return (-g / l) * sin(theta);
}


int main(int argc, char *argv[])
{
    double theta, thetanext, u, unext, ku1, ku2, ku3, ku4, kt1, kt2, kt3, kt4;
    if (SDL_Init(SDL_INIT_VIDEO) != 0)
        return 1;

    atexit(SDL_Quit);
    SDL_Surface *screen = SDL_SetVideoMode(width, height, 0, SDL_DOUBLEBUF);
    if (screen == NULL)
        return 2;
    //putting inital values to the function
    u = u0;
    theta = theta0;

    while(true)
    {
        SDL_Event event;
        while(SDL_PollEvent(&event))
        {
            if(event.type == SDL_QUIT)
                return 0;
        }

        double x = xoffset + l * sin(theta);
        double y = yoffset + l * cos(theta);

        SDL_LockSurface(screen);

        //string hanging position
        draw_circle(screen, xoffset, yoffset, 10, 0x0000ff00);
        fill_circle(screen, xoffset, yoffset, 10, 0x0000ff00);

        //draw string
        draw_line(screen, xoffset, yoffset, x, y, 0xff3366ff);

        //draw bob's current position
        fill_circle(screen, (int)x, (int)y, r, 0xff004400);
        draw_circle(screen, (int)x, (int)y, r, 0xff3366ff);

        SDL_Delay(300);
        //SDL_FreeSurface(screen);
        SDL_Flip(screen);

        //Numerical integration of equation 1
        kt1 = thetadot(u);
        kt2 = thetadot(u) + 0.5 * h * kt1;
        kt3 = thetadot(u) + 0.5 * h * kt2;
        kt4 = thetadot(u) + h * kt3;
        thetanext = thetadot(u) + (h / 6) * (kt1 + 2 * kt2 + 2 * kt3 + kt4);

        //Numerical integration of equation 2
        ku1 = udot(theta);
        ku2 = udot(theta) + 0.5 * h * ku1;
        ku3 = udot(theta) + 0.5 * h * ku2;
        ku4 = udot(theta) + h * ku3;
        unext = udot(theta) + (h / 6) * (ku1 + 2 * ku2 + 2 * ku3 + ku4);

        //updating values
        u = unext;
        theta = thetanext;
    }
    return 0;
}

And the output is coming as follows

Can anybody let me know where it went wrong?

Answer 1

You can convert your equation to a first order system by introducing a new function $u = \theta'$. The system satisfied by $u$ and $\theta$ is:

\begin{align*} \theta' &= u \\ u' &= - \frac{g}{\ell} \sin \theta. \end{align*}

You can solve this first order system using Runge-Kutta.

Answer 2

Multiply both sides of the equation by $d\theta/dt$ to get

$$\frac{d\theta}{dt} \frac{d^2 \theta}{dt^2} + \frac{g}{\ell} \frac{d\theta}{dt} \sin{\theta} = 0$$

or

$$\frac12 \frac{d}{dt} \left ( \frac{d\theta}{dt}\right )^2-\frac{g}{\ell} \frac{d}{dt} \cos{\theta} = 0 $$

Therefore

$$\left ( \frac{d\theta}{dt}\right )^2 - \frac{2 g}{\ell} \cos{\theta} = C$$

where $C$ is a constant of integration. This equation is integrable:

$$\int_0^{\theta} \frac{d\theta'}{\sqrt{C+(2 g/\ell) \cos{\theta'}}} = t$$

assuming that $\theta(0)=0$. The integral on the LHS is expressible in terms of a elliptic integral, which must then be inverted to get $\theta(t)$.