1 /* ode-initval/odeiv2.h 2 * 3 * Copyright (C) 1996, 1997, 1998, 1999, 2000 Gerard Jungman 4 * 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. 9 * 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. 14 * 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. 18 */ 19 20 /* Author: G. Jungman */ 21 /* Modified by Tuomo Keskitalo */ 22 23 #ifndef __GSL_ODEIV2_H__ 24 #define __GSL_ODEIV2_H__ 25 26 #include <stdio.h> 27 #include <stdlib.h> 28 #include <gsl/gsl_types.h> 29 30 #undef __BEGIN_DECLS 31 #undef __END_DECLS 32 #ifdef __cplusplus 33 # define __BEGIN_DECLS extern "C" { 34 # define __END_DECLS } 35 #else 36 # define __BEGIN_DECLS /* empty */ 37 # define __END_DECLS /* empty */ 38 #endif 39 40 __BEGIN_DECLS 41 /* Description of a system of ODEs. 42 * 43 * y' = f(t,y) = dydt(t, y) 44 * 45 * The system is specified by giving the right-hand-side 46 * of the equation and possibly a jacobian function. 47 * 48 * Some methods require the jacobian function, which calculates 49 * the matrix dfdy and the vector dfdt. The matrix dfdy conforms 50 * to the GSL standard, being a continuous range of floating point 51 * values, in row-order. 52 * 53 * As with GSL function objects, user-supplied parameter 54 * data is also present. 55 */ 56 typedef struct 57 { 58 int (*function) (double t, const double y[], double dydt[], void *params); 59 int (*jacobian) (double t, const double y[], double *dfdy, double dfdt[], 60 void *params); 61 size_t dimension; 62 void *params; 63 } 64 gsl_odeiv2_system; 65 66 /* Function evaluation macros */ 67 68 #define GSL_ODEIV_FN_EVAL(S,t,y,f) (*((S)->function))(t,y,f,(S)->params) 69 #define GSL_ODEIV_JA_EVAL(S,t,y,dfdy,dfdt) (*((S)->jacobian))(t,y,dfdy,dfdt,(S)->params) 70 71 /* Type definitions */ 72 73 typedef struct gsl_odeiv2_step_struct gsl_odeiv2_step; 74 typedef struct gsl_odeiv2_control_struct gsl_odeiv2_control; 75 typedef struct gsl_odeiv2_evolve_struct gsl_odeiv2_evolve; 76 typedef struct gsl_odeiv2_driver_struct gsl_odeiv2_driver; 77 78 /* Stepper object 79 * 80 * Opaque object for stepping an ODE system from t to t+h. 81 * In general the object has some state which facilitates 82 * iterating the stepping operation. 83 */ 84 85 typedef struct 86 { 87 const char *name; 88 int can_use_dydt_in; 89 int gives_exact_dydt_out; 90 void *(*alloc) (size_t dim); 91 int (*apply) (void *state, size_t dim, double t, double h, double y[], 92 double yerr[], const double dydt_in[], double dydt_out[], 93 const gsl_odeiv2_system * dydt); 94 int (*set_driver) (void *state, const gsl_odeiv2_driver * d); 95 int (*reset) (void *state, size_t dim); 96 unsigned int (*order) (void *state); 97 void (*free) (void *state); 98 } 99 gsl_odeiv2_step_type; 100 101 struct gsl_odeiv2_step_struct 102 { 103 const gsl_odeiv2_step_type *type; 104 size_t dimension; 105 void *state; 106 }; 107 108 /* Available stepper types */ 109 110 GSL_VAR const gsl_odeiv2_step_type *gsl_odeiv2_step_rk2; 111 GSL_VAR const gsl_odeiv2_step_type *gsl_odeiv2_step_rk4; 112 GSL_VAR const gsl_odeiv2_step_type *gsl_odeiv2_step_rkf45; 113 GSL_VAR const gsl_odeiv2_step_type *gsl_odeiv2_step_rkck; 114 GSL_VAR const gsl_odeiv2_step_type *gsl_odeiv2_step_rk8pd; 115 GSL_VAR const gsl_odeiv2_step_type *gsl_odeiv2_step_rk2imp; 116 GSL_VAR const gsl_odeiv2_step_type *gsl_odeiv2_step_rk4imp; 117 GSL_VAR const gsl_odeiv2_step_type *gsl_odeiv2_step_bsimp; 118 GSL_VAR const gsl_odeiv2_step_type *gsl_odeiv2_step_rk1imp; 119 GSL_VAR const gsl_odeiv2_step_type *gsl_odeiv2_step_msadams; 120 GSL_VAR const gsl_odeiv2_step_type *gsl_odeiv2_step_msbdf; 121 122 /* Stepper object methods */ 123 124 gsl_odeiv2_step *gsl_odeiv2_step_alloc (const gsl_odeiv2_step_type * T, 125 size_t dim); 126 int gsl_odeiv2_step_reset (gsl_odeiv2_step * s); 127 void gsl_odeiv2_step_free (gsl_odeiv2_step * s); 128 const char *gsl_odeiv2_step_name (const gsl_odeiv2_step * s); 129 unsigned int gsl_odeiv2_step_order (const gsl_odeiv2_step * s); 130 int gsl_odeiv2_step_apply (gsl_odeiv2_step * s, double t, double h, 131 double y[], double yerr[], const double dydt_in[], 132 double dydt_out[], const gsl_odeiv2_system * dydt); 133 int gsl_odeiv2_step_set_driver (gsl_odeiv2_step * s, 134 const gsl_odeiv2_driver * d); 135 136 /* Step size control object. */ 137 138 typedef struct 139 { 140 const char *name; 141 void *(*alloc) (void); 142 int (*init) (void *state, double eps_abs, double eps_rel, double a_y, 143 double a_dydt); 144 int (*hadjust) (void *state, size_t dim, unsigned int ord, const double y[], 145 const double yerr[], const double yp[], double *h); 146 int (*errlevel) (void *state, const double y, const double dydt, 147 const double h, const size_t ind, double *errlev); 148 int (*set_driver) (void *state, const gsl_odeiv2_driver * d); 149 void (*free) (void *state); 150 } 151 gsl_odeiv2_control_type; 152 153 struct gsl_odeiv2_control_struct 154 { 155 const gsl_odeiv2_control_type *type; 156 void *state; 157 }; 158 159 /* Possible return values for an hadjust() evolution method */ 160 161 #define GSL_ODEIV_HADJ_INC 1 /* step was increased */ 162 #define GSL_ODEIV_HADJ_NIL 0 /* step unchanged */ 163 #define GSL_ODEIV_HADJ_DEC (-1) /* step decreased */ 164 165 /* General step size control methods. 166 * 167 * The hadjust() method controls the adjustment of 168 * step size given the result of a step and the error. 169 * Valid hadjust() methods must return one of the codes below. 170 * errlevel function calculates the desired error level D0. 171 * 172 * The general data can be used by specializations 173 * to store state and control their heuristics. 174 */ 175 176 gsl_odeiv2_control *gsl_odeiv2_control_alloc (const gsl_odeiv2_control_type * 177 T); 178 int gsl_odeiv2_control_init (gsl_odeiv2_control * c, double eps_abs, 179 double eps_rel, double a_y, double a_dydt); 180 void gsl_odeiv2_control_free (gsl_odeiv2_control * c); 181 int gsl_odeiv2_control_hadjust (gsl_odeiv2_control * c, gsl_odeiv2_step * s, 182 const double y[], const double yerr[], 183 const double dydt[], double *h); 184 const char *gsl_odeiv2_control_name (const gsl_odeiv2_control * c); 185 int gsl_odeiv2_control_errlevel (gsl_odeiv2_control * c, const double y, 186 const double dydt, const double h, 187 const size_t ind, double *errlev); 188 int gsl_odeiv2_control_set_driver (gsl_odeiv2_control * c, 189 const gsl_odeiv2_driver * d); 190 191 /* Available control object constructors. 192 * 193 * The standard control object is a four parameter heuristic 194 * defined as follows: 195 * D0 = eps_abs + eps_rel * (a_y |y| + a_dydt h |y'|) 196 * D1 = |yerr| 197 * q = consistency order of method (q=4 for 4(5) embedded RK) 198 * S = safety factor (0.9 say) 199 * 200 * / (D0/D1)^(1/(q+1)) D0 >= D1 201 * h_NEW = S h_OLD * | 202 * \ (D0/D1)^(1/q) D0 < D1 203 * 204 * This encompasses all the standard error scaling methods. 205 * 206 * The y method is the standard method with a_y=1, a_dydt=0. 207 * The yp method is the standard method with a_y=0, a_dydt=1. 208 */ 209 210 gsl_odeiv2_control *gsl_odeiv2_control_standard_new (double eps_abs, 211 double eps_rel, 212 double a_y, 213 double a_dydt); 214 gsl_odeiv2_control *gsl_odeiv2_control_y_new (double eps_abs, double eps_rel); 215 gsl_odeiv2_control *gsl_odeiv2_control_yp_new (double eps_abs, 216 double eps_rel); 217 218 /* This controller computes errors using different absolute errors for 219 * each component 220 * 221 * D0 = eps_abs * scale_abs[i] + eps_rel * (a_y |y| + a_dydt h |y'|) 222 */ 223 224 gsl_odeiv2_control *gsl_odeiv2_control_scaled_new (double eps_abs, 225 double eps_rel, double a_y, 226 double a_dydt, 227 const double scale_abs[], 228 size_t dim); 229 230 /* Evolution object */ 231 232 struct gsl_odeiv2_evolve_struct 233 { 234 size_t dimension; 235 double *y0; 236 double *yerr; 237 double *dydt_in; 238 double *dydt_out; 239 double last_step; 240 unsigned long int count; 241 unsigned long int failed_steps; 242 const gsl_odeiv2_driver *driver; 243 }; 244 245 /* Evolution object methods */ 246 247 gsl_odeiv2_evolve *gsl_odeiv2_evolve_alloc (size_t dim); 248 int gsl_odeiv2_evolve_apply (gsl_odeiv2_evolve * e, gsl_odeiv2_control * con, 249 gsl_odeiv2_step * step, 250 const gsl_odeiv2_system * dydt, double *t, 251 double t1, double *h, double y[]); 252 int gsl_odeiv2_evolve_apply_fixed_step (gsl_odeiv2_evolve * e, 253 gsl_odeiv2_control * con, 254 gsl_odeiv2_step * step, 255 const gsl_odeiv2_system * dydt, 256 double *t, const double h0, 257 double y[]); 258 int gsl_odeiv2_evolve_reset (gsl_odeiv2_evolve * e); 259 void gsl_odeiv2_evolve_free (gsl_odeiv2_evolve * e); 260 int gsl_odeiv2_evolve_set_driver (gsl_odeiv2_evolve * e, 261 const gsl_odeiv2_driver * d); 262 263 /* Driver object 264 * 265 * This is a high level wrapper for step, control and 266 * evolve objects. 267 */ 268 269 struct gsl_odeiv2_driver_struct 270 { 271 const gsl_odeiv2_system *sys; /* ODE system */ 272 gsl_odeiv2_step *s; /* stepper object */ 273 gsl_odeiv2_control *c; /* control object */ 274 gsl_odeiv2_evolve *e; /* evolve object */ 275 double h; /* step size */ 276 double hmin; /* minimum step size allowed */ 277 double hmax; /* maximum step size allowed */ 278 unsigned long int n; /* number of steps taken */ 279 unsigned long int nmax; /* Maximum number of steps allowed */ 280 }; 281 282 /* Driver object methods */ 283 284 gsl_odeiv2_driver *gsl_odeiv2_driver_alloc_y_new (const gsl_odeiv2_system * 285 sys, 286 const gsl_odeiv2_step_type * 287 T, const double hstart, 288 const double epsabs, 289 const double epsrel); 290 gsl_odeiv2_driver *gsl_odeiv2_driver_alloc_yp_new (const gsl_odeiv2_system * 291 sys, 292 const gsl_odeiv2_step_type 293 * T, const double hstart, 294 const double epsabs, 295 const double epsrel); 296 gsl_odeiv2_driver *gsl_odeiv2_driver_alloc_scaled_new (const gsl_odeiv2_system 297 * sys, 298 const 299 gsl_odeiv2_step_type * 300 T, const double hstart, 301 const double epsabs, 302 const double epsrel, 303 const double a_y, 304 const double a_dydt, 305 const double 306 scale_abs[]); 307 gsl_odeiv2_driver *gsl_odeiv2_driver_alloc_standard_new (const 308 gsl_odeiv2_system * 309 sys, 310 const 311 gsl_odeiv2_step_type 312 * T, 313 const double hstart, 314 const double epsabs, 315 const double epsrel, 316 const double a_y, 317 const double a_dydt); 318 int gsl_odeiv2_driver_set_hmin (gsl_odeiv2_driver * d, const double hmin); 319 int gsl_odeiv2_driver_set_hmax (gsl_odeiv2_driver * d, const double hmax); 320 int gsl_odeiv2_driver_set_nmax (gsl_odeiv2_driver * d, 321 const unsigned long int nmax); 322 int gsl_odeiv2_driver_apply (gsl_odeiv2_driver * d, double *t, 323 const double t1, double y[]); 324 int gsl_odeiv2_driver_apply_fixed_step (gsl_odeiv2_driver * d, double *t, 325 const double h, 326 const unsigned long int n, 327 double y[]); 328 int gsl_odeiv2_driver_reset (gsl_odeiv2_driver * d); 329 int gsl_odeiv2_driver_reset_hstart (gsl_odeiv2_driver * d, const double hstart); 330 void gsl_odeiv2_driver_free (gsl_odeiv2_driver * state); 331 332 __END_DECLS 333 #endif /* __GSL_ODEIV2_H__ */ 334