# How do I simulate a simple pendulum?

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
{
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
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?

• Hmm, I haven't read the code carefully, but this should definitely work so most likely there's some straightforward bug. Have you tried using just Euler's method? That should give reasonable results, though not as accurate as Runge-Kutta. I'd get that working first. Also have you considered doing this first in Matlab or python? It's often easier to debug code in one of those languages. The Matlab plotting functions are convenient. Nov 22, 2013 at 1:19
• @littleO I couldn't trace it out. I know matlab, but I don't know how to use graphics. Nov 22, 2013 at 2:12
• What if you just plot $\theta$ as a function of $t$? Does that look correct? I can't tell if the problem is just with your graphics, or if the solution is actually incorrect. Nov 22, 2013 at 7:22
• From your display I can't tell what you think is wrong. Do you expect it to come exactly back to the bottom, or does it only go to one side, or what? Apr 8, 2014 at 13:06
• @RossMillikan I got it correct. There was a problem with my numerical integration code. Apr 11, 2014 at 9:25

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.

• Nice, So I need to select initial values for theta and u. Right? Nov 21, 2013 at 11:07
• Yeah, that's right. To get a unique solution, we normally specify the values $\theta(t_0)$ and $\theta'(t_0)$ for some $t_0$. Nov 21, 2013 at 11:09
• could you please see my updated question? Nov 21, 2013 at 13:55

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

• Thanks Ron, but it seems to be complicated,:-) Nov 21, 2013 at 11:14
• @noufal: Perhaps, but 1) this leads to an analytical form, 2) this leads to a potentially simple Runge-Kutta scheme if you wish to evaluate the integral numerically, and 3) you can also taylor expand the cosine in the integrand to get a useful, analytical approximation out to however many terms you wish. This is all well known. Nov 21, 2013 at 11:16
• could you please see my question update? Nov 21, 2013 at 13:55

Yours:

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


Corrected:

 kt1=thetadot(u);
kt2=thetadot(u + 0.5 * h * kt1);
kt3=thetadot(u + 0.5 * h * kt2);
thetanext = theta + (h / 6.0) * (k1 + 2 * k2 + 2 * k3 + k4);


Similarly for the Udot. Yours:

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


Corrected:

 ku1=udot(theta);
ku2=udot(theta + 0.5 * h * k1);
ku3=udot(theta + 0.5 * h * k2);
ku4=udot(theta + h * k3);
unext = u + (h / 6.0) * (k1 + 2 * k2 + 2 * k3 + k4);


I don't know how to use SDL, but I evaluated this using OpenGL's GLUT libraries on a linux box. I called the file pendulum.cc and ran the following in a terminal

$g++ -o app pendulum.cc -lglut -lGLU -lGL$  ./app


Here is the full source:

#include <iostream>
#include <math.h>       /* sin */
#include <GL/glut.h>
#include <fstream>
const int g=30;
const double l = 1.5;
double u = 4;
double theta = 1;
const double h = 0.01;
double translatex=0;
double translatey=0;
GLuint texture_background;
using namespace std;

//simplified equations 1
{
return u;
}

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

/*---------------------------------------------*/
/* Runge kutta for the theta dot:              */
/* where thetadot = u                      */
/*                                             */
/*---------------------------------------------*/
double RKt (double u, double theta)
{
double thetanext, k1, k2, k3, k4;
k2=thetadot(theta+ 0.5 * h * k1, u + 0.5 * h * k1);
k3=thetadot(theta+ 0.5 * h * k2, u + 0.5 * h * k2);
k4=thetadot(theta+h * k3, u + h * k3);
thetanext = theta + (h / 6.0) * (k1 + 2 * k2 + 2 * k3 + k4);
return thetanext;
}

/*---------------------------------------------*/
/* Runge kutta for the u     dot:              */
/* where udot = d theta/dt = -g/l sin(theta)   */
/*                                             */
/*---------------------------------------------*/
double RKu (double u, double theta)
{
double unext, k1, k2, k3, k4;
k1=udot(theta,u);
k2=udot(theta + 0.5 * h * k1, u + 0.5 * h * k1);
k3=udot(theta + 0.5 * h * k2, u + 0.5 * h * k2);
k4=udot(theta + h * k3, u + h * k3);
unext = u + (h / 6.0) * (k1 + 2 * k2 + 2 * k3 + k4);
return unext;
}

static void display(void)
{

const int width = glutGet(GLUT_WINDOW_WIDTH);
const int height = glutGet(GLUT_WINDOW_HEIGHT);
const float ar = ( float ) width / ( float ) height;

glClear( GL_COLOR_BUFFER_BIT | GL_DEPTH_BUFFER_BIT );

glViewport( 0, 0, width, height );
glMatrixMode( GL_PROJECTION );
glOrtho(-ar, ar, -1, 1, -1, 1);

glMatrixMode( GL_MODELVIEW );
//Background
if(texture_background) {
glBindTexture( GL_TEXTURE_2D, texture_background );
glEnable( GL_TEXTURE_2D );
glTexCoord2f( 1.0, 1.0 );
glVertex2f( -1.0, 1.0 );
glTexCoord2f( 0.0, 1.0 );
glVertex2f( 1.0, 1.0 );
glTexCoord2f( 0.0, 0.0 );
glVertex2f( 1.0, -1.0 );
glTexCoord2f( 1.0, 0.0 );
glVertex2f( -1.0, -1.0 );
glEnd();
glDisable( GL_TEXTURE_2D );
}

//TEXTURE 1

glPushMatrix();

glTranslatef(translatex, translatey, 0);
glBegin( GL_TRIANGLE_FAN );

/*---------------------------------------------*/
/*                                             */
/* This is the only part i messed with for the */
/* size of the pendulum box.                   */
/*---------------------------------------------*/

glVertex2f( 0.1f, 0.1f );  //center
glVertex2f( -0.1f, 0.1f );  //center
glVertex2f( -0.1f, -0.1f );  //center
glVertex2f( 0.1f, -0.1f );  //center

glEnd(  );
glPopMatrix(  );
glDisable( GL_TEXTURE_2D );

glutSwapBuffers();

}

/*---------------------------------------------*/
/*                                             */
/* This is how opengl takes keyboard commands  */
/*      I don't need this, but i could \        */
/*      use this to change the speed           */
/*---------------------------------------------*/
static void key(
unsigned char key,
int a,
int b )
{
glutPostRedisplay();
}

/*---------------------------------------------*/
/*                                             */
/* This is the function that is run between  */
/*  in the loop. I put all the numerical        */
/*   estimatations here.                        */
/*---------------------------------------------*/
static void idle(
void )
{
double oldtheta = theta;
u = RKu(u,theta);
theta = RKt(u,theta);
translatex=l*sin(theta)-sin(oldtheta);
translatey=-l*cos(theta)+cos(oldtheta);
// cout<<translatex;
glutPostRedisplay(  );
}

int main(int argc, char *argv[])
{

glutInit( &argc, argv );

glutInitWindowSize( 640, 640 );
glutInitWindowPosition( 50, 50 );
pre/* glutInitDisplayMode must be called before calling glutCreateWindow
* GLUT, like OpenGL is stateful */
glutInitDisplayMode( GLUT_DOUBLE | GLUT_RGBA | GLUT_DEPTH );
glutCreateWindow( "Pendulum" );

glutDisplayFunc( display );
glutKeyboardFunc( key );
glutIdleFunc( idle );

glutMainLoop(  );
return EXIT_SUCCESS;

}


• Note that while you corrected the arguments passed to thetadot for RK (resp. to udot), the variables in the final evaluation of thetanext (resp. of unext) have been misnamed k* rather than kt* (resp. k* rather than ku*). Jan 23, 2015 at 1:31

I have done it in Android. Please have a look at http://som-itsolutions.blogspot.in/2012/06/android-graphics-and-animation-pendulum.html