diff options
Diffstat (limited to 'gsl-1.9/ode-initval/rk2.c')
| -rw-r--r-- | gsl-1.9/ode-initval/rk2.c | 250 |
1 files changed, 250 insertions, 0 deletions
diff --git a/gsl-1.9/ode-initval/rk2.c b/gsl-1.9/ode-initval/rk2.c new file mode 100644 index 0000000..e7696d1 --- /dev/null +++ b/gsl-1.9/ode-initval/rk2.c @@ -0,0 +1,250 @@ +/* ode-initval/rk2.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. + */ + +/* Runge-Kutta 2(3), Euler-Cauchy */ + +/* Author: G. Jungman + */ + +/* Reference: Abramowitz & Stegun, section 25.5. Runge-Kutta 2nd (25.5.7) + and 3rd (25.5.8) order methods */ + +#include <config.h> +#include <stdlib.h> +#include <string.h> +#include <gsl/gsl_errno.h> +#include <gsl/gsl_odeiv.h> + +#include "odeiv_util.h" + +typedef struct +{ + double *k1; + double *k2; + double *k3; + double *ytmp; +} +rk2_state_t; + +static void * +rk2_alloc (size_t dim) +{ + rk2_state_t *state = (rk2_state_t *) malloc (sizeof (rk2_state_t)); + + if (state == 0) + { + GSL_ERROR_NULL ("failed to allocate space for rk2_state", GSL_ENOMEM); + } + + state->k1 = (double *) malloc (dim * sizeof (double)); + + if (state->k1 == 0) + { + free (state); + GSL_ERROR_NULL ("failed to allocate space for k1", GSL_ENOMEM); + } + + state->k2 = (double *) malloc (dim * sizeof (double)); + + if (state->k2 == 0) + { + free (state->k1); + free (state); + GSL_ERROR_NULL ("failed to allocate space for k2", GSL_ENOMEM); + } + + state->k3 = (double *) malloc (dim * sizeof (double)); + + if (state->k3 == 0) + { + free (state->k2); + free (state->k1); + free (state); + GSL_ERROR_NULL ("failed to allocate space for k2", GSL_ENOMEM); + } + + state->ytmp = (double *) malloc (dim * sizeof (double)); + + if (state->ytmp == 0) + { + free (state->k3); + free (state->k2); + free (state->k1); + free (state); + GSL_ERROR_NULL ("failed to allocate space for k2", GSL_ENOMEM); + } + + return state; +} + + +static int +rk2_apply (void *vstate, + size_t dim, + double t, + double h, + double y[], + double yerr[], + const double dydt_in[], + double dydt_out[], + const gsl_odeiv_system * sys) +{ + rk2_state_t *state = (rk2_state_t *) vstate; + + size_t i; + + double *const k1 = state->k1; + double *const k2 = state->k2; + double *const k3 = state->k3; + double *const ytmp = state->ytmp; + + /* k1 step */ + /* k1 = f(t,y) */ + + if (dydt_in != NULL) + { + DBL_MEMCPY (k1, dydt_in, dim); + } + else + { + int s = GSL_ODEIV_FN_EVAL (sys, t, y, k1); + + if (s != GSL_SUCCESS) + { + return s; + } + } + + /* k2 step */ + /* k2 = f(t + 0.5*h, y + 0.5*k1) */ + + for (i = 0; i < dim; i++) + { + ytmp[i] = y[i] + 0.5 * h * k1[i]; + } + + { + int s = GSL_ODEIV_FN_EVAL (sys, t + 0.5 * h, ytmp, k2); + + if (s != GSL_SUCCESS) + { + return s; + } + } + + /* k3 step */ + /* for 3rd order estimates, is used for error estimation + k3 = f(t + h, y - k1 + 2*k2) */ + + for (i = 0; i < dim; i++) + { + ytmp[i] = y[i] + h * (-k1[i] + 2.0 * k2[i]); + } + + { + int s = GSL_ODEIV_FN_EVAL (sys, t + h, ytmp, k3); + + if (s != GSL_SUCCESS) + { + return s; + } + } + + /* final sum */ + + for (i = 0; i < dim; i++) + { + /* Save original values if derivative evaluation below fails */ + ytmp[i] = y[i]; + + { + const double ksum3 = (k1[i] + 4.0 * k2[i] + k3[i]) / 6.0; + y[i] += h * ksum3; + } + } + + /* Derivatives at output */ + + if (dydt_out != NULL) + { + int s = GSL_ODEIV_FN_EVAL (sys, t + h, y, dydt_out); + + if (s != GSL_SUCCESS) + { + /* Restore original values */ + DBL_MEMCPY (y, ytmp, dim); + + return s; + } + } + + /* Error estimation */ + + for (i = 0; i < dim; i++) + { + const double ksum3 = (k1[i] + 4.0 * k2[i] + k3[i]) / 6.0; + yerr[i] = h * (k2[i] - ksum3); + } + + return GSL_SUCCESS; +} + +static int +rk2_reset (void *vstate, size_t dim) +{ + rk2_state_t *state = (rk2_state_t *) vstate; + + DBL_ZERO_MEMSET (state->k1, dim); + DBL_ZERO_MEMSET (state->k2, dim); + DBL_ZERO_MEMSET (state->k3, dim); + DBL_ZERO_MEMSET (state->ytmp, dim); + + return GSL_SUCCESS; +} + +static unsigned int +rk2_order (void *vstate) +{ + rk2_state_t *state = (rk2_state_t *) vstate; + state = 0; /* prevent warnings about unused parameters */ + return 2; +} + +static void +rk2_free (void *vstate) +{ + rk2_state_t *state = (rk2_state_t *) vstate; + free (state->k1); + free (state->k2); + free (state->k3); + free (state->ytmp); + free (state); +} + +static const gsl_odeiv_step_type rk2_type = { "rk2", /* name */ + 1, /* can use dydt_in */ + 1, /* gives exact dydt_out */ + &rk2_alloc, + &rk2_apply, + &rk2_reset, + &rk2_order, + &rk2_free +}; + +const gsl_odeiv_step_type *gsl_odeiv_step_rk2 = &rk2_type; |