cgsldifferential-equationsnonlinear-functions

How to consult GSL ODE with many parameters and harmonic functions


I'm working on non-linear differential equation using GSL. The thing is I'm quite new on C stuffs. I just adapted the sample on GNU site into the equation I'm interested in right now.

This is the equation:

d2x/dt2 + r*dx/dy + cos(x) + v*cos(2*x+0.4) E1*sin(wt) + E2*sin(2*w*t+a) = 0

What I am stuck is I have no idea how to plug in multiple parameters in the codes. Moreover, I don't know how to employ cosine or sine function in this code.

I tried to figure out this problem, by searching on Google all the way. I couldn't find any thing that helps me.

#include <stdio.h>
#include <gsl/gsl_errno.h>
#include <math.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_odeiv2.h>

int func (double t, const double x[], double y[], void *params)
{
    double r = *(double *)params;
    double v = *(double *)params;
    double w = *(double *)params;
    double E1 = *(double *)params;
    double E2  = *(double *)params;
    double a  = *(double *)params;
    y[0] = x[1];
    y[1] = -r*x[1] - cos(x[0]) - v*cos(2*x[0]+0.4) - E1*sin(w*t) - E2*sin(2*w*t+a);
    return GSL_SUCCESS;
}

int jac (double t, const double x[], double *dydx, double dydt[], void *params)
{
    double r = *(double *)params;
    double v = *(double *)params;
    double w = *(double *)params;
    double E1 = *(double *)params;
    double E2  = *(double *)params;
    double a  = *(double *)params;
    gsl_matrix_view dydx_mat = gsl_matrix_view_array (dydx, 2, 2);
    gsl_matrix * m = &dydx_mat.matrix;
    gsl_matrix_set (m, 0, 0, 0.0);
    gsl_matrix_set (m, 0, 1, 1.0);
    gsl_matrix_set (m, 1, 0, sin(x[0]) + 2*v*sin(2*x[0]+0.4));
    gsl_matrix_set (m, 1, 1, -r);
    dydt[0] = 0.0;
    dydt[1] = 0.0;
    return GSL_SUCCESS;
}

int main (void)
{
    double r = 0.0;
    double v = 0.0;
    double w = 2.4;
    double E1 = -2.3;
    double E2 = 0;
    double a = 0.7;
    gsl_odeiv2_system sys = {func, jac, 2, &r, &v, &w, &E1, &E2, &a};

    gsl_odeiv2_driver *d = gsl_odeiv2_driver_alloc_x_new (&sys, gsl_odeiv2_step_rk8pd, 1e-6, 1e-6, 0.0);
    int i;
    double t = 0.0, t1 = 10000;
    double x[2] = {0.0, 0.0};

    for (i = 1 ; i<=10000; i++)
        {
            double ti = i*t1/10000;
            int status = gsl_odeiv2_driver_apply (d, &t, ti, x);

            if (status != GSL_SUCCESS)
                {
                    printf("error, return value%d\n", status);
                    break;
                }
            printf("%.5e %.5e %.5e\n", t, x[0], x[1]);
        }

    gsl_odeiv2_driver_free (d);
    return 0;
}

Solution

  • The params argument is a pointer (address / memory location) to some arbitrary data structure. In the example from the GSL documentation, their equation contained only one parameter, which means it's okay to just pass the address of a double-precision number.

    However, for your problem, you need to access 6 different parameters. You can't access every parameter with the same address!

    /* this doesn't work! */
    double r  = *(double *)params;
    double v  = *(double *)params;
    double w  = *(double *)params;
    double E1 = *(double *)params;
    double E2 = *(double *)params;
    double a  = *(double *)params;
    

    Since all the addresses are the same, you are referring to the same number. To remedy this, you can either: store all the parameters in an array of length 6, or store them in a predefined data structure. The latter approach is more readable so I will demonstrate that.

    First define a data type to specify what parameters you will store:

    struct param_type {
        double r;
        double v;
        double w;
        double E1;
        double E2;
        double a;
    };
    

    Now, create a structure of this type in the main function and store the actual values of the parameter:

    struct param_type my_params = {r, v, w, E1, E2, a};
    

    When defining the system, you store a pointer to that struct param_type:

    gsl_odeiv2_system sys = {func, jac, 2, &my_params};
    

    To use the parameter inside func and jac, you simply cast the params argument from a generic pointer (void *) to a pointer for your specific data type (struct param_type *):

    struct param_type *my_params_pointer = params;
    

    (Note that in C++ this must be written with an explicit cast.) Finally, you can access the parameters via:

    double r  = my_params_pointer->r;
    double v  = my_params_pointer->v;
    double w  = my_params_pointer->w;
    double E1 = my_params_pointer->E1;
    double E2 = my_params_pointer->E2;
    double a  = my_params_pointer->a;
    

    The arrow -> is used here instead of the dot . because my_params_pointer is a pointer and needs to be dereferenced before use.