[Help-gsl] passing parameters in gsl_odeiv

[Tomás Revilla]

Hi,

I am currently struggling on how to pass parameters as
vectors or matrices (and structured data too) for
function definition in gsl_odeiv_... The van der Pol
example (only 1 parameter) is too simple, and I
googled the internet for examples using many
parameters with NULL success. Actually, I solved it by
just declaring the parameters globally and setting:
   
gsl_odeiv_system sys = {func, NULL, 3, NULL};

So func() can make use of global parameters. But I
really want to pass the whole par[3][3] matrix from
main() to func() in the following (and more
elaborated) case:

/* Nonlinear 3-species competition
 *
 *   dN1/dt = f1(N1,N2,N3) = N1 * (1 - N1 - alpha*N2 -
beta*N3)
 *   dN2/dt = f2(N1,N2,N3) = N2 * (1 - beta*N1 - N2 -
alpha*N3)
 *   dN3/dt = f3(N1,N2,N3) = N3 * (1 - alpha*N1 -
beta*N2 - N3))
 *
 */

#include <stdio.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_odeiv.h>

#define ALPHA 0.8
#define BETA  1.2
#define N1 1.0
#define N2 0.8
#define N3 0.2

int func(double t, const double y[], double f[], void
*params)
{
        /* HOWTO RECOVER MY par[3][3] MATRIX? */
        /* Which is the rigth way?            */

      par .... ?

        f[0] = y[0]*(1.0 - par[0][0]*y[0] - par[0][1]*y[1] -
par[0][2]*y[2]); 
        f[1] = y[1]*(1.0 - par[1][0]*y[0] - par[1][1]*y[1] -
par[1][2]*y[2]);
        f[2] = y[2]*(1.0 - par[2][0]*y[0] - par[2][1]*y[1] -
par[2][2]*y[2]);
        
        return GSL_SUCCESS;
}

int main(void)
{
      /* These are the parameters */
        double par[3][3] =
        {
                {1.0, ALPHA, BETA},
                {BETA, 1.0, ALPHA},
                {ALPHA, BETA, 1.0}
        };

        const gsl_odeiv_step_type *T = gsl_odeiv_step_rk4;
        gsl_odeiv_step *s = gsl_odeiv_step_alloc(T,3);
        gsl_odeiv_control *c = gsl_odeiv_control_y_new(1e-6,
0.0);
        gsl_odeiv_evolve *e = gsl_odeiv_evolve_alloc(3);

        gsl_odeiv_system sys = {func, NULL, 3, &par};

        double t = 0.0, t_end = 3000.0;
        double h = 1e-6;
        double y[3] = {N1, N2, N3};
        int tt;

        /* Evolution loop, at fixed intervals */
        for(tt = 0; tt <= 3000; tt++)
        {
                double ti = tt * 1.0;
                while(t < ti)
                {
                        int status = gsl_odeiv_evolve_apply
                                (e, c, s, &sys, &t, ti, &h, y);
                        
                        if(status != GSL_SUCCESS)
                                break;
                }
                
                printf("%.5e\t%.5e\t%.5e\t%.5e\n", t, y[0], y[1],
y[2]);
        }

        gsl_odeiv_evolve_free(e);
        gsl_odeiv_control_free(c);
        gsl_odeiv_step_free(s);

        return 0;
}


Thanks in advance,


     Tomas Revilla


PS: I am new with GSL, does anyone have a tutorial or
something, apart from the tech manual?


__________________________________________________
Correo Yahoo!
Espacio para todos tus mensajes, antivirus y antispam ¡gratis! 
Regístrate ya - http://correo.espanol.yahoo.com/