Parent Directory | Revision Log

Revision **1445** -
(**show annotations**)
(**download**)
(**as text**)

*Sat May 26 14:43:48 2007 UTC*
(17 years, 4 months ago)
by *jpye*

File MIME type: text/x-chdr

File size: 9724 byte(s)

File MIME type: text/x-chdr

File size: 9724 byte(s)

Starting some work on a DOPRI5 integrator.

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 splitted 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); |

203 | typedef void SolTrait(long nr, double xold, double x, double* y, unsigned n, int* irtrn); |

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 | ); |

231 | |

232 | extern double contd5( |

233 | unsigned ii, /* index of desired component */ |

234 | double x /* approximation at x */ |

235 | ); |

236 | |

237 | extern long nfcnRead (void); /* encapsulation of statistical data */ |

238 | extern long nstepRead (void); |

239 | extern long naccptRead (void); |

240 | extern long nrejctRead (void); |

241 | extern double hRead (void); |

242 | extern double xRead (void); |

243 |

john.pye@anu.edu.au | ViewVC Help |

Powered by ViewVC 1.1.22 |