|
Packit |
67cb25 |
/* ode-initval/evolve.c
|
|
Packit |
67cb25 |
*
|
|
Packit |
67cb25 |
* Copyright (C) 1996, 1997, 1998, 1999, 2000 Gerard Jungman
|
|
Packit |
67cb25 |
*
|
|
Packit |
67cb25 |
* This program is free software; you can redistribute it and/or modify
|
|
Packit |
67cb25 |
* it under the terms of the GNU General Public License as published by
|
|
Packit |
67cb25 |
* the Free Software Foundation; either version 3 of the License, or (at
|
|
Packit |
67cb25 |
* your option) any later version.
|
|
Packit |
67cb25 |
*
|
|
Packit |
67cb25 |
* This program is distributed in the hope that it will be useful, but
|
|
Packit |
67cb25 |
* WITHOUT ANY WARRANTY; without even the implied warranty of
|
|
Packit |
67cb25 |
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
|
|
Packit |
67cb25 |
* General Public License for more details.
|
|
Packit |
67cb25 |
*
|
|
Packit |
67cb25 |
* You should have received a copy of the GNU General Public License
|
|
Packit |
67cb25 |
* along with this program; if not, write to the Free Software
|
|
Packit |
67cb25 |
* Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
|
|
Packit |
67cb25 |
*/
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* Author: G. Jungman
|
|
Packit |
67cb25 |
*/
|
|
Packit |
67cb25 |
#include <config.h>
|
|
Packit |
67cb25 |
#include <string.h>
|
|
Packit |
67cb25 |
#include <stdlib.h>
|
|
Packit |
67cb25 |
#include <gsl/gsl_math.h>
|
|
Packit |
67cb25 |
#include <gsl/gsl_errno.h>
|
|
Packit |
67cb25 |
#include <gsl/gsl_odeiv.h>
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
#include "odeiv_util.h"
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
gsl_odeiv_evolve *
|
|
Packit |
67cb25 |
gsl_odeiv_evolve_alloc (size_t dim)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
gsl_odeiv_evolve *e =
|
|
Packit |
67cb25 |
(gsl_odeiv_evolve *) malloc (sizeof (gsl_odeiv_evolve));
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
if (e == 0)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
GSL_ERROR_NULL ("failed to allocate space for evolve struct",
|
|
Packit |
67cb25 |
GSL_ENOMEM);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
e->y0 = (double *) malloc (dim * sizeof (double));
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
if (e->y0 == 0)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
free (e);
|
|
Packit |
67cb25 |
GSL_ERROR_NULL ("failed to allocate space for y0", GSL_ENOMEM);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
e->yerr = (double *) malloc (dim * sizeof (double));
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
if (e->yerr == 0)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
free (e->y0);
|
|
Packit |
67cb25 |
free (e);
|
|
Packit |
67cb25 |
GSL_ERROR_NULL ("failed to allocate space for yerr", GSL_ENOMEM);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
e->dydt_in = (double *) malloc (dim * sizeof (double));
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
if (e->dydt_in == 0)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
free (e->yerr);
|
|
Packit |
67cb25 |
free (e->y0);
|
|
Packit |
67cb25 |
free (e);
|
|
Packit |
67cb25 |
GSL_ERROR_NULL ("failed to allocate space for dydt_in", GSL_ENOMEM);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
e->dydt_out = (double *) malloc (dim * sizeof (double));
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
if (e->dydt_out == 0)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
free (e->dydt_in);
|
|
Packit |
67cb25 |
free (e->yerr);
|
|
Packit |
67cb25 |
free (e->y0);
|
|
Packit |
67cb25 |
free (e);
|
|
Packit |
67cb25 |
GSL_ERROR_NULL ("failed to allocate space for dydt_out", GSL_ENOMEM);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
e->dimension = dim;
|
|
Packit |
67cb25 |
e->count = 0;
|
|
Packit |
67cb25 |
e->failed_steps = 0;
|
|
Packit |
67cb25 |
e->last_step = 0.0;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
return e;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
int
|
|
Packit |
67cb25 |
gsl_odeiv_evolve_reset (gsl_odeiv_evolve * e)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
e->count = 0;
|
|
Packit |
67cb25 |
e->failed_steps = 0;
|
|
Packit |
67cb25 |
e->last_step = 0.0;
|
|
Packit |
67cb25 |
return GSL_SUCCESS;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
void
|
|
Packit |
67cb25 |
gsl_odeiv_evolve_free (gsl_odeiv_evolve * e)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
RETURN_IF_NULL (e);
|
|
Packit |
67cb25 |
free (e->dydt_out);
|
|
Packit |
67cb25 |
free (e->dydt_in);
|
|
Packit |
67cb25 |
free (e->yerr);
|
|
Packit |
67cb25 |
free (e->y0);
|
|
Packit |
67cb25 |
free (e);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* Evolution framework method.
|
|
Packit |
67cb25 |
*
|
|
Packit |
67cb25 |
* Uses an adaptive step control object
|
|
Packit |
67cb25 |
*/
|
|
Packit |
67cb25 |
int
|
|
Packit |
67cb25 |
gsl_odeiv_evolve_apply (gsl_odeiv_evolve * e,
|
|
Packit |
67cb25 |
gsl_odeiv_control * con,
|
|
Packit |
67cb25 |
gsl_odeiv_step * step,
|
|
Packit |
67cb25 |
const gsl_odeiv_system * dydt,
|
|
Packit |
67cb25 |
double *t, double t1, double *h, double y[])
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
const double t0 = *t;
|
|
Packit |
67cb25 |
double h0 = *h;
|
|
Packit |
67cb25 |
int step_status;
|
|
Packit |
67cb25 |
int final_step = 0;
|
|
Packit |
67cb25 |
double dt = t1 - t0; /* remaining time, possibly less than h */
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
if (e->dimension != step->dimension)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
GSL_ERROR ("step dimension must match evolution size", GSL_EINVAL);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
if ((dt < 0.0 && h0 > 0.0) || (dt > 0.0 && h0 < 0.0))
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
GSL_ERROR ("step direction must match interval direction", GSL_EINVAL);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* No need to copy if we cannot control the step size. */
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
if (con != NULL)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
DBL_MEMCPY (e->y0, y, e->dimension);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* Calculate initial dydt once if the method can benefit. */
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
if (step->type->can_use_dydt_in)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
int status = GSL_ODEIV_FN_EVAL (dydt, t0, y, e->dydt_in);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
if (status)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
return status;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
try_step:
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
if ((dt >= 0.0 && h0 > dt) || (dt < 0.0 && h0 < dt))
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
h0 = dt;
|
|
Packit |
67cb25 |
final_step = 1;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
else
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
final_step = 0;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
if (step->type->can_use_dydt_in)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
step_status =
|
|
Packit |
67cb25 |
gsl_odeiv_step_apply (step, t0, h0, y, e->yerr, e->dydt_in,
|
|
Packit |
67cb25 |
e->dydt_out, dydt);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
else
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
step_status =
|
|
Packit |
67cb25 |
gsl_odeiv_step_apply (step, t0, h0, y, e->yerr, NULL, e->dydt_out,
|
|
Packit |
67cb25 |
dydt);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* Check for stepper internal failure */
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
if (step_status != GSL_SUCCESS)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
*h = h0; /* notify user of step-size which caused the failure */
|
|
Packit |
67cb25 |
*t = t0; /* restore original t value */
|
|
Packit |
67cb25 |
return step_status;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
e->count++;
|
|
Packit |
67cb25 |
e->last_step = h0;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
if (final_step)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
*t = t1;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
else
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
*t = t0 + h0;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
if (con != NULL)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
/* Check error and attempt to adjust the step. */
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
double h_old = h0;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
const int hadjust_status
|
|
Packit |
67cb25 |
= gsl_odeiv_control_hadjust (con, step, y, e->yerr, e->dydt_out, &h0;;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
if (hadjust_status == GSL_ODEIV_HADJ_DEC)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
/* Check that the reported status is correct (i.e. an actual
|
|
Packit |
67cb25 |
decrease in h0 occured) and the suggested h0 will change
|
|
Packit |
67cb25 |
the time by at least 1 ulp */
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
double t_curr = GSL_COERCE_DBL(*t);
|
|
Packit |
67cb25 |
double t_next = GSL_COERCE_DBL((*t) + h0);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
if (fabs(h0) < fabs(h_old) && t_next != t_curr)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
/* Step was decreased. Undo step, and try again with new h0. */
|
|
Packit |
67cb25 |
DBL_MEMCPY (y, e->y0, dydt->dimension);
|
|
Packit |
67cb25 |
e->failed_steps++;
|
|
Packit |
67cb25 |
goto try_step;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
else
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
h0 = h_old; /* keep current step size */
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
*h = h0; /* suggest step size for next time-step */
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
return step_status;
|
|
Packit |
67cb25 |
}
|