/* ode-initval2/driver.c
*
* Copyright (C) 2009, 2010 Tuomo Keskitalo
*
* This program is free software; you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
* the Free Software Foundation; either version 3 of the License, or (at
* your option) any later version.
*
* This program is distributed in the hope that it will be useful, but
* WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
* General Public License for more details.
*
* You should have received a copy of the GNU General Public License
* along with this program; if not, write to the Free Software
* Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
*/
/* Driver routine for odeiv2. This is a wrapper for low level GSL
functions that allows a simple interface to step, control and
evolve layers.
*/
#include <config.h>
#include <math.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_odeiv2.h>
#include <gsl/gsl_machine.h>
static gsl_odeiv2_driver *
driver_alloc (const gsl_odeiv2_system * sys, const double hstart,
const gsl_odeiv2_step_type * T)
{
/* Allocates and initializes an ODE driver system. Step and evolve
objects are allocated here, but control object is allocated in
another function.
*/
gsl_odeiv2_driver *state;
if (sys == NULL)
{
GSL_ERROR_NULL ("gsl_odeiv2_system must be defined", GSL_EINVAL);
}
state = (gsl_odeiv2_driver *) calloc (1, sizeof (gsl_odeiv2_driver));
if (state == NULL)
{
GSL_ERROR_NULL ("failed to allocate space for driver state",
GSL_ENOMEM);
}
{
const size_t dim = sys->dimension;
if (dim == 0)
{
gsl_odeiv2_driver_free(state);
GSL_ERROR_NULL
("gsl_odeiv2_system dimension must be a positive integer",
GSL_EINVAL);
}
state->sys = sys;
state->s = gsl_odeiv2_step_alloc (T, dim);
if (state->s == NULL)
{
gsl_odeiv2_driver_free(state);
GSL_ERROR_NULL ("failed to allocate step object", GSL_ENOMEM);
}
state->e = gsl_odeiv2_evolve_alloc (dim);
}
if (state->e == NULL)
{
gsl_odeiv2_driver_free(state);
GSL_ERROR_NULL ("failed to allocate evolve object", GSL_ENOMEM);
}
if (hstart > 0.0 || hstart < 0.0)
{
state->h = hstart;
}
else
{
gsl_odeiv2_driver_free(state);
GSL_ERROR_NULL ("invalid hstart", GSL_EINVAL);
}
state->h = hstart;
state->hmin = 0.0;
state->hmax = GSL_DBL_MAX;
state->nmax = 0;
state->n = 0;
state->c = NULL;
return state;
}
int
gsl_odeiv2_driver_set_hmin (gsl_odeiv2_driver * d, const double hmin)
{
/* Sets minimum allowed step size fabs(hmin) for driver. It is
required that hmin <= fabs(h) <= hmax. */
if ((fabs (hmin) > fabs (d->h)) || (fabs (hmin) > d->hmax))
{
GSL_ERROR_NULL ("hmin <= fabs(h) <= hmax required", GSL_EINVAL);
}
d->hmin = fabs (hmin);
return GSL_SUCCESS;
}
int
gsl_odeiv2_driver_set_hmax (gsl_odeiv2_driver * d, const double hmax)
{
/* Sets maximum allowed step size fabs(hmax) for driver. It is
required that hmin <= fabs(h) <= hmax. */
if ((fabs (hmax) < fabs (d->h)) || (fabs (hmax) < d->hmin))
{
GSL_ERROR_NULL ("hmin <= fabs(h) <= hmax required", GSL_EINVAL);
}
if (hmax > 0.0 || hmax < 0.0)
{
d->hmax = fabs (hmax);
}
else
{
GSL_ERROR_NULL ("invalid hmax", GSL_EINVAL);
}
return GSL_SUCCESS;
}
int
gsl_odeiv2_driver_set_nmax (gsl_odeiv2_driver * d,
const unsigned long int nmax)
{
/* Sets maximum number of allowed steps (nmax) for driver */
d->nmax = nmax;
return GSL_SUCCESS;
}
gsl_odeiv2_driver *
gsl_odeiv2_driver_alloc_y_new (const gsl_odeiv2_system * sys,
const gsl_odeiv2_step_type * T,
const double hstart,
const double epsabs, const double epsrel)
{
/* Initializes an ODE driver system with control object of type y_new. */
gsl_odeiv2_driver *state = driver_alloc (sys, hstart, T);
if (state == NULL)
{
GSL_ERROR_NULL ("failed to allocate driver object", GSL_ENOMEM);
}
if (epsabs >= 0.0 && epsrel >= 0.0)
{
state->c = gsl_odeiv2_control_y_new (epsabs, epsrel);
if (state->c == NULL)
{
gsl_odeiv2_driver_free (state);
GSL_ERROR_NULL ("failed to allocate control object", GSL_ENOMEM);
}
}
else
{
gsl_odeiv2_driver_free (state);
GSL_ERROR_NULL ("epsabs and epsrel must be positive", GSL_EINVAL);
}
/* Distribute pointer to driver object */
gsl_odeiv2_step_set_driver (state->s, state);
gsl_odeiv2_evolve_set_driver (state->e, state);
gsl_odeiv2_control_set_driver (state->c, state);
return state;
}
gsl_odeiv2_driver *
gsl_odeiv2_driver_alloc_yp_new (const gsl_odeiv2_system * sys,
const gsl_odeiv2_step_type * T,
const double hstart,
const double epsabs, const double epsrel)
{
/* Initializes an ODE driver system with control object of type yp_new. */
gsl_odeiv2_driver *state = driver_alloc (sys, hstart, T);
if (state == NULL)
{
GSL_ERROR_NULL ("failed to allocate driver object", GSL_ENOMEM);
}
if (epsabs >= 0.0 && epsrel >= 0.0)
{
state->c = gsl_odeiv2_control_yp_new (epsabs, epsrel);
if (state->c == NULL)
{
gsl_odeiv2_driver_free (state);
GSL_ERROR_NULL ("failed to allocate control object", GSL_ENOMEM);
}
}
else
{
gsl_odeiv2_driver_free (state);
GSL_ERROR_NULL ("epsabs and epsrel must be positive", GSL_EINVAL);
}
/* Distribute pointer to driver object */
gsl_odeiv2_step_set_driver (state->s, state);
gsl_odeiv2_evolve_set_driver (state->e, state);
gsl_odeiv2_control_set_driver (state->c, state);
return state;
}
gsl_odeiv2_driver *
gsl_odeiv2_driver_alloc_standard_new (const gsl_odeiv2_system * sys,
const gsl_odeiv2_step_type * T,
const double hstart,
const double epsabs,
const double epsrel, const double a_y,
const double a_dydt)
{
/* Initializes an ODE driver system with control object of type
standard_new.
*/
gsl_odeiv2_driver *state = driver_alloc (sys, hstart, T);
if (state == NULL)
{
GSL_ERROR_NULL ("failed to allocate driver object", GSL_ENOMEM);
}
if (epsabs >= 0.0 && epsrel >= 0.0)
{
state->c =
gsl_odeiv2_control_standard_new (epsabs, epsrel, a_y, a_dydt);
if (state->c == NULL)
{
gsl_odeiv2_driver_free (state);
GSL_ERROR_NULL ("failed to allocate control object", GSL_ENOMEM);
}
}
else
{
gsl_odeiv2_driver_free (state);
GSL_ERROR_NULL ("epsabs and epsrel must be positive", GSL_EINVAL);
}
/* Distribute pointer to driver object */
gsl_odeiv2_step_set_driver (state->s, state);
gsl_odeiv2_evolve_set_driver (state->e, state);
gsl_odeiv2_control_set_driver (state->c, state);
return state;
}
gsl_odeiv2_driver *
gsl_odeiv2_driver_alloc_scaled_new (const gsl_odeiv2_system * sys,
const gsl_odeiv2_step_type * T,
const double hstart,
const double epsabs, const double epsrel,
const double a_y, const double a_dydt,
const double scale_abs[])
{
/* Initializes an ODE driver system with control object of type
scaled_new.
*/
gsl_odeiv2_driver *state = driver_alloc (sys, hstart, T);
if (state == NULL)
{
GSL_ERROR_NULL ("failed to allocate driver object", GSL_ENOMEM);
}
if (epsabs >= 0.0 && epsrel >= 0.0)
{
state->c = gsl_odeiv2_control_scaled_new (epsabs, epsrel, a_y, a_dydt,
scale_abs, sys->dimension);
if (state->c == NULL)
{
gsl_odeiv2_driver_free (state);
GSL_ERROR_NULL ("failed to allocate control object", GSL_ENOMEM);
}
}
else
{
gsl_odeiv2_driver_free (state);
GSL_ERROR_NULL ("epsabs and epsrel must be positive", GSL_EINVAL);
}
/* Distribute pointer to driver object */
gsl_odeiv2_step_set_driver (state->s, state);
gsl_odeiv2_evolve_set_driver (state->e, state);
gsl_odeiv2_control_set_driver (state->c, state);
return state;
}
int
gsl_odeiv2_driver_apply (gsl_odeiv2_driver * d, double *t,
const double t1, double y[])
{
/* Main driver function that evolves the system from t to t1. In
beginning vector y contains the values of dependent variables at
t. This function returns values at t=t1 in y. In case of
unrecoverable error, y and t contains the values after the last
successful step.
*/
int sign = 0;
d->n = 0;
/* Determine integration direction sign */
if (d->h > 0.0)
{
sign = 1;
}
else
{
sign = -1;
}
/* Check that t, t1 and step direction are sensible */
if (sign * (t1 - *t) < 0.0)
{
GSL_ERROR_NULL
("integration limits and/or step direction not consistent",
GSL_EINVAL);
}
/* Evolution loop */
while (sign * (t1 - *t) > 0.0)
{
int s = gsl_odeiv2_evolve_apply (d->e, d->c, d->s, d->sys,
t, t1, &(d->h), y);
if (s != GSL_SUCCESS)
{
return s;
}
/* Check for maximum allowed steps */
if ((d->nmax > 0) && (d->n > d->nmax))
{
return GSL_EMAXITER;
}
/* Set step size if maximum size is exceeded */
if (fabs (d->h) > d->hmax)
{
d->h = sign * d->hmax;
}
/* Check for too small step size */
if (fabs (d->h) < d->hmin)
{
return GSL_ENOPROG;
}
d->n++;
}
return GSL_SUCCESS;
}
int
gsl_odeiv2_driver_apply_fixed_step (gsl_odeiv2_driver * d, double *t,
const double h, const unsigned long int n,
double y[])
{
/* Alternative driver function that evolves the system from t using
* n steps of size h. In the beginning vector y contains the values
* of dependent variables at t. This function returns values at t =
* t + n * h in y. In case of an unrecoverable error, y and t
* contains the values after the last successful step.
*/
unsigned long int i;
d->n = 0;
/* Evolution loop */
for (i = 0; i < n; i++)
{
int s = gsl_odeiv2_evolve_apply_fixed_step (d->e, d->c, d->s, d->sys,
t, h, y);
if (s != GSL_SUCCESS)
{
return s;
}
d->n++;
}
return GSL_SUCCESS;
}
int
gsl_odeiv2_driver_reset (gsl_odeiv2_driver * d)
{
/* Reset the driver. Resets evolve and step objects. */
{
int s = gsl_odeiv2_evolve_reset (d->e);
if (s != GSL_SUCCESS)
{
return s;
}
}
{
int s = gsl_odeiv2_step_reset (d->s);
if (s != GSL_SUCCESS)
{
return s;
}
}
return GSL_SUCCESS;
}
int
gsl_odeiv2_driver_reset_hstart (gsl_odeiv2_driver * d, const double hstart)
{
/* Resets current driver and sets initial step size to hstart */
gsl_odeiv2_driver_reset (d);
if ((d->hmin > fabs (hstart)) || (fabs (hstart) > d->hmax))
{
GSL_ERROR_NULL ("hmin <= fabs(h) <= hmax required", GSL_EINVAL);
}
if (hstart > 0.0 || hstart < 0.0)
{
d->h = hstart;
}
else
{
GSL_ERROR_NULL ("invalid hstart", GSL_EINVAL);
}
return GSL_SUCCESS;
}
void
gsl_odeiv2_driver_free (gsl_odeiv2_driver * state)
{
if (state->c)
gsl_odeiv2_control_free (state->c);
if (state->e)
gsl_odeiv2_evolve_free (state->e);
if (state->s)
gsl_odeiv2_step_free (state->s);
free (state);
}