c - Double pendulum solution using GSL

Question

I am learning to use GSL to solve ODE. I wanted to solve double pendulum problem using GSL ODE functions. I decided to use this equations:

(source: http://www.physics.usyd.edu.au/~wheat/dpend_html/)

My code:

#include <iostream>
#include <cmath>
#include "gsl/gsl_errno.h"
#include "gsl/gsl_matrix.h"
#include "gsl/gsl_odeiv2.h"
#include "constants.h"

double L1;
double L2;
double M1;
double M2;
double T_START;
double T_END;
double S1_ANGLE;
double S2_ANGLE;
double V1_INIT;
double V2_INIT;

int func(double t, const double y[], double f[], void *params) {

    /*
     * y[0] = theta_2
     * y[1] = omega_1
     * y[2] = theta_1
     * y[3] = omega_2
     */

    f[0] = y[1];
    double del = y[2] - y[1];
    double den1 = (M1 + M2) * L1 - M2 * L1 * cos(del) * cos(del);
    f[1] = (M2 * L1 * y[1] * y[1] * sin(del) * cos(del)
            + M2 * G * sin(y[2]) * cos(del) + M2 * L2 * y[3] * y[3] * sin(del)
            - (M1 + M2) * G * sin(y[0])) / den1;

    f[2] = y[3];
    double den2 = (L2 / L1) * den1;
    f[3] = (-M2 * L2 * y[3] * y[3] * sin(del) * cos(del)
            + (M1 + M2) * G * sin(y[0]) * cos(del)
            - (M1 + M2) * L1 * y[1] * y[1] * sin(del)
            - (M1 + M2) * G * sin(y[2])) / den2;

    return GSL_SUCCESS;
}


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

    /*
     * Arguments list:
     * 1 - length of pendulum 1
     * 2 - length of pendulum 2
     * 3 - mass of pendulum 1
     * 4 - mass of pendulum 2
     * 5 - start time (seconds)
     * 6 - end time (seconds)
     * 7 - initial angle of 1 pendulum (degrees)
     * 8 - initial angle od 2 pendulum
     * 9 - initial angular velocity of 1 pendulum (deegres per second)
     * 10 - initial angular velocity of 2 pendulum
     */

    if (argc != 11) {
        printf("Wrong number of arguments... 
");
        exit(1);
    }

    //Attribution of arguments
    L1 = atof(argv[1]);
    L2 = atof(argv[2]);
    M1 = atof(argv[3]);
    M2 = atof(argv[4]);
    T_START = atof(argv[5]);
    T_END = atof(argv[6]);
    S1_ANGLE = atof(argv[7]);
    S2_ANGLE = atof(argv[8]);
    V1_INIT = atof(argv[9]);
    V2_INIT = atof(argv[10]);

    //converting to radians
    S1_ANGLE=S1_ANGLE*PI/180.0;
    S2_ANGLE=S2_ANGLE*PI/180.0;
    V1_INIT=V1_INIT*PI/180.0;
    V2_INIT=V2_INIT*PI/180.0;
    printf("L1:%f
L2: %f
M1 :%f
M2:%f
T_START:%f
T_END:%f
S1_ANGLE: %f 
S2_ANGLE: %f
V1_INIT: %f 
V2_INIT: %f 
",
    L1,L2,M1,M2,T_START,T_END,S1_ANGLE,S2_ANGLE,V1_INIT,V2_INIT);


    gsl_odeiv2_system sys = {func, NULL, 4, NULL};
    gsl_odeiv2_driver *d =
            gsl_odeiv2_driver_alloc_y_new(&sys, gsl_odeiv2_step_rk4, 1e-6, 1e-6, 0.0);



    double y[4] = {S2_ANGLE,V1_INIT,S1_ANGLE,V2_INIT};
    double t = T_START;
    for (int i = 1; i <= 100; i++) {
        double ti = i * (T_END - T_START) / 100.0;
        int status = gsl_odeiv2_driver_apply(d, &t, ti, y);
        printf("%.5e %.5e %.5e %.5e %.5e 
", t, y[0], y[1],y[2],y[3]);
    }


    return 0;
}

It really does not matter what parameters I input, I always get an error:

gsl: driver.c:354: ERROR: integration limits and/or step direction not consistent Default GSL error handler invoked.

I does not understand what causes my error.

Answer 1

You could also use C++ with odeint library if you are interested. For the double pendulum system. For matrices I use Eigen and for solving ODEs I use odeint , therefore this is the code for your problem.

#include <iostream>
#include <cmath>
#include <Eigen/Dense>
#include <boost/math/constants/constants.hpp>
#include <boost/numeric/odeint.hpp>
#include <iomanip>

using namespace std;
using namespace boost::numeric::odeint;


typedef std::vector< double > state_type;


void equations(const state_type &y, state_type &dy, double x)
{
    const double m1(0.5), m2(0.5),
                 L1(0.1), L2(0.1),
                 g(9.81);

    Eigen::MatrixXd M(2, 2), C(2, 1), B(2,1);

    /*
      Theta 1 =  y[0]
     dTheta 1 =  y[1] = dy[0]
    ddTheta 1 = dy[1]

      Theta 2 =  y[2]
     dTheta 2 =  y[3] = dy[2]
    ddTheta 2 = dy[3]
    */
    double delta(y[0] - y[2]);

    M <<      (m1 + m2)*L1, m2*L2*cos(delta),
          m2*L1*cos(delta),            m2*L2;


    C <<  m2*L1*L2*y[3]*y[3]*sin(delta) + g*(m1 + m2)*sin(y[0]),
         -m2*L1*y[1]*y[1]*sin(delta)    +        m2*g*sin(y[2]);


    //#####################( ODE Equations )################################
    dy[0] = y[1];
    dy[2] = y[3];

    B = M.inverse() * (-C);


    dy[1] = B(0,0);
    dy[3] = B(1,0);

}


int main(int argc, char **argv)
{
    const double dt = 0.01;
    runge_kutta_dopri5 < state_type > stepper;
    double pi = boost::math::constants::pi<double>();

    state_type y(4); 
    // initial values
    y[0] = pi/3.0;  //  Theta 1 
    y[1] = 0.0;     // dTheta 1
    y[2] = pi/4.0; //   Theta 2
    y[3] = 0.0;    //  dTheta 2

    for (double t(0.0); t <= 5; t += dt){

        cout << fixed << setprecision(2);
        std::cout << "t: " << t << "  Theta 1: " << y[0] << "   dTheta 1: " << y[1]
                  << "  Theta 2: " << y[2] << "  dTheta 2: " << y[3] << std::endl;


        stepper.do_step(equations, y, t, dt);
    }

}

The results are