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
58typedef 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}
65gsl_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
78typedef 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}
89gsl_odeiv_step_type;
90
91typedef struct {
92 const gsl_odeiv_step_type * type;
93 size_t dimension;
94 void * state;
95}
96gsl_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
111GSL_VAR const gsl_odeiv_step_type *gsl_odeiv_step_rk2;
112GSL_VAR const gsl_odeiv_step_type *gsl_odeiv_step_rk4;
113GSL_VAR const gsl_odeiv_step_type *gsl_odeiv_step_rkf45;
114GSL_VAR const gsl_odeiv_step_type *gsl_odeiv_step_rkck;
115GSL_VAR const gsl_odeiv_step_type *gsl_odeiv_step_rk8pd;
116GSL_VAR const gsl_odeiv_step_type *gsl_odeiv_step_rk2imp;
117GSL_VAR const gsl_odeiv_step_type *gsl_odeiv_step_rk2simp;
118GSL_VAR const gsl_odeiv_step_type *gsl_odeiv_step_rk4imp;
119GSL_VAR const gsl_odeiv_step_type *gsl_odeiv_step_bsimp;
120GSL_VAR const gsl_odeiv_step_type *gsl_odeiv_step_gear1;
121GSL_VAR const gsl_odeiv_step_type *gsl_odeiv_step_gear2;
122
123
124/* Constructor for specialized stepper objects.
125 */
126gsl_odeiv_step * gsl_odeiv_step_alloc(const gsl_odeiv_step_type * T, size_t dim);
127int gsl_odeiv_step_reset(gsl_odeiv_step * s);
128void gsl_odeiv_step_free(gsl_odeiv_step * s);
129
130/* General stepper object methods.
131 */
132const char * gsl_odeiv_step_name(const gsl_odeiv_step * s);
133unsigned int gsl_odeiv_step_order(const gsl_odeiv_step * s);
134
135int gsl_odeiv_step_apply(gsl_odeiv_step * s, 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
147typedef 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}
155gsl_odeiv_control_type;
156
157typedef struct
158{
159 const gsl_odeiv_control_type * type;
160 void * state;
161}
162gsl_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
170gsl_odeiv_control * gsl_odeiv_control_alloc(const gsl_odeiv_control_type * T);
171int gsl_odeiv_control_init(gsl_odeiv_control * c, double eps_abs, double eps_rel, double a_y, double a_dydt);
172void gsl_odeiv_control_free(gsl_odeiv_control * c);
173int gsl_odeiv_control_hadjust (gsl_odeiv_control * c, gsl_odeiv_step * s, const double y[], const double yerr[], const double dydt[], double * h);
174const 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
195gsl_odeiv_control * gsl_odeiv_control_standard_new(double eps_abs, double eps_rel, double a_y, double a_dydt);
196gsl_odeiv_control * gsl_odeiv_control_y_new(double eps_abs, double eps_rel);
197gsl_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 */
204gsl_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 */
208typedef 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}
218gsl_odeiv_evolve;
219
220/* Evolution object methods.
221 */
222gsl_odeiv_evolve * gsl_odeiv_evolve_alloc(size_t dim);
223int gsl_odeiv_evolve_apply(gsl_odeiv_evolve * e, gsl_odeiv_control * con, gsl_odeiv_step * step, const gsl_odeiv_system * dydt, double * t, double t1, double * h, double y[]);
224int gsl_odeiv_evolve_reset(gsl_odeiv_evolve * e);
225void gsl_odeiv_evolve_free(gsl_odeiv_evolve * e);
226
227
228__END_DECLS
229
230#endif /* __GSL_ODEIV_H__ */
231