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$
  • $\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 '13 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 '13 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 '13 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$ – Ross Millikan Apr 8 '14 at 13:06
  • $\begingroup$ @RossMillikan I got it correct. There was a problem with my numerical integration code. $\endgroup$ – noufal Apr 11 '14 at 9:25
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$
  • $\begingroup$ Nice, So I need to select initial values for theta and u. Right? $\endgroup$ – noufal Nov 21 '13 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 '13 at 11:09
  • $\begingroup$ could you please see my updated question? $\endgroup$ – noufal Nov 21 '13 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$
  • $\begingroup$ Thanks Ron, but it seems to be complicated,:-) $\endgroup$ – noufal Nov 21 '13 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 '13 at 11:16
  • $\begingroup$ could you please see my question update? $\endgroup$ – noufal Nov 21 '13 at 13:55
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$
0
$\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
    $\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 '15 at 1:31

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.