include/gsl/gsl_odeiv2.h

0001 /* ode-initval/odeiv2.h
0002  *
0003  * Copyright (C) 1996, 1997, 1998, 1999, 2000 Gerard Jungman
0004  *
0005  * This program is free software; you can redistribute it and/or modify
0006  * it under the terms of the GNU General Public License as published by
0007  * the Free Software Foundation; either version 3 of the License, or (at
0008  * your option) any later version.
0009  *
0010  * This program is distributed in the hope that it will be useful, but
0011  * WITHOUT ANY WARRANTY; without even the implied warranty of
0012  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
0013  * General Public License for more details.
0014  *
0015  * You should have received a copy of the GNU General Public License
0016  * along with this program; if not, write to the Free Software
0017  * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
0018  */
0019
0020 /* Author:  G. Jungman */
0021 /* Modified by Tuomo Keskitalo */
0022
0023 #ifndef __GSL_ODEIV2_H__
0024 #define __GSL_ODEIV2_H__
0025
0026 #include <stdio.h>
0027 #include <stdlib.h>
0028 #include <gsl/gsl_types.h>
0029
0030 #undef __BEGIN_DECLS
0031 #undef __END_DECLS
0032 #ifdef __cplusplus
0033 # define __BEGIN_DECLS extern "C" {
0034 # define __END_DECLS }
0035 #else
0036 # define __BEGIN_DECLS          /* empty */
0037 # define __END_DECLS            /* empty */
0038 #endif
0039
0040 __BEGIN_DECLS
0041 /* Description of a system of ODEs.
0042  *
0043  * y' = f(t,y) = dydt(t, y)
0044  *
0045  * The system is specified by giving the right-hand-side
0046  * of the equation and possibly a jacobian function.
0047  *
0048  * Some methods require the jacobian function, which calculates
0049  * the matrix dfdy and the vector dfdt. The matrix dfdy conforms
0050  * to the GSL standard, being a continuous range of floating point
0051  * values, in row-order.
0052  *
0053  * As with GSL function objects, user-supplied parameter
0054  * data is also present.
0055  */
0056   typedef struct
0057 {
0058   int (*function) (double t, const double y[], double dydt[], void *params);
0059   int (*jacobian) (double t, const double y[], double *dfdy, double dfdt[],
0060                    void *params);
0061   size_t dimension;
0062   void *params;
0063 }
0064 gsl_odeiv2_system;
0065
0066 /* Function evaluation macros */
0067
0068 #define GSL_ODEIV_FN_EVAL(S,t,y,f)  (*((S)->function))(t,y,f,(S)->params)
0069 #define GSL_ODEIV_JA_EVAL(S,t,y,dfdy,dfdt)  (*((S)->jacobian))(t,y,dfdy,dfdt,(S)->params)
0070
0071 /* Type definitions */
0072
0073 typedef struct gsl_odeiv2_step_struct gsl_odeiv2_step;
0074 typedef struct gsl_odeiv2_control_struct gsl_odeiv2_control;
0075 typedef struct gsl_odeiv2_evolve_struct gsl_odeiv2_evolve;
0076 typedef struct gsl_odeiv2_driver_struct gsl_odeiv2_driver;
0077
0078 /* Stepper object
0079  *
0080  * Opaque object for stepping an ODE system from t to t+h.
0081  * In general the object has some state which facilitates
0082  * iterating the stepping operation.
0083  */
0084
0085 typedef struct
0086 {
0087   const char *name;
0088   int can_use_dydt_in;
0089   int gives_exact_dydt_out;
0090   void *(*alloc) (size_t dim);
0091   int (*apply) (void *state, size_t dim, double t, double h, double y[],
0092                 double yerr[], const double dydt_in[], double dydt_out[],
0093                 const gsl_odeiv2_system * dydt);
0094   int (*set_driver) (void *state, const gsl_odeiv2_driver * d);
0095   int (*reset) (void *state, size_t dim);
0096   unsigned int (*order) (void *state);
0097   void (*free) (void *state);
0098 }
0099 gsl_odeiv2_step_type;
0100
0101 struct gsl_odeiv2_step_struct
0102 {
0103   const gsl_odeiv2_step_type *type;
0104   size_t dimension;
0105   void *state;
0106 };
0107
0108 /* Available stepper types */
0109
0110 GSL_VAR const gsl_odeiv2_step_type *gsl_odeiv2_step_rk2;
0111 GSL_VAR const gsl_odeiv2_step_type *gsl_odeiv2_step_rk4;
0112 GSL_VAR const gsl_odeiv2_step_type *gsl_odeiv2_step_rkf45;
0113 GSL_VAR const gsl_odeiv2_step_type *gsl_odeiv2_step_rkck;
0114 GSL_VAR const gsl_odeiv2_step_type *gsl_odeiv2_step_rk8pd;
0115 GSL_VAR const gsl_odeiv2_step_type *gsl_odeiv2_step_rk2imp;
0116 GSL_VAR const gsl_odeiv2_step_type *gsl_odeiv2_step_rk4imp;
0117 GSL_VAR const gsl_odeiv2_step_type *gsl_odeiv2_step_bsimp;
0118 GSL_VAR const gsl_odeiv2_step_type *gsl_odeiv2_step_rk1imp;
0119 GSL_VAR const gsl_odeiv2_step_type *gsl_odeiv2_step_msadams;
0120 GSL_VAR const gsl_odeiv2_step_type *gsl_odeiv2_step_msbdf;
0121
0122 /* Stepper object methods */
0123
0124 gsl_odeiv2_step *gsl_odeiv2_step_alloc (const gsl_odeiv2_step_type * T,
0125                                         size_t dim);
0126 int gsl_odeiv2_step_reset (gsl_odeiv2_step * s);
0127 void gsl_odeiv2_step_free (gsl_odeiv2_step * s);
0128 const char *gsl_odeiv2_step_name (const gsl_odeiv2_step * s);
0129 unsigned int gsl_odeiv2_step_order (const gsl_odeiv2_step * s);
0130 int gsl_odeiv2_step_apply (gsl_odeiv2_step * s, double t, double h,
0131                            double y[], double yerr[], const double dydt_in[],
0132                            double dydt_out[], const gsl_odeiv2_system * dydt);
0133 int gsl_odeiv2_step_set_driver (gsl_odeiv2_step * s,
0134                                 const gsl_odeiv2_driver * d);
0135
0136 /* Step size control object. */
0137
0138 typedef struct
0139 {
0140   const char *name;
0141   void *(*alloc) (void);
0142   int (*init) (void *state, double eps_abs, double eps_rel, double a_y,
0143                double a_dydt);
0144   int (*hadjust) (void *state, size_t dim, unsigned int ord, const double y[],
0145                   const double yerr[], const double yp[], double *h);
0146   int (*errlevel) (void *state, const double y, const double dydt,
0147                    const double h, const size_t ind, double *errlev);
0148   int (*set_driver) (void *state, const gsl_odeiv2_driver * d);
0149   void (*free) (void *state);
0150 }
0151 gsl_odeiv2_control_type;
0152
0153 struct gsl_odeiv2_control_struct
0154 {
0155   const gsl_odeiv2_control_type *type;
0156   void *state;
0157 };
0158
0159 /* Possible return values for an hadjust() evolution method */
0160
0161 #define GSL_ODEIV_HADJ_INC   1  /* step was increased */
0162 #define GSL_ODEIV_HADJ_NIL   0  /* step unchanged     */
0163 #define GSL_ODEIV_HADJ_DEC (-1) /* step decreased     */
0164
0165 /* General step size control methods.
0166  *
0167  * The hadjust() method controls the adjustment of
0168  * step size given the result of a step and the error.
0169  * Valid hadjust() methods must return one of the codes below.
0170  * errlevel function calculates the desired error level D0.
0171  *
0172  * The general data can be used by specializations
0173  * to store state and control their heuristics.
0174  */
0175
0176 gsl_odeiv2_control *gsl_odeiv2_control_alloc (const gsl_odeiv2_control_type *
0177                                               T);
0178 int gsl_odeiv2_control_init (gsl_odeiv2_control * c, double eps_abs,
0179                              double eps_rel, double a_y, double a_dydt);
0180 void gsl_odeiv2_control_free (gsl_odeiv2_control * c);
0181 int gsl_odeiv2_control_hadjust (gsl_odeiv2_control * c, gsl_odeiv2_step * s,
0182                                 const double y[], const double yerr[],
0183                                 const double dydt[], double *h);
0184 const char *gsl_odeiv2_control_name (const gsl_odeiv2_control * c);
0185 int gsl_odeiv2_control_errlevel (gsl_odeiv2_control * c, const double y,
0186                                  const double dydt, const double h,
0187                                  const size_t ind, double *errlev);
0188 int gsl_odeiv2_control_set_driver (gsl_odeiv2_control * c,
0189                                    const gsl_odeiv2_driver * d);
0190
0191 /* Available control object constructors.
0192  *
0193  * The standard control object is a four parameter heuristic
0194  * defined as follows:
0195  *    D0 = eps_abs + eps_rel * (a_y |y| + a_dydt h |y'|)
0196  *    D1 = |yerr|
0197  *    q  = consistency order of method (q=4 for 4(5) embedded RK)
0198  *    S  = safety factor (0.9 say)
0199  *
0200  *                      /  (D0/D1)^(1/(q+1))  D0 >= D1
0201  *    h_NEW = S h_OLD * |
0202  *                      \  (D0/D1)^(1/q)      D0 < D1
0203  *
0204  * This encompasses all the standard error scaling methods.
0205  *
0206  * The y method is the standard method with a_y=1, a_dydt=0.
0207  * The yp method is the standard method with a_y=0, a_dydt=1.
0208  */
0209
0210 gsl_odeiv2_control *gsl_odeiv2_control_standard_new (double eps_abs,
0211                                                      double eps_rel,
0212                                                      double a_y,
0213                                                      double a_dydt);
0214 gsl_odeiv2_control *gsl_odeiv2_control_y_new (double eps_abs, double eps_rel);
0215 gsl_odeiv2_control *gsl_odeiv2_control_yp_new (double eps_abs,
0216                                                double eps_rel);
0217
0218 /* This controller computes errors using different absolute errors for
0219  * each component
0220  *
0221  *    D0 = eps_abs * scale_abs[i] + eps_rel * (a_y |y| + a_dydt h |y'|)
0222  */
0223
0224 gsl_odeiv2_control *gsl_odeiv2_control_scaled_new (double eps_abs,
0225                                                    double eps_rel, double a_y,
0226                                                    double a_dydt,
0227                                                    const double scale_abs[],
0228                                                    size_t dim);
0229
0230 /* Evolution object */
0231
0232 struct gsl_odeiv2_evolve_struct
0233 {
0234   size_t dimension;
0235   double *y0;
0236   double *yerr;
0237   double *dydt_in;
0238   double *dydt_out;
0239   double last_step;
0240   unsigned long int count;
0241   unsigned long int failed_steps;
0242   const gsl_odeiv2_driver *driver;
0243 };
0244
0245 /* Evolution object methods */
0246
0247 gsl_odeiv2_evolve *gsl_odeiv2_evolve_alloc (size_t dim);
0248 int gsl_odeiv2_evolve_apply (gsl_odeiv2_evolve * e, gsl_odeiv2_control * con,
0249                              gsl_odeiv2_step * step,
0250                              const gsl_odeiv2_system * dydt, double *t,
0251                              double t1, double *h, double y[]);
0252 int gsl_odeiv2_evolve_apply_fixed_step (gsl_odeiv2_evolve * e,
0253                                         gsl_odeiv2_control * con,
0254                                         gsl_odeiv2_step * step,
0255                                         const gsl_odeiv2_system * dydt,
0256                                         double *t, const double h0,
0257                                         double y[]);
0258 int gsl_odeiv2_evolve_reset (gsl_odeiv2_evolve * e);
0259 void gsl_odeiv2_evolve_free (gsl_odeiv2_evolve * e);
0260 int gsl_odeiv2_evolve_set_driver (gsl_odeiv2_evolve * e,
0261                                   const gsl_odeiv2_driver * d);
0262
0263 /* Driver object
0264  *
0265  * This is a high level wrapper for step, control and
0266  * evolve objects.
0267  */
0268
0269 struct gsl_odeiv2_driver_struct
0270 {
0271   const gsl_odeiv2_system *sys; /* ODE system */
0272   gsl_odeiv2_step *s;           /* stepper object */
0273   gsl_odeiv2_control *c;        /* control object */
0274   gsl_odeiv2_evolve *e;         /* evolve object */
0275   double h;                     /* step size */
0276   double hmin;                  /* minimum step size allowed */
0277   double hmax;                  /* maximum step size allowed */
0278   unsigned long int n;          /* number of steps taken */
0279   unsigned long int nmax;       /* Maximum number of steps allowed */
0280 };
0281
0282 /* Driver object methods */
0283
0284 gsl_odeiv2_driver *gsl_odeiv2_driver_alloc_y_new (const gsl_odeiv2_system *
0285                                                   sys,
0286                                                   const gsl_odeiv2_step_type *
0287                                                   T, const double hstart,
0288                                                   const double epsabs,
0289                                                   const double epsrel);
0290 gsl_odeiv2_driver *gsl_odeiv2_driver_alloc_yp_new (const gsl_odeiv2_system *
0291                                                    sys,
0292                                                    const gsl_odeiv2_step_type
0293                                                    * T, const double hstart,
0294                                                    const double epsabs,
0295                                                    const double epsrel);
0296 gsl_odeiv2_driver *gsl_odeiv2_driver_alloc_scaled_new (const gsl_odeiv2_system
0297                                                        * sys,
0298                                                        const
0299                                                        gsl_odeiv2_step_type *
0300                                                        T, const double hstart,
0301                                                        const double epsabs,
0302                                                        const double epsrel,
0303                                                        const double a_y,
0304                                                        const double a_dydt,
0305                                                        const double
0306                                                        scale_abs[]);
0307 gsl_odeiv2_driver *gsl_odeiv2_driver_alloc_standard_new (const
0308                                                          gsl_odeiv2_system *
0309                                                          sys,
0310                                                          const
0311                                                          gsl_odeiv2_step_type
0312                                                          * T,
0313                                                          const double hstart,
0314                                                          const double epsabs,
0315                                                          const double epsrel,
0316                                                          const double a_y,
0317                                                          const double a_dydt);
0318 int gsl_odeiv2_driver_set_hmin (gsl_odeiv2_driver * d, const double hmin);
0319 int gsl_odeiv2_driver_set_hmax (gsl_odeiv2_driver * d, const double hmax);
0320 int gsl_odeiv2_driver_set_nmax (gsl_odeiv2_driver * d,
0321                                 const unsigned long int nmax);
0322 int gsl_odeiv2_driver_apply (gsl_odeiv2_driver * d, double *t,
0323                              const double t1, double y[]);
0324 int gsl_odeiv2_driver_apply_fixed_step (gsl_odeiv2_driver * d, double *t,
0325                                         const double h,
0326                                         const unsigned long int n,
0327                                         double y[]);
0328 int gsl_odeiv2_driver_reset (gsl_odeiv2_driver * d);
0329 int gsl_odeiv2_driver_reset_hstart (gsl_odeiv2_driver * d, const double hstart);
0330 void gsl_odeiv2_driver_free (gsl_odeiv2_driver * state);
0331
0332 __END_DECLS
0333 #endif /* __GSL_ODEIV2_H__ */