diff options
Diffstat (limited to 'gsl-1.9/ode-initval/cstd.c')
| -rw-r--r-- | gsl-1.9/ode-initval/cstd.c | 184 |
1 files changed, 184 insertions, 0 deletions
diff --git a/gsl-1.9/ode-initval/cstd.c b/gsl-1.9/ode-initval/cstd.c new file mode 100644 index 0000000..e4cf4a9 --- /dev/null +++ b/gsl-1.9/ode-initval/cstd.c @@ -0,0 +1,184 @@ +/* ode-initval/cstd.c + * + * Copyright (C) 1996, 1997, 1998, 1999, 2000 Gerard Jungman + * + * 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 2 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. + */ + +#include <config.h> +#include <stdlib.h> +#include <math.h> +#include <gsl/gsl_errno.h> +#include <gsl/gsl_math.h> +#include <gsl/gsl_odeiv.h> + +typedef struct +{ + double eps_abs; + double eps_rel; + double a_y; + double a_dydt; +} +std_control_state_t; + +static void * +std_control_alloc (void) +{ + std_control_state_t * s = + (std_control_state_t *) malloc (sizeof(std_control_state_t)); + + if (s == 0) + { + GSL_ERROR_NULL ("failed to allocate space for std_control_state", + GSL_ENOMEM); + } + + return s; +} + +static int +std_control_init (void * vstate, + double eps_abs, double eps_rel, double a_y, double a_dydt) +{ + std_control_state_t * s = (std_control_state_t *) vstate; + + if (eps_abs < 0) + { + GSL_ERROR ("eps_abs is negative", GSL_EINVAL); + } + else if (eps_rel < 0) + { + GSL_ERROR ("eps_rel is negative", GSL_EINVAL); + } + else if (a_y < 0) + { + GSL_ERROR ("a_y is negative", GSL_EINVAL); + } + else if (a_dydt < 0) + { + GSL_ERROR ("a_dydt is negative", GSL_EINVAL); + } + + s->eps_rel = eps_rel; + s->eps_abs = eps_abs; + s->a_y = a_y; + s->a_dydt = a_dydt; + + return GSL_SUCCESS; +} + +static int +std_control_hadjust(void * vstate, size_t dim, unsigned int ord, const double y[], const double yerr[], const double yp[], double * h) +{ + std_control_state_t *state = (std_control_state_t *) vstate; + + const double eps_abs = state->eps_abs; + const double eps_rel = state->eps_rel; + const double a_y = state->a_y; + const double a_dydt = state->a_dydt; + + const double S = 0.9; + const double h_old = *h; + + double rmax = DBL_MIN; + size_t i; + + for(i=0; i<dim; i++) { + const double D0 = + eps_rel * (a_y * fabs(y[i]) + a_dydt * fabs(h_old * yp[i])) + eps_abs; + const double r = fabs(yerr[i]) / fabs(D0); + rmax = GSL_MAX_DBL(r, rmax); + } + + if(rmax > 1.1) { + /* decrease step, no more than factor of 5, but a fraction S more + than scaling suggests (for better accuracy) */ + double r = S / pow(rmax, 1.0/ord); + + if (r < 0.2) + r = 0.2; + + *h = r * h_old; + + return GSL_ODEIV_HADJ_DEC; + } + else if(rmax < 0.5) { + /* increase step, no more than factor of 5 */ + double r = S / pow(rmax, 1.0/(ord+1.0)); + + if (r > 5.0) + r = 5.0; + + if (r < 1.0) /* don't allow any decrease caused by S<1 */ + r = 1.0; + + *h = r * h_old; + + return GSL_ODEIV_HADJ_INC; + } + else { + /* no change */ + return GSL_ODEIV_HADJ_NIL; + } +} + + +static void +std_control_free (void * vstate) +{ + std_control_state_t *state = (std_control_state_t *) vstate; + free (state); +} + +static const gsl_odeiv_control_type std_control_type = +{"standard", /* name */ + &std_control_alloc, + &std_control_init, + &std_control_hadjust, + &std_control_free}; + +const gsl_odeiv_control_type *gsl_odeiv_control_standard = &std_control_type; + + +gsl_odeiv_control * +gsl_odeiv_control_standard_new(double eps_abs, double eps_rel, + double a_y, double a_dydt) +{ + gsl_odeiv_control * c = + gsl_odeiv_control_alloc (gsl_odeiv_control_standard); + + int status = gsl_odeiv_control_init (c, eps_abs, eps_rel, a_y, a_dydt); + + if (status != GSL_SUCCESS) + { + gsl_odeiv_control_free (c); + GSL_ERROR_NULL ("error trying to initialize control", status); + } + + return c; +} + +gsl_odeiv_control * +gsl_odeiv_control_y_new(double eps_abs, double eps_rel) +{ + return gsl_odeiv_control_standard_new(eps_abs, eps_rel, 1.0, 0.0); +} + + +gsl_odeiv_control * +gsl_odeiv_control_yp_new(double eps_abs, double eps_rel) +{ + return gsl_odeiv_control_standard_new(eps_abs, eps_rel, 0.0, 1.0); +} |