Re: [Help-gsl] passing parameters in gsl_odeiv

[Pau Cervera Badia]

In the following code I'm just introducing a dummy parameter alpha that 
is just 1 (and so I'm not modifying the van der Pol eq.).

It works. I don't know if there's a better solution. I would like to 
know, in that case.

Here it is the code:

-------------------- begining of code ----------------------------------

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

#define SNAPSHOTS 100
#define EPS_ABS 1e-6

struct parametres_{
 double mu;
 double alpha;
};

int
func (double t, const double y[], double f[], void *params)
{
 struct parametres_ parametres = *(struct parametres_ *)params;
 double mu = parametres.mu;
 double alpha = parametres.alpha;
 f[0] = alpha * y[1];
 f[1] = -y[0] - mu * y[1] * (y[0] * y[0] - 1);

return GSL_SUCCESS;
}

int
jac (double t, const double y[], double *dfdy,
    double dfdt[], void * params)
{
 struct parametres_ parametres = *(struct parametres_ *)params;
 double mu = parametres.mu;
 double alpha = parametres.alpha;
 gsl_matrix_view dfdy_mat = 
   gsl_matrix_view_array (dfdy, 2, 2);
 gsl_matrix *m = &dfdy_mat.matrix;
 gsl_matrix_set (m, 0, 0, 0.0);
 gsl_matrix_set (m, 0, 1, 1.0);
 gsl_matrix_set (m, 1, 0, -2.0 * mu * y[0] * y[1] - 1.0);
 gsl_matrix_set (m, 1, 1, - mu * (y[0] * y[0] - 1.0));
 dfdt[0] = 0.0;
 dfdt[1] = 0.0;

 return GSL_SUCCESS;
}

int
main (void)
{
 static double eps_abs = EPS_ABS;
 const gsl_odeiv_step_type *T = gsl_odeiv_step_rk8pd;

 gsl_odeiv_step *s = gsl_odeiv_step_alloc (T, 2);
 gsl_odeiv_control *c = gsl_odeiv_control_y_new(eps_abs, 0.0);
 gsl_odeiv_evolve *e = gsl_odeiv_evolve_alloc (2);

 struct parametres_ parametres;
 parametres.mu = 1.;
 parametres.alpha = 10.;

 gsl_odeiv_system sys = {func, jac, 2, &parametres};

 double t = 0.0, t1 = 1000.0;
 double h = 1e-6;
 double y[2] = {1.0, 0.0};

 int i;
 for (i = 0; i <= SNAPSHOTS; i++) {
   double ti = i * t1/SNAPSHOTS;

   while (t < ti) {
     int status = gsl_odeiv_evolve_apply (e, c, s,
                                           &sys,
                                           &t, t1,
                                           &h, y);


     if (status != GSL_SUCCESS)
        break;
   }

   printf("%.5e %.5e %.5e\n", t, y[0], y[1]);
 }

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

 return 0;
}

--------------------------end of code --------------------------------



Tomás Revilla wrote:

> 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/ 
>
>
> _______________________________________________
> Help-gsl mailing list
> address@hidden
> http://lists.gnu.org/mailman/listinfo/help-gsl
>
>
>  

--
Pau Cervera i Badia       e-mail: address@hidden

***************************************************
Departament de Física Fonamental
Universitat de Barcelona
Martí i Franqués, 1
Planta 3, despatx 346 bis
08028 Barcelona
Spain

tel: +34 934 921 155
***************************************************
"Most of the time I don't have much fun. The rest of
the time I don't have any fun at all." WA