htswanalysis/MACS/lib/gsl/gsl-1.11/ode-initval/gsl_odeiv.h

   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/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__ */