help-gsl
[Top][All Lists]

## [Help-gsl] passing parameters in gsl_odeiv

 From: Tomás Revilla Subject: [Help-gsl] passing parameters in gsl_odeiv Date: Wed, 15 Jun 2005 08:40:28 -0500 (CDT)

```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;
}

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/

```