GSL (GNU Scientific Library)

Hace ya un tiempo que estoy utilizando la librería GSL y la verdad que me ha parecido una librería fantástica con una buena documentación. Aún así, sobre las ODE voy a dar un par de apuntes.

Código algo más legible

Vamos a tomar el mismo ejemplo que en la documentación, un oscilador de Van der Pol con la ecuación

Para hacer el código más legible lo haremos más cercano a nuestras ecuaciones

int
func (double t,
      const double y[],
      double yp[],
      void* params)
{
  (void)(t);  /* Evitar warning no usado */
  double mu = *(double*)params;

  double u = y[0];
  double v = y[1];

  double ut = v;
  double vt = -u + mu*v*(1 - u*u);

  yp[0] = ut;
  yp[1] = vt;

  return GSL_SUCCESS;
}

Parece una tontería… pero viendo la definición ‘ut’ y ‘vt’ vemos claramente que es exactamente nuestro sistema de ecuaciones. Déjate de arreglos con índices y usa las variables de tus ecuaciones.

Es habitual referirnos a las diferentes derivadas (respecto del tiempo) como: ‘ut’ para la primera derivada, ‘utt’ para la segunda, etc Si fuera respecto de ‘x’ el aspecto sería: ‘ux’, ‘uxx’, etc… También es habitual encontrar ‘up’, ‘upp’, etc para ‘u prima’, ‘u prima prima’, etc…

Por cierto… aunque solo se usa un parámetro, me resulta perturbador el uso que se hace de ‘params’. Yo personalmente prefiero utilizar el método genérico para el caso de tener varios parámetros y así utilizar siempre la misma forma de uso homogeneizando el código. Pero esto ya es un gusto personal.

¿Qué hace la jacobiana en el ejemplo?

Pues… nada… porque en el ejemplo se usa un método explícito como es ‘gsl_odeiv2_step_rk8pd’ por tanto olvídate de la jacobiana si haces uso de un método explícito.

¿Qué pasa si usamos un método implícito?

Esta vez sí, tendrás que hacer uso de la jacobiana pero si usas el driver no te lies con

…This stepper requires the access to the driver object via gsl_odeiv2_step_set_driver().

No hace falta que crees un ‘stepper’ (después lo tendrías que borrar tú) y se lo pases al ‘driver’. El ‘driver’ tiene acceso a todo lo necesario, así que podemos actuar de la forma habitual

...
T = gsl_odeiv2_step_msadams;
d = gsl_odeiv2_driver_alloc_y_new(&sys, T, hInitial, epsAbs, epsRel);
...
status = gsl_odeiv2_driver_apply(d, &t0, ti, y);
...
gsl_odeiv2_driver_free(d);
...