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 * s); |
133 | unsigned int gsl_odeiv_step_order(const gsl_odeiv_step * s); |
134 | |
135 | int 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 | |
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 * e, 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 * e); |
225 | void gsl_odeiv_evolve_free(gsl_odeiv_evolve * e); |
226 | |
227 | |
228 | __END_DECLS |
229 | |
230 | #endif /* __GSL_ODEIV_H__ */ |
231 | |