使用GSL双摆解决方案 [英] Double pendulum solution using GSL

问题描述

我学习使用GSL解决ODE。我想用GSL ODE函数来解决双摆问题。我决定用这个公式：

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:

的（来源： HTTP：//www.physics.usyd。 edu.au/~wheat/dpend_html/ ）的

我的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... \n");
        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\nL2: %f\nM1 :%f\nM2:%f\nT_START:%f\nT_END:%f\nS1_ANGLE: %f \nS2_ANGLE: %f\nV1_INIT: %f \nV2_INIT: %f \n",
    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 \n", t, y[0], y[1],y[2],y[3]);
    }


    return 0;
}

这真的不要紧，参数I输入，我总是得到一个错误：

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

GSL：driver.c：354：ERROR：积分限制和/或步骤的方向不
  一致默认GSL错误处理程序调用。

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.

推荐答案

这真的是论点是的顺序。在引用源是在自然命令，为你的混音的原因是不是真的清楚。赋予名称的一些子前pressions的ODE功能也可以写成

It really is the order of arguments in y. In the cited source it was in the "natural" order, the reason for your mix is not really clear. Giving names to some sub-expressions, the ODE function might also be written as

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


    double th1 = y[0], w1 = y[1];
    double th2 = y[2], w2 = y[3];

    f[0] = w1; // dot theta_1 = omega_1
    f[2] = w2; // dot theta_2 = omega_2

    double del = th2 - th1;
    double den = (M1 + M2) - M2 * cos(del) * cos(del);
    double Lwws1 = L1 * (w1*w1) * sin(del);
    double Lwws2 = L2 * (w2*w2) * sin(del);
    double Gs1 = G*sin(th1), Gs2 = G*sin(th2);

    f[1] = (M2 * (Lwws1 + Gs2) * cos(del) + M2 * Lwws2 - (M1 + M2) * Gs1) / (L1*den);

    f[3] = (-M2 * Lwws2 * cos(del) + (M1 + M2) * ( Gs1 * cos(del) - Lwws1 - Gs2) / (L2*den);

    return GSL_SUCCESS;
}

当然，这假定初始点向量被定义为

Of course, this assumes that the initial point vector is defined as

    double y[4] = {S1_ANGLE,V1_INIT,S2_ANGLE,V2_INIT};

和该结果的帧间pretation对应于

and that the interpretation of the results corresponds to that.

这篇关于使用GSL双摆解决方案的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！