2
$\begingroup$

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

Screenshot of output


Can anybody let me know where it went wrong?

$\endgroup$
6
  • $\begingroup$ 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. $\endgroup$
    – littleO
    Nov 22, 2013 at 1:19
  • $\begingroup$ @littleO I couldn't trace it out. I know matlab, but I don't know how to use graphics. $\endgroup$
    – noufal
    Nov 22, 2013 at 2:12
  • $\begingroup$ 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. $\endgroup$
    – littleO
    Nov 22, 2013 at 7:22
  • $\begingroup$ 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? $\endgroup$ Apr 8, 2014 at 13:06
  • $\begingroup$ @RossMillikan I got it correct. There was a problem with my numerical integration code. $\endgroup$
    – noufal
    Apr 11, 2014 at 9:25

4 Answers 4

2
$\begingroup$

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.

$\endgroup$
3
  • $\begingroup$ Nice, So I need to select initial values for theta and u. Right? $\endgroup$
    – noufal
    Nov 21, 2013 at 11:07
  • $\begingroup$ 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$. $\endgroup$
    – littleO
    Nov 21, 2013 at 11:09
  • $\begingroup$ could you please see my updated question? $\endgroup$
    – noufal
    Nov 21, 2013 at 13:55
1
$\begingroup$

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

$\endgroup$
3
  • $\begingroup$ Thanks Ron, but it seems to be complicated,:-) $\endgroup$
    – noufal
    Nov 21, 2013 at 11:14
  • $\begingroup$ @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. $\endgroup$
    – Ron Gordon
    Nov 21, 2013 at 11:16
  • $\begingroup$ could you please see my question update? $\endgroup$
    – noufal
    Nov 21, 2013 at 13:55
1
$\begingroup$

You have made mistakes in your implementation of the Runge-Kutta Method. Please see the differences.

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);
 kt4=thetadot(u + h * kt3);
 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
double thetadot(double theta, double u)
{
    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;   
     k1=thetadot(theta,u);
     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 );
    glLoadIdentity();
    glOrtho(-ar, ar, -1, 1, -1, 1);

    glMatrixMode( GL_MODELVIEW );
    glLoadIdentity();
//Background
    if(texture_background) {
        glBindTexture( GL_TEXTURE_2D, texture_background );
        glEnable( GL_TEXTURE_2D );
        glBegin( GL_QUADS );
        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;


}

$\endgroup$
1
  • 1
    $\begingroup$ 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*). $\endgroup$
    – hardmath
    Jan 23, 2015 at 1:31
0
$\begingroup$

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

You can also download the source code from https://github.com/sommukhopadhyay/pendulumsimulation

Hope this helps you.

$\endgroup$

You must log in to answer this question.

Not the answer you're looking for? Browse other questions tagged .