| Geant4 Cross Reference |
1 // 1
2 // *******************************************
3 // * License and Disclaimer
4 // *
5 // * The Geant4 software is copyright of th
6 // * the Geant4 Collaboration. It is provided
7 // * conditions of the Geant4 Software License
8 // * LICENSE and available at http://cern.ch/
9 // * include a list of copyright holders.
10 // *
11 // * Neither the authors of this software syst
12 // * institutes,nor the agencies providing fin
13 // * work make any representation or warran
14 // * regarding this software system or assum
15 // * use. Please see the license in the file
16 // * for the full disclaimer and the limitatio
17 // *
18 // * This code implementation is the result
19 // * technical work of the GEANT4 collaboratio
20 // * By using, copying, modifying or distri
21 // * any work based on the software) you ag
22 // * use in resulting scientific publicati
23 // * acceptance of all terms of the Geant4 Sof
24 // *******************************************
25 //
26 // G4DormandPrince745 implementation
27 //
28 // DormandPrince7 - 5(4) non-FSAL
29 // definition of the stepper() method that eva
30 // field propagation.
31 // The coefficients and the algorithm have bee
32 //
33 // J. R. Dormand and P. J. Prince, "A family o
34 // Journal of computational and applied Math.,
35 //
36 // The Butcher table of the Dormand-Prince-7-4
37 //
38 // 0 |
39 // 1/5 | 1/5
40 // 3/10| 3/40 9/40
41 // 4/5 | 44/45 56/15 32/9
42 // 8/9 | 19372/6561 25360/2187 64448/6561
43 // 1 | 9017/3168 355/33 46732/5247
44 // 1 | 35/384 0 500/1113
45 // ----------------------------------------
46 // 35/384 0 500/1113
47 // 5179/57600 0 7571/16695
48 //
49 // Created: Somnath Banerjee, Google Summer of
50 // Supervision: John Apostolakis, CERN
51 // -------------------------------------------
52
53 #include "G4DormandPrince745.hh"
54 #include "G4LineSection.hh"
55
56 #include <cstring>
57
58 using namespace field_utils;
59
60 // Name of this steppers
61 const G4String& G4DormandPrince745::StepperTyp
62 {
63 static G4String _stepperType("G4DormandPrinc
64 return _stepperType;
65 }
66
67 // Description of this steppers - plus details
68 const G4String& G4DormandPrince745::StepperDes
69 {
70 static G4String _stepperDescription(
71 "Embedeed 5th order Runge-Kutta stepper -
72 return _stepperDescription;
73 }
74
75 G4DormandPrince745::G4DormandPrince745(G4Equat
76 G4int n
77 : G4MagIntegratorStepper(equation, noInteg
78 {
79 }
80
81 void G4DormandPrince745::Stepper(const G4doubl
82 const G4doubl
83 G4doubl
84 G4doubl
85 G4doubl
86 G4doubl
87 {
88 Stepper(yInput, dydx, hstep, yOutput, yError
89 field_utils::copy(dydxOutput, ak7);
90 }
91
92 // Stepper
93 //
94 // Passing in the value of yInput[],the first
95 // Giving back yOut and yErr arrays for output
96 //
97 void G4DormandPrince745::Stepper(const G4doubl
98 const G4doubl
99 G4doubl
100 G4doubl
101 G4doubl
102 {
103 // The various constants defined on the ba
104 //
105 const G4double b21 = 0.2,
106 b31 = 3.0 / 40.0, b32 = 9.0
107 b41 = 44.0 / 45.0, b42 = -5
108
109 b51 = 19372.0 / 6561.0, b52 = -25360.0
110 b54 = -212.0 / 729.0,
111
112 b61 = 9017.0 / 3168.0 , b62 = -355.0
113 b63 = 46732.0 / 5247.0, b64 = 49.0 / 1
114 b65 = -5103.0 / 18656.0,
115
116 b71 = 35.0 / 384.0, b72 = 0.,
117 b73 = 500.0 / 1113.0, b74 = 125.0 / 19
118 b75 = -2187.0 / 6784.0, b76 = 11.0 / 8
119
120 //Sum of columns, sum(bij) = ei
121 // e1 = 0. ,
122 // e2 = 1.0/5.0 ,
123 // e3 = 3.0/10.0 ,
124 // e4 = 4.0/5.0 ,
125 // e5 = 8.0/9.0 ,
126 // e6 = 1.0 ,
127 // e7 = 1.0 ,
128
129 // Difference between the higher and the l
130 // b7j are the coefficients of higher orde
131
132 dc1 = -(b71 - 5179.0 / 57600.0),
133 dc2 = -(b72 - .0),
134 dc3 = -(b73 - 7571.0 / 16695.0),
135 dc4 = -(b74 - 393.0 / 640.0),
136 dc5 = -(b75 + 92097.0 / 339200.0),
137 dc6 = -(b76 - 187.0 / 2100.0),
138 dc7 = -(- 1.0 / 40.0);
139
140 const G4int numberOfVariables = GetNumberO
141 State yTemp = {0., 0., 0., 0., 0., 0., 0.,
142
143 // The number of variables to be integrate
144 //
145 yOut[7] = yTemp[7] = yInput[7];
146
147 // Saving yInput because yInput and yOut
148 //
149 for(G4int i = 0; i < numberOfVariables; ++
150 {
151 fyIn[i] = yInput[i];
152 }
153 // RightHandSide(yIn, dydx); // Not don
154
155 for(G4int i = 0; i < numberOfVariables; ++
156 {
157 yTemp[i] = fyIn[i] + b21 * hstep * dyd
158 }
159 RightHandSide(yTemp, ak2); //
160
161 for(G4int i = 0; i < numberOfVariables; ++
162 {
163 yTemp[i] = fyIn[i] + hstep * (b31 * dy
164 }
165 RightHandSide(yTemp, ak3); //
166
167 for(G4int i = 0; i < numberOfVariables; ++
168 {
169 yTemp[i] = fyIn[i] + hstep * (
170 b41 * dydx[i] + b42 * ak2[i] + b43
171 }
172 RightHandSide(yTemp, ak4); //
173
174 for(G4int i = 0; i < numberOfVariables; ++
175 {
176 yTemp[i] = fyIn[i] + hstep * (
177 b51 * dydx[i] + b52 * ak2[i] + b53
178 }
179 RightHandSide(yTemp, ak5); //
180
181 for(G4int i = 0; i < numberOfVariables; ++
182 {
183 yTemp[i] = fyIn[i] + hstep * (
184 b61 * dydx[i] + b62 * ak2[i] +
185 b63 * ak3[i] + b64 * ak4[i] + b65
186 }
187 RightHandSide(yTemp, ak6); //
188
189 for(G4int i = 0; i < numberOfVariables; ++
190 {
191 yOut[i] = fyIn[i] + hstep * (
192 b71 * dydx[i] + b72 * ak2[i] + b73
193 b74 * ak4[i] + b75 * ak5[i] + b76
194 }
195 RightHandSide(yOut, ak7); //
196
197 for(G4int i = 0; i < numberOfVariables; ++
198 {
199 yErr[i] = hstep * (
200 dc1 * dydx[i] + dc2 * ak2[i] +
201 dc3 * ak3[i] + dc4 * ak4[i] +
202 dc5 * ak5[i] + dc6 * ak6[i] + dc7
203 ) + 1.5e-18;
204
205 // Store Input and Final values, for p
206 //
207 fyOut[i] = yOut[i];
208 fdydxIn[i] = dydx[i];
209 }
210
211 fLastStepLength = hstep;
212 }
213
214 G4double G4DormandPrince745::DistChord() const
215 {
216 // Coefficients were taken from Some Pract
217 // by Lawrence F. Shampine, page 149, c*
218 //
219 const G4double hf1 = 6025192743.0 / 300855
220 hf3 = 51252292925.0 / 65400
221 hf4 = - 2691868925.0 / 4512
222 hf5 = 187940372067.0 / 1594
223 hf6 = - 1776094331.0 / 1974
224 hf7 = 11237099.0 / 23504338
225
226 G4ThreeVector mid;
227
228 for(G4int i = 0; i < 3; ++i)
229 {
230 mid[i] = fyIn[i] + 0.5 * fLastStepLeng
231 hf1 * fdydxIn[i] + hf3 * ak3[i] +
232 hf4 * ak4[i] + hf5 * ak5[i] + hf6
233 }
234
235 const G4ThreeVector begin = makeVector(fyI
236 const G4ThreeVector end = makeVector(fyOut
237
238 return G4LineSection::Distline(mid, begin,
239 }
240
241 // The lower (4th) order interpolant given by
242 // J. R. Dormand and P. J. Prince, "Run
243 // Computers & Mathematics with Applica
244 // pp. 1007-1017, 1986.
245 //
246 void G4DormandPrince745::
247 Interpolate4thOrder(G4double yOut[], G4double
248 {
249 const G4int numberOfVariables = GetNumberO
250
251 const G4double tau2 = tau * tau,
252 tau3 = tau * tau2,
253 tau4 = tau2 * tau2;
254
255 const G4double bf1 = 1.0 / 11282082432.0 *
256 157015080.0 * tau4 - 13107642775.0 * t
257 32272833064.0 * tau + 11282082432.0);
258
259 const G4double bf3 = - 100.0 / 32700410799
260 15701508.0 * tau3 - 914128567.0 * tau2
261 1323431896.0);
262
263 const G4double bf4 = 25.0 / 5641041216.0 *
264 94209048.0 * tau3 - 1518414297.0 * tau
265 889289856.0);
266
267 const G4double bf5 = - 2187.0 / 1993167896
268 52338360.0 * tau3 - 451824525.0 * tau2
269 259006536.0);
270
271 const G4double bf6 = 11.0 / 2467955532.0 *
272 106151040.0 * tau3 - 661884105.0 * tau
273 946554244.0 * tau - 361440756.0);
274
275 const G4double bf7 = 1.0 / 29380423.0 * ta
276 8293050.0 * tau2 - 82437520.0 * tau +
277
278 for(G4int i = 0; i < numberOfVariables; ++
279 {
280 yOut[i] = fyIn[i] + fLastStepLength *
281 bf1 * fdydxIn[i] + bf3 * ak3[i] +
282 bf5 * ak5[i] + bf6 * ak6[i] + bf7
283 }
284 }
285
286 // Following interpolant of order 5 was given
287 // T. S. Baker, J. R. Dormand, J. P. Gi
288 // "Continuous approximation with embed
289 // Applied Numerical Mathematics, vol.
290 //
291 // Calculating the extra stages for the interp
292 //
293 void G4DormandPrince745::SetupInterpolation5th
294 {
295 // Coefficients for the additional stages
296 //
297 const G4double b81 = 6245.0 / 62208.0,
298 b82 = 0.0,
299 b83 = 8875.0 / 103032.0,
300 b84 = -125.0 / 1728.0,
301 b85 = 801.0 / 13568.0,
302 b86 = -13519.0 / 368064.0,
303 b87 = 11105.0 / 368064.0,
304
305 b91 = 632855.0 / 4478976.0,
306 b92 = 0.0,
307 b93 = 4146875.0 / 6491016.0
308 b94 = 5490625.0 /14183424.0
309 b95 = -15975.0 / 108544.0,
310 b96 = 8295925.0 / 220286304
311 b97 = -1779595.0 / 62938944
312 b98 = -805.0 / 4104.0;
313
314 const G4int numberOfVariables = GetNumberO
315 State yTemp = {0., 0., 0., 0., 0., 0., 0.,
316
317 // Evaluate the extra stages
318 //
319 for(G4int i = 0; i < numberOfVariables; ++
320 {
321 yTemp[i] = fyIn[i] + fLastStepLength *
322 b81 * fdydxIn[i] + b82 * ak2[i] +
323 b84 * ak4[i] + b85 * ak5[i] + b86
324 b87 * ak7[i]
325 );
326 }
327 RightHandSide(yTemp, ak8); // 8th
328
329 for(G4int i = 0; i < numberOfVariables; ++
330 {
331 yTemp[i] = fyIn[i] + fLastStepLength *
332 b91 * fdydxIn[i] + b92 * ak2[i] +
333 b94 * ak4[i] + b95 * ak5[i] + b96
334 b97 * ak7[i] + b98 * ak8[i]
335 );
336 }
337 RightHandSide(yTemp, ak9); // 9th
338 }
339
340 // Calculating the interpolated result yOut wi
341 //
342 void G4DormandPrince745::
343 Interpolate5thOrder(G4double yOut[], G4double
344 {
345 // Define the coefficients for the polynom
346 //
347 G4double bi[10][5];
348
349 // COEFFICIENTS OF bi[1]
350 bi[1][0] = 1.0,
351 bi[1][1] = -38039.0 / 7040.0,
352 bi[1][2] = 125923.0 / 10560.0,
353 bi[1][3] = -19683.0 / 1760.0,
354 bi[1][4] = 3303.0 / 880.0,
355 // --------------------------------------
356 //
357 // COEFFICIENTS OF bi[2]
358 bi[2][0] = 0.0,
359 bi[2][1] = 0.0,
360 bi[2][2] = 0.0,
361 bi[2][3] = 0.0,
362 bi[2][4] = 0.0,
363 // --------------------------------------
364 //
365 // COEFFICIENTS OF bi[3]
366 bi[3][0] = 0.0,
367 bi[3][1] = -12500.0 / 4081.0,
368 bi[3][2] = 205000.0 / 12243.0,
369 bi[3][3] = -90000.0 / 4081.0,
370 bi[3][4] = 36000.0 / 4081.0,
371 // --------------------------------------
372 //
373 // COEFFICIENTS OF bi[4]
374 bi[4][0] = 0.0,
375 bi[4][1] = -3125.0 / 704.0,
376 bi[4][2] = 25625.0 / 1056.0,
377 bi[4][3] = -5625.0 / 176.0,
378 bi[4][4] = 1125.0 / 88.0,
379 // --------------------------------------
380 //
381 // COEFFICIENTS OF bi[5]
382 bi[5][0] = 0.0,
383 bi[5][1] = 164025.0 / 74624.0,
384 bi[5][2] = -448335.0 / 37312.0,
385 bi[5][3] = 295245.0 / 18656.0,
386 bi[5][4] = -59049.0 / 9328.0,
387 // --------------------------------------
388 //
389 // COEFFICIENTS OF bi[6]
390 bi[6][0] = 0.0,
391 bi[6][1] = -25.0 / 28.0,
392 bi[6][2] = 205.0 / 42.0,
393 bi[6][3] = -45.0 / 7.0,
394 bi[6][4] = 18.0 / 7.0,
395 // --------------------------------------
396 //
397 // COEFFICIENTS OF bi[7]
398 bi[7][0] = 0.0,
399 bi[7][1] = -2.0 / 11.0,
400 bi[7][2] = 73.0 / 55.0,
401 bi[7][3] = -171.0 / 55.0,
402 bi[7][4] = 108.0 / 55.0,
403 // --------------------------------------
404 //
405 // COEFFICIENTS OF bi[8]
406 bi[8][0] = 0.0,
407 bi[8][1] = 189.0 / 22.0,
408 bi[8][2] = -1593.0 / 55.0,
409 bi[8][3] = 3537.0 / 110.0,
410 bi[8][4] = -648.0 / 55.0,
411 // --------------------------------------
412 //
413 // COEFFICIENTS OF bi[9]
414 bi[9][0] = 0.0,
415 bi[9][1] = 351.0 / 110.0,
416 bi[9][2] = -999.0 / 55.0,
417 bi[9][3] = 2943.0 / 110.0,
418 bi[9][4] = -648.0 / 55.0;
419 // --------------------------------------
420
421 // Calculating the polynomials
422
423 G4double b[10];
424 std::memset(b, 0.0, sizeof(b));
425
426 G4double tauPower = 1.0;
427 for(G4int j = 0; j <= 4; ++j)
428 {
429 for(G4int iStage = 1; iStage <= 9; ++iS
430 {
431 b[iStage] += bi[iStage][j] * tauPowe
432 }
433 tauPower *= tau;
434 }
435
436 const G4int numberOfVariables = GetNumberO
437 const G4double stepLen = fLastStepLength *
438 for(G4int i = 0; i < numberOfVariables; ++
439 {
440 yOut[i] = fyIn[i] + stepLen * (
441 b[1] * fdydxIn[i] + b[2] * ak2[i]
442 b[4] * ak4[i] + b[5] * ak5[i] + b[
443 b[7] * ak7[i] + b[8] * ak8[i] + b[
444 );
445 }
446 }
447