1 /* ode-initval/gsl_odeiv.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 */ 22 #ifndef __GSL_ODEIV_H__ 23 #define __GSL_ODEIV_H__ 24 25 #include <stdio.h> 26 #include <stdlib.h> 27 #include "gsl_types.h" 28 29 #undef __BEGIN_DECLS 30 #undef __END_DECLS 31 #ifdef __cplusplus 32 # define __BEGIN_DECLS extern "C" { 33 # define __END_DECLS } 34 #else 35 # define __BEGIN_DECLS /* empty */ 36 # define __END_DECLS /* empty */ 37 #endif 38 39 __BEGIN_DECLS 40 41 42 /* Description of a system of ODEs. 43 * 44 * y' = f(t,y) = dydt(t, y) 45 * 46 * The system is specified by giving the right-hand-side 47 * of the equation and possibly a jacobian function. 48 * 49 * Some methods require the jacobian function, which calculates 50 * the matrix dfdy and the vector dfdt. The matrix dfdy conforms 51 * to the GSL standard, being a continuous range of floating point 52 * values, in row-order. 53 * 54 * As with GSL function objects, user-supplied parameter 55 * data is also present. 56 */ 57 58 typedef struct 59 { 60 int (* function) (double t, const double y[], double dydt[], void * params); 61 int (* jacobian) (double t, const double y[], double * dfdy, double dfdt[], void * params); 62 size_t dimension; 63 void * params; 64 } 65 gsl_odeiv_system; 66 67 #define GSL_ODEIV_FN_EVAL(S,t,y,f) (*((S)->function))(t,y,f,(S)->params) 68 #define GSL_ODEIV_JA_EVAL(S,t,y,dfdy,dfdt) (*((S)->jacobian))(t,y,dfdy,dfdt,(S)->params) 69 70 71 /* General stepper object. 72 * 73 * Opaque object for stepping an ODE system from t to t+h. 74 * In general the object has some state which facilitates 75 * iterating the stepping operation. 76 */ 77 78 typedef struct 79 { 80 const char * name; 81 int can_use_dydt_in; 82 int gives_exact_dydt_out; 83 void * (*alloc) (size_t dim); 84 int (*apply) (void * state, size_t dim, double t, double h, double y[], double yerr[], const double dydt_in[], double dydt_out[], const gsl_odeiv_system * dydt); 85 int (*reset) (void * state, size_t dim); 86 unsigned int (*order) (void * state); 87 void (*free) (void * state); 88 } 89 gsl_odeiv_step_type; 90 91 typedef struct { 92 const gsl_odeiv_step_type * type; 93 size_t dimension; 94 void * state; 95 } 96 gsl_odeiv_step; 97 98 99 /* Available stepper types. 100 * 101 * rk2 : embedded 2nd(3rd) Runge-Kutta 102 * rk4 : 4th order (classical) Runge-Kutta 103 * rkck : embedded 4th(5th) Runge-Kutta, Cash-Karp 104 * rk8pd : embedded 8th(9th) Runge-Kutta, Prince-Dormand 105 * rk2imp : implicit 2nd order Runge-Kutta at Gaussian points 106 * rk4imp : implicit 4th order Runge-Kutta at Gaussian points 107 * gear1 : M=1 implicit Gear method 108 * gear2 : M=2 implicit Gear method 109 */ 110 111 GSL_VAR const gsl_odeiv_step_type *gsl_odeiv_step_rk2; 112 GSL_VAR const gsl_odeiv_step_type *gsl_odeiv_step_rk4; 113 GSL_VAR const gsl_odeiv_step_type *gsl_odeiv_step_rkf45; 114 GSL_VAR const gsl_odeiv_step_type *gsl_odeiv_step_rkck; 115 GSL_VAR const gsl_odeiv_step_type *gsl_odeiv_step_rk8pd; 116 GSL_VAR const gsl_odeiv_step_type *gsl_odeiv_step_rk2imp; 117 GSL_VAR const gsl_odeiv_step_type *gsl_odeiv_step_rk2simp; 118 GSL_VAR const gsl_odeiv_step_type *gsl_odeiv_step_rk4imp; 119 GSL_VAR const gsl_odeiv_step_type *gsl_odeiv_step_bsimp; 120 GSL_VAR const gsl_odeiv_step_type *gsl_odeiv_step_gear1; 121 GSL_VAR const gsl_odeiv_step_type *gsl_odeiv_step_gear2; 122 123 124 /* Constructor for specialized stepper objects. 125 */ 126 gsl_odeiv_step * gsl_odeiv_step_alloc(const gsl_odeiv_step_type * T, size_t dim); 127 int gsl_odeiv_step_reset(gsl_odeiv_step * s); 128 void gsl_odeiv_step_free(gsl_odeiv_step * s); 129 130 /* General stepper object methods. 131 */ 132 const char * gsl_odeiv_step_name(const gsl_odeiv_step *); 133 unsigned int gsl_odeiv_step_order(const gsl_odeiv_step * s); 134 135 int gsl_odeiv_step_apply(gsl_odeiv_step *, double t, double h, double y[], double yerr[], const double dydt_in[], double dydt_out[], const gsl_odeiv_system * dydt); 136 137 /* General step size control object. 138 * 139 * The hadjust() method controls the adjustment of 140 * step size given the result of a step and the error. 141 * Valid hadjust() methods must return one of the codes below. 142 * 143 * The general data can be used by specializations 144 * to store state and control their heuristics. 145 */ 146 147 typedef struct 148 { 149 const char * name; 150 void * (*alloc) (void); 151 int (*init) (void * state, double eps_abs, double eps_rel, double a_y, double a_dydt); 152 int (*hadjust) (void * state, size_t dim, unsigned int ord, const double y[], const double yerr[], const double yp[], double * h); 153 void (*free) (void * state); 154 } 155 gsl_odeiv_control_type; 156 157 typedef struct 158 { 159 const gsl_odeiv_control_type * type; 160 void * state; 161 } 162 gsl_odeiv_control; 163 164 /* Possible return values for an hadjust() evolution method. 165 */ 166 #define GSL_ODEIV_HADJ_INC 1 /* step was increased */ 167 #define GSL_ODEIV_HADJ_NIL 0 /* step unchanged */ 168 #define GSL_ODEIV_HADJ_DEC (-1) /* step decreased */ 169 170 gsl_odeiv_control * gsl_odeiv_control_alloc(const gsl_odeiv_control_type * T); 171 int gsl_odeiv_control_init(gsl_odeiv_control * c, double eps_abs, double eps_rel, double a_y, double a_dydt); 172 void gsl_odeiv_control_free(gsl_odeiv_control * c); 173 int gsl_odeiv_control_hadjust (gsl_odeiv_control * c, gsl_odeiv_step * s, const double y[], const double yerr[], const double dydt[], double * h); 174 const char * gsl_odeiv_control_name(const gsl_odeiv_control * c); 175 176 /* Available control object constructors. 177 * 178 * The standard control object is a four parameter heuristic 179 * defined as follows: 180 * D0 = eps_abs + eps_rel * (a_y |y| + a_dydt h |y'|) 181 * D1 = |yerr| 182 * q = consistency order of method (q=4 for 4(5) embedded RK) 183 * S = safety factor (0.9 say) 184 * 185 * / (D0/D1)^(1/(q+1)) D0 >= D1 186 * h_NEW = S h_OLD * | 187 * \ (D0/D1)^(1/q) D0 < D1 188 * 189 * This encompasses all the standard error scaling methods. 190 * 191 * The y method is the standard method with a_y=1, a_dydt=0. 192 * The yp method is the standard method with a_y=0, a_dydt=1. 193 */ 194 195 gsl_odeiv_control * gsl_odeiv_control_standard_new(double eps_abs, double eps_rel, double a_y, double a_dydt); 196 gsl_odeiv_control * gsl_odeiv_control_y_new(double eps_abs, double eps_rel); 197 gsl_odeiv_control * gsl_odeiv_control_yp_new(double eps_abs, double eps_rel); 198 199 /* This controller computes errors using different absolute errors for 200 * each component 201 * 202 * D0 = eps_abs * scale_abs[i] + eps_rel * (a_y |y| + a_dydt h |y'|) 203 */ 204 gsl_odeiv_control * gsl_odeiv_control_scaled_new(double eps_abs, double eps_rel, double a_y, double a_dydt, const double scale_abs[], size_t dim); 205 206 /* General evolution object. 207 */ 208 typedef struct { 209 size_t dimension; 210 double * y0; 211 double * yerr; 212 double * dydt_in; 213 double * dydt_out; 214 double last_step; 215 unsigned long int count; 216 unsigned long int failed_steps; 217 } 218 gsl_odeiv_evolve; 219 220 /* Evolution object methods. 221 */ 222 gsl_odeiv_evolve * gsl_odeiv_evolve_alloc(size_t dim); 223 int gsl_odeiv_evolve_apply(gsl_odeiv_evolve *, gsl_odeiv_control * con, gsl_odeiv_step * step, const gsl_odeiv_system * dydt, double * t, double t1, double * h, double y[]); 224 int gsl_odeiv_evolve_reset(gsl_odeiv_evolve *); 225 void gsl_odeiv_evolve_free(gsl_odeiv_evolve *); 226 227 228 __END_DECLS 229 230 #endif /* __GSL_ODEIV_H__ */ 231