[Top][All Lists]
[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]
[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;
}
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] passing parameters in gsl_odeiv,
Tomás Revilla <=