3 * Copyright (C) 1996, 1997, 1998, 1999, 2000 Gerard Jungman
5 * This program is free software; you can redistribute it and/or modify
6 * it under the terms of the GNU General Public License as published by
7 * the Free Software Foundation; either version 3 of the License, or (at
8 * your option) any later version.
10 * This program is distributed in the hope that it will be useful, but
11 * WITHOUT ANY WARRANTY; without even the implied warranty of
12 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
13 * General Public License for more details.
15 * You should have received a copy of the GNU General Public License
16 * along with this program; if not, write to the Free Software
17 * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
20 /* Gear 1. This is the implicit Euler a.k.a backward Euler method. */
25 /* Error estimation by step doubling, see eg. Ascher, U.M., Petzold,
26 L.R., Computer methods for ordinary differential and
27 differential-algebraic equations, SIAM, Philadelphia, 1998.
28 The method is also described in eg. this reference.
34 #include <gsl/gsl_math.h>
35 #include <gsl/gsl_errno.h>
36 #include <gsl/gsl_odeiv.h>
38 #include "odeiv_util.h"
50 gear1_alloc (size_t dim)
52 gear1_state_t *state = (gear1_state_t *) malloc (sizeof (gear1_state_t));
56 GSL_ERROR_NULL ("failed to allocate space for gear1_state", GSL_ENOMEM);
59 state->k = (double *) malloc (dim * sizeof (double));
64 GSL_ERROR_NULL ("failed to allocate space for k", GSL_ENOMEM);
67 state->y0 = (double *) malloc (dim * sizeof (double));
73 GSL_ERROR_NULL ("failed to allocate space for y0", GSL_ENOMEM);
76 state->y0_orig = (double *) malloc (dim * sizeof (double));
78 if (state->y0_orig == 0)
83 GSL_ERROR_NULL ("failed to allocate space for y0_orig", GSL_ENOMEM);
86 state->y_onestep = (double *) malloc (dim * sizeof (double));
88 if (state->y_onestep == 0)
90 free (state->y0_orig);
94 GSL_ERROR_NULL ("failed to allocate space for y_onestep", GSL_ENOMEM);
101 gear1_step (double *y, gear1_state_t *state,
102 const double h, const double t,
103 const size_t dim, const gsl_odeiv_system *sys)
105 /* Makes an implicit Euler advance with step size h.
106 y0 is the initial values of variables y.
108 The implicit matrix equations to solve are:
110 k = y0 + h * f(t + h, k)
112 y = y0 + h * f(t + h, k)
115 const int iter_steps = 3;
118 double *y0 = state->y0;
119 double *k = state->k;
121 /* Iterative solution of k = y0 + h * f(t + h, k)
123 Note: This method does not check for convergence of the
127 for (nu = 0; nu < iter_steps; nu++)
129 int s = GSL_ODEIV_FN_EVAL(sys, t + h, y, k);
131 if (s != GSL_SUCCESS)
136 for (i=0; i<dim; i++)
138 y[i] = y0[i] + h * k[i];
146 gear1_apply(void * vstate,
152 const double dydt_in[],
154 const gsl_odeiv_system * sys)
156 gear1_state_t *state = (gear1_state_t *) vstate;
160 double *y0 = state->y0;
161 double *y0_orig = state->y0_orig;
162 double *y_onestep = state->y_onestep;
165 DBL_MEMCPY(y0, y, dim);
167 /* Save initial values for possible failures */
168 DBL_MEMCPY (y0_orig, y, dim);
170 /* First traverse h with one step (save to y_onestep) */
171 DBL_MEMCPY (y_onestep, y, dim);
174 int s = gear1_step (y_onestep, state, h, t, dim, sys);
176 if (s != GSL_SUCCESS)
182 /* Then with two steps with half step length (save to y) */
184 int s = gear1_step (y, state, h / 2.0, t, dim, sys);
186 if (s != GSL_SUCCESS)
188 /* Restore original y vector */
189 DBL_MEMCPY (y, y0_orig, dim);
194 DBL_MEMCPY (y0, y, dim);
197 int s = gear1_step (y, state, h / 2.0, t + h / 2.0, dim, sys);
199 if (s != GSL_SUCCESS)
201 /* Restore original y vector */
202 DBL_MEMCPY (y, y0_orig, dim);
209 if (dydt_out != NULL)
211 int s = GSL_ODEIV_FN_EVAL (sys, t + h, y, dydt_out);
213 if (s != GSL_SUCCESS)
215 /* Restore original y vector */
216 DBL_MEMCPY (y, y0_orig, dim);
221 /* Error estimation */
223 for (i = 0; i < dim; i++)
225 yerr[i] = 4.0 * (y[i] - y_onestep[i]);
232 gear1_reset (void *vstate, size_t dim)
234 gear1_state_t *state = (gear1_state_t *) vstate;
236 DBL_ZERO_MEMSET (state->y_onestep, dim);
237 DBL_ZERO_MEMSET (state->y0_orig, dim);
238 DBL_ZERO_MEMSET (state->y0, dim);
239 DBL_ZERO_MEMSET (state->k, dim);
244 gear1_order (void *vstate)
246 gear1_state_t *state = (gear1_state_t *) vstate;
247 state = 0; /* prevent warnings about unused parameters */
252 gear1_free (void *vstate)
254 gear1_state_t *state = (gear1_state_t *) vstate;
255 free (state->y_onestep);
256 free (state->y0_orig);
262 static const gsl_odeiv_step_type gear1_type = { "gear1", /* name */
263 0, /* can use dydt_in? */
264 1, /* gives exact dydt_out? */
272 const gsl_odeiv_step_type *gsl_odeiv_step_gear1 = &gear1_type;