/[ascend]/trunk/solvers/dopri5/dopri5.h
ViewVC logotype

Contents of /trunk/solvers/dopri5/dopri5.h

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1733 - (show annotations) (download) (as text)
Wed Feb 6 00:18:31 2008 UTC (12 years ago) by jpye
File MIME type: text/x-chdr
File size: 9805 byte(s)
In process of moving DOPRI5 solver to 'solvers' dir.
1 /*
2 This is a modified form of the DOPRI5 code, originally written in Fortran 77
3 by E Hairer and G Wanner, then translated to C by J Colinge.
4 Changes have been made to suite ASCEND requirements. Minor only.
5 -- John Pye, May 2007.
6 *//*
7 DOPRI5
8 Copyright (c) 2004, Ernst Hairer
9
10 Redistribution and use in source and binary forms, with or without
11 modification, are permitted provided that the following conditions are
12 met:
13
14 - Redistributions of source code must retain the above copyright
15 notice, this list of conditions and the following disclaimer.
16
17 - Redistributions in binary form must reproduce the above copyright
18 notice, this list of conditions and the following disclaimer in the
19 documentation and/or other materials provided with the distribution.
20
21 THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS “AS
22 IS” AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED
23 TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A
24 PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE REGENTS OR
25 CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
26 EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
27 PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
28 PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
29 LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
30 NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
31 SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
32 *//** @file
33 This code computes the numerical solution of a system of first order ordinary
34 differential equations y'=f(x,y). It uses an explicit Runge-Kutta method of
35 order (4)5 due to Dormand & Prince with step size control and dense output.
36
37 Authors : E. Hairer & G. Wanner
38 Universite de Geneve, dept. de Mathematiques
39 CH-1211 GENEVE 4, SWITZERLAND
40 E-mail : HAIRER@DIVSUN.UNIGE.CH, WANNER@DIVSUN.UNIGE.CH
41
42 The code is described in : E. Hairer, S.P. Norsett and G. Wanner, Solving
43 ordinary differential equations I, nonstiff problems, 2nd edition,
44 Springer Series in Computational Mathematics, Springer-Verlag (1993).
45
46 Version of April 28, 1994.
47
48 Remarks about the C version : this version allocates memory by itself, the
49 iwork array (among the initial FORTRAN parameters) has been split into
50 independant initial parameters, the statistical variables and last step size
51 and x have been encapsulated in the module and are now accessible through
52 dedicated functions, the variable names have been kept to maintain a kind
53 of reading compatibility between the C and FORTRAN codes; adaptation made by
54 J.Colinge (COLINGE@DIVSUN.UNIGE.CH).
55
56 INPUT PARAMETERS
57 ----------------
58
59 n Dimension of the system (n < UINT_MAX).
60
61 fcn A pointer the the function definig the differential equation, this
62 function must have the following prototype
63
64 void fcn (unsigned n, double x, double *y, double *f)
65
66 where the array f will be filled with the function result.
67
68 x Initial x value.
69
70 *y Initial y values (double y[n]).
71
72 xend Final x value (xend-x may be positive or negative).
73
74 *rtoler Relative and absolute error tolerances. They can be both scalars or
75 *atoler vectors of length n (in the scalar case pass the addresses of
76 variables where you have placed the tolerance values).
77
78 itoler Switch for atoler and rtoler :
79 itoler=0 : both atoler and rtoler are scalars, the code keeps
80 roughly the local error of y[i] below
81 rtoler*abs(y[i])+atoler.
82 itoler=1 : both rtoler and atoler are vectors, the code keeps
83 the local error of y[i] below
84 rtoler[i]*abs(y[i])+atoler[i].
85
86 solout A pointer to the output function called during integration.
87 If iout >= 1, it is called after every successful step. If iout = 0,
88 pass a pointer equal to NULL. solout must must have the following
89 prototype
90
91 solout (long nr, double xold, double x, double* y, unsigned n, int* irtrn)
92
93 where y is the solution the at nr-th grid point x, xold is the
94 previous grid point and irtrn serves to interrupt the integration
95 (if set to a negative value).
96
97 Continuous output : during the calls to solout, a continuous solution
98 for the interval (xold,x) is available through the function
99
100 contd5(i,s)
101
102 which provides an approximation to the i-th component of the solution
103 at the point s (s must lie in the interval (xold,x)).
104
105 iout Switch for calling solout :
106 iout=0 : no call,
107 iout=1 : solout only used for output,
108 iout=2 : dense output is performed in solout (in this case nrdens
109 must be greater than 0).
110
111 fileout A pointer to the stream used for messages, if you do not want any
112 message, just pass NULL.
113
114 icont An array containing the indexes of components for which dense
115 output is required. If no dense output is required, pass NULL.
116
117 licont The number of cells in icont.
118
119
120 Sophisticated setting of parameters
121 -----------------------------------
122
123 Several parameters have a default value (if set to 0) but, to better
124 adapt the code to your problem, you can specify particular initial
125 values.
126
127 uround The rounding unit, default 2.3E-16 (this default value can be
128 replaced in the code by DBL_EPSILON providing float.h defines it
129 in your system).
130
131 safe Safety factor in the step size prediction, default 0.9.
132
133 fac1 Parameters for step size selection; the new step size is chosen
134 fac2 subject to the restriction fac1 <= hnew/hold <= fac2.
135 Default values are fac1=0.2 and fac2=10.0.
136
137 beta The "beta" for stabilized step size control (see section IV.2 of our
138 book). Larger values for beta ( <= 0.1 ) make the step size control
139 more stable. dopri5 needs a larger beta than Higham & Hall. Negative
140 initial value provoke beta=0; default beta=0.04.
141
142 hmax Maximal step size, default xend-x.
143
144 h Initial step size, default is a guess computed by the function hinit.
145
146 nmax Maximal number of allowed steps, default 100000.
147
148 meth Switch for the choice of the method coefficients; at the moment the
149 only possibility and default value are 1.
150
151 nstiff Test for stiffness is activated when the current step number is a
152 multiple of nstiff. A negative value means no test and the default
153 is 1000.
154
155 nrdens Number of components for which dense outpout is required, default 0.
156 For 0 < nrdens < n, the components have to be specified in icont[0],
157 icont[1], ... icont[nrdens-1]. Note that if nrdens=0 or nrdens=n, no
158 icont is needed, pass NULL.
159
160
161 Memory requirements
162 -------------------
163
164 The function dopri5 allocates dynamically 8*n doubles for the method
165 stages, 5*nrdens doubles for the interpolation if dense output is
166 performed and n unsigned if 0 < nrdens < n.
167
168
169
170 OUTPUT PARAMETERS
171 -----------------
172
173 y numerical solution at x=xRead() (see below).
174
175 dopri5 returns the following values
176
177 1 : computation successful,
178 2 : computation successful interrupted by solout,
179 -1 : input is not consistent,
180 -2 : larger nmax is needed,
181 -3 : step size becomes too small,
182 -4 : the problem is probably stff (interrupted).
183
184
185 Several functions provide access to different values :
186
187 xRead x value for which the solution has been computed (x=xend after
188 successful return).
189
190 hRead Predicted step size of the last accepted step (useful for a
191 subsequent call to dopri5).
192
193 nstepRead Number of used steps.
194 naccptRead Number of accepted steps.
195 nrejctRead Number of rejected steps.
196 nfcnRead Number of function calls.
197 */
198
199 #include <stdio.h>
200 #include <limits.h>
201
202 typedef void FcnEqDiff(unsigned n, double x, double *y, double *f, void *user_data);
203 typedef void SolTrait(long nr, double xold, double x, double* y, unsigned n, int* irtrn, void *user_data);
204
205 extern int dopri5(
206 unsigned n, /* dimension of the system <= UINT_MAX-1*/
207 FcnEqDiff *fcn, /* function computing the value of f(x,y) */
208 double x, /* initial x-value */
209 double* y, /* initial values for y */
210 double xend, /* final x-value (xend-x may be positive or negative) */
211 double* rtoler, /* relative error tolerance */
212 double* atoler, /* absolute error tolerance */
213 int itoler, /* switch for rtoler and atoler */
214 SolTrait *solout, /* function providing the numerical solution during integration */
215 int iout, /* switch for calling solout */
216 FILE* fileout, /* messages stream */
217 double uround, /* rounding unit */
218 double safe, /* safety factor */
219 double fac1, /* parameters for step size selection */
220 double fac2,
221 double beta, /* for stabilized step size control */
222 double hmax, /* maximal step size */
223 double h, /* initial step size */
224 long nmax, /* maximal number of allowed steps */
225 int meth, /* switch for the choice of the coefficients */
226 long nstiff, /* test for stiffness */
227 unsigned nrdens, /* number of components for which dense outpout is required */
228 unsigned* icont, /* indexes of components for which dense output is required, >= nrdens */
229 unsigned licont /* declared length of icon */
230 ,void *user_data /* passed to solout and fcn */
231 );
232
233 extern double contd5(
234 unsigned ii, /* index of desired component */
235 double x /* approximation at x */
236 );
237
238 extern long nfcnRead (void); /* encapsulation of statistical data */
239 extern long nstepRead (void);
240 extern long naccptRead (void);
241 extern long nrejctRead (void);
242 extern double hRead (void);
243 extern double xRead (void);
244

john.pye@anu.edu.au
ViewVC Help
Powered by ViewVC 1.1.22