1 /*
2 * Licensed to the Hipparchus project under one or more
3 * contributor license agreements. See the NOTICE file distributed with
4 * this work for additional information regarding copyright ownership.
5 * The Hipparchus project licenses this file to You under the Apache License, Version 2.0
6 * (the "License"); you may not use this file except in compliance with
7 * the License. You may obtain a copy of the License at
8 *
9 * https://www.apache.org/licenses/LICENSE-2.0
10 *
11 * Unless required by applicable law or agreed to in writing, software
12 * distributed under the License is distributed on an "AS IS" BASIS,
13 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
14 * See the License for the specific language governing permissions and
15 * limitations under the License.
16 */
17
18 package org.hipparchus.ode.nonstiff;
19
20 import org.hipparchus.ode.EquationsMapper;
21 import org.hipparchus.ode.ODEStateAndDerivative;
22 import org.hipparchus.util.FastMath;
23
24
25 /**
26 * This class implements the 8(5,3) Dormand-Prince integrator for Ordinary
27 * Differential Equations.
28 *
29 * <p>This integrator is an embedded Runge-Kutta integrator
30 * of order 8(5,3) used in local extrapolation mode (i.e. the solution
31 * is computed using the high order formula) with stepsize control
32 * (and automatic step initialization) and continuous output. This
33 * method uses 12 functions evaluations per step for integration and 4
34 * evaluations for interpolation. However, since the first
35 * interpolation evaluation is the same as the first integration
36 * evaluation of the next step, we have included it in the integrator
37 * rather than in the interpolator and specified the method was an
38 * <i>fsal</i>. Hence, despite we have 13 stages here, the cost is
39 * really 12 evaluations per step even if no interpolation is done,
40 * and the overcost of interpolation is only 3 evaluations.</p>
41 *
42 * <p>This method is based on an 8(6) method by Dormand and Prince
43 * (i.e. order 8 for the integration and order 6 for error estimation)
44 * modified by Hairer and Wanner to use a 5th order error estimator
45 * with 3rd order correction. This modification was introduced because
46 * the original method failed in some cases (wrong steps can be
47 * accepted when step size is too large, for example in the
48 * Brusselator problem) and also had <i>severe difficulties when
49 * applied to problems with discontinuities</i>. This modification is
50 * explained in the second edition of the first volume (Nonstiff
51 * Problems) of the reference book by Hairer, Norsett and Wanner:
52 * <i>Solving Ordinary Differential Equations</i> (Springer-Verlag,
53 * ISBN 3-540-56670-8).</p>
54 *
55 */
56
57 public class DormandPrince853Integrator extends EmbeddedRungeKuttaIntegrator {
58
59 /** Name of integration scheme. */
60 public static final String METHOD_NAME = "Dormand-Prince 8 (5, 3)";
61
62 /** First error weights array, element 1. */
63 static final double E1_01 = 116092271.0 / 8848465920.0;
64
65 // elements 2 to 5 are zero, so they are neither stored nor used
66
67 /** First error weights array, element 6. */
68 static final double E1_06 = -1871647.0 / 1527680.0;
69
70 /** First error weights array, element 7. */
71 static final double E1_07 = -69799717.0 / 140793660.0;
72
73 /** First error weights array, element 8. */
74 static final double E1_08 = 1230164450203.0 / 739113984000.0;
75
76 /** First error weights array, element 9. */
77 static final double E1_09 = -1980813971228885.0 / 5654156025964544.0;
78
79 /** First error weights array, element 10. */
80 static final double E1_10 = 464500805.0 / 1389975552.0;
81
82 /** First error weights array, element 11. */
83 static final double E1_11 = 1606764981773.0 / 19613062656000.0;
84
85 /** First error weights array, element 12. */
86 static final double E1_12 = -137909.0 / 6168960.0;
87
88
89 /** Second error weights array, element 1. */
90 static final double E2_01 = -364463.0 / 1920240.0;
91
92 // elements 2 to 5 are zero, so they are neither stored nor used
93
94 /** Second error weights array, element 6. */
95 static final double E2_06 = 3399327.0 / 763840.0;
96
97 /** Second error weights array, element 7. */
98 static final double E2_07 = 66578432.0 / 35198415.0;
99
100 /** Second error weights array, element 8. */
101 static final double E2_08 = -1674902723.0 / 288716400.0;
102
103 /** Second error weights array, element 9. */
104 static final double E2_09 = -74684743568175.0 / 176692375811392.0;
105
106 /** Second error weights array, element 10. */
107 static final double E2_10 = -734375.0 / 4826304.0;
108
109 /** Second error weights array, element 11. */
110 static final double E2_11 = 171414593.0 / 851261400.0;
111
112 /** Second error weights array, element 12. */
113 static final double E2_12 = 69869.0 / 3084480.0;
114
115 /** Simple constructor.
116 * Build a fifth order Dormand-Prince integrator with the given step bounds
117 * @param minStep minimal step (sign is irrelevant, regardless of
118 * integration direction, forward or backward), the last step can
119 * be smaller than this
120 * @param maxStep maximal step (sign is irrelevant, regardless of
121 * integration direction, forward or backward), the last step can
122 * be smaller than this
123 * @param scalAbsoluteTolerance allowed absolute error
124 * @param scalRelativeTolerance allowed relative error
125 */
126 public DormandPrince853Integrator(final double minStep, final double maxStep,
127 final double scalAbsoluteTolerance,
128 final double scalRelativeTolerance) {
129 super(METHOD_NAME, 12,
130 minStep, maxStep, scalAbsoluteTolerance, scalRelativeTolerance);
131 }
132
133 /** Simple constructor.
134 * Build a fifth order Dormand-Prince integrator with the given step bounds
135 * @param minStep minimal step (sign is irrelevant, regardless of
136 * integration direction, forward or backward), the last step can
137 * be smaller than this
138 * @param maxStep maximal step (sign is irrelevant, regardless of
139 * integration direction, forward or backward), the last step can
140 * be smaller than this
141 * @param vecAbsoluteTolerance allowed absolute error
142 * @param vecRelativeTolerance allowed relative error
143 */
144 public DormandPrince853Integrator(final double minStep, final double maxStep,
145 final double[] vecAbsoluteTolerance,
146 final double[] vecRelativeTolerance) {
147 super(METHOD_NAME, 12,
148 minStep, maxStep, vecAbsoluteTolerance, vecRelativeTolerance);
149 }
150
151 /** {@inheritDoc} */
152 @Override
153 public double[] getC() {
154 final double sqrt6 = FastMath.sqrt(6.0);
155 return new double[] {
156 (12.0 - 2.0 * sqrt6) / 135.0,
157 (6.0 - sqrt6) / 45.0,
158 (6.0 - sqrt6) / 30.0,
159 (6.0 + sqrt6) / 30.0,
160 1.0/3.0,
161 1.0/4.0,
162 4.0/13.0,
163 127.0/195.0,
164 3.0/5.0,
165 6.0/7.0,
166 1.0,
167 1.0,
168 1.0/10.0,
169 1.0/5.0,
170 7.0/9.0
171 };
172 }
173
174 /** {@inheritDoc} */
175 @Override
176 public double[][] getA() {
177 final double sqrt6 = FastMath.sqrt(6.0);
178 return new double[][] {
179 {
180 (12.0 - 2.0 * sqrt6) / 135.0
181 }, {
182 (6.0 - sqrt6) / 180.0,
183 (6.0 - sqrt6) / 60.0
184 }, {
185 (6.0 - sqrt6) / 120.0,
186 0.0,
187 (6.0 - sqrt6) / 40.0
188 }, {
189 (462.0 + 107.0 * sqrt6) / 3000.0,
190 0.0,
191 (-402.0 - 197.0 * sqrt6) / 1000.0,
192 (168.0 + 73.0 * sqrt6) / 375.0
193 }, {
194 1.0 / 27.0,
195 0.0,
196 0.0,
197 (16.0 + sqrt6) / 108.0,
198 (16.0 - sqrt6) / 108.0
199 }, {
200 19.0 / 512.0,
201 0.0,
202 0.0,
203 (118.0 + 23.0 * sqrt6) / 1024.0,
204 (118.0 - 23.0 * sqrt6) / 1024.0,
205 -9.0 / 512.0
206 }, {
207 13772.0 / 371293.0,
208 0.0,
209 0.0,
210 (51544.0 + 4784.0 * sqrt6) / 371293.0,
211 (51544.0 - 4784.0 * sqrt6) / 371293.0,
212 -5688.0 / 371293.0,
213 3072.0 / 371293.0
214 }, {
215 58656157643.0 / 93983540625.0,
216 0.0,
217 0.0,
218 (-1324889724104.0 - 318801444819.0 * sqrt6) / 626556937500.0,
219 (-1324889724104.0 + 318801444819.0 * sqrt6) / 626556937500.0,
220 96044563816.0 / 3480871875.0,
221 5682451879168.0 / 281950621875.0,
222 -165125654.0 / 3796875.0
223 }, {
224 8909899.0 / 18653125.0,
225 0.0,
226 0.0,
227 (-4521408.0 - 1137963.0 * sqrt6) / 2937500.0,
228 (-4521408.0 + 1137963.0 * sqrt6) / 2937500.0,
229 96663078.0 / 4553125.0,
230 2107245056.0 / 137915625.0,
231 -4913652016.0 / 147609375.0,
232 -78894270.0 / 3880452869.0
233 }, {
234 -20401265806.0 / 21769653311.0,
235 0.0,
236 0.0,
237 (354216.0 + 94326.0 * sqrt6) / 112847.0,
238 (354216.0 - 94326.0 * sqrt6) / 112847.0,
239 -43306765128.0 / 5313852383.0,
240 -20866708358144.0 / 1126708119789.0,
241 14886003438020.0 / 654632330667.0,
242 35290686222309375.0 / 14152473387134411.0,
243 -1477884375.0 / 485066827.0
244 }, {
245 39815761.0 / 17514443.0,
246 0.0,
247 0.0,
248 (-3457480.0 - 960905.0 * sqrt6) / 551636.0,
249 (-3457480.0 + 960905.0 * sqrt6) / 551636.0,
250 -844554132.0 / 47026969.0,
251 8444996352.0 / 302158619.0,
252 -2509602342.0 / 877790785.0,
253 -28388795297996250.0 / 3199510091356783.0,
254 226716250.0 / 18341897.0,
255 1371316744.0 / 2131383595.0
256 }, {
257 // the following stage is both for interpolation and the first stage in next step
258 // (the coefficients are identical to the B array)
259 104257.0/1920240.0,
260 0.0,
261 0.0,
262 0.0,
263 0.0,
264 3399327.0/763840.0,
265 66578432.0/35198415.0,
266 -1674902723.0/288716400.0,
267 54980371265625.0/176692375811392.0,
268 -734375.0/4826304.0,
269 171414593.0/851261400.0,
270 137909.0/3084480.0
271 }, {
272 // the following stages are for interpolation only
273 13481885573.0 / 240030000000.0,
274 0.0,
275 0.0,
276 0.0,
277 0.0,
278 0.0,
279 139418837528.0 / 549975234375.0,
280 -11108320068443.0 / 45111937500000.0,
281 -1769651421925959.0 / 14249385146080000.0,
282 57799439.0 / 377055000.0,
283 793322643029.0 / 96734250000000.0,
284 1458939311.0 / 192780000000.0,
285 -4149.0 / 500000.0
286 }, {
287 1595561272731.0 / 50120273500000.0,
288 0.0,
289 0.0,
290 0.0,
291 0.0,
292 975183916491.0 / 34457688031250.0,
293 38492013932672.0 / 718912673015625.0,
294 -1114881286517557.0 / 20298710767500000.0,
295 0.0,
296 0.0,
297 -2538710946863.0 / 23431227861250000.0,
298 8824659001.0 / 23066716781250.0,
299 -11518334563.0 / 33831184612500.0,
300 1912306948.0 / 13532473845.0
301 }, {
302 -13613986967.0 / 31741908048.0,
303 0.0,
304 0.0,
305 0.0,
306 0.0,
307 -4755612631.0 / 1012344804.0,
308 42939257944576.0 / 5588559685701.0,
309 77881972900277.0 / 19140370552944.0,
310 22719829234375.0 / 63689648654052.0,
311 0.0,
312 0.0,
313 0.0,
314 -1199007803.0 / 857031517296.0,
315 157882067000.0 / 53564469831.0,
316 -290468882375.0 / 31741908048.0
317 }
318 };
319 }
320
321 /** {@inheritDoc} */
322 @Override
323 public double[] getB() {
324 return new double[] {
325 104257.0/1920240.0,
326 0.0,
327 0.0,
328 0.0,
329 0.0,
330 3399327.0/763840.0,
331 66578432.0/35198415.0,
332 -1674902723.0/288716400.0,
333 54980371265625.0/176692375811392.0,
334 -734375.0/4826304.0,
335 171414593.0/851261400.0,
336 137909.0/3084480.0,
337 0.0,
338 0.0,
339 0.0,
340 0.0
341 };
342 }
343
344 /** {@inheritDoc} */
345 @Override
346 protected DormandPrince853StateInterpolator
347 createInterpolator(final boolean forward, double[][] yDotK,
348 final ODEStateAndDerivative globalPreviousState,
349 final ODEStateAndDerivative globalCurrentState,
350 final EquationsMapper mapper) {
351 return new DormandPrince853StateInterpolator(forward, yDotK,
352 globalPreviousState, globalCurrentState,
353 globalPreviousState, globalCurrentState,
354 mapper);
355 }
356
357 /** {@inheritDoc} */
358 @Override
359 public int getOrder() {
360 return 8;
361 }
362
363 /** {@inheritDoc} */
364 @Override
365 protected double estimateError(final double[][] yDotK,
366 final double[] y0, final double[] y1,
367 final double h) {
368
369 final StepsizeHelper helper = getStepSizeHelper();
370 double error1 = 0;
371 double error2 = 0;
372
373 for (int j = 0; j < helper.getMainSetDimension(); ++j) {
374 final double errSum1 = E1_01 * yDotK[0][j] + E1_06 * yDotK[5][j] +
375 E1_07 * yDotK[6][j] + E1_08 * yDotK[7][j] +
376 E1_09 * yDotK[8][j] + E1_10 * yDotK[9][j] +
377 E1_11 * yDotK[10][j] + E1_12 * yDotK[11][j];
378 final double errSum2 = E2_01 * yDotK[0][j] + E2_06 * yDotK[5][j] +
379 E2_07 * yDotK[6][j] + E2_08 * yDotK[7][j] +
380 E2_09 * yDotK[8][j] + E2_10 * yDotK[9][j] +
381 E2_11 * yDotK[10][j] + E2_12 * yDotK[11][j];
382
383 final double tol = helper.getTolerance(j, FastMath.max(FastMath.abs(y0[j]), FastMath.abs(y1[j])));
384 final double ratio1 = errSum1 / tol;
385 error1 += ratio1 * ratio1;
386 final double ratio2 = errSum2 / tol;
387 error2 += ratio2 * ratio2;
388 }
389
390 double den = error1 + 0.01 * error2;
391 if (den <= 0.0) {
392 den = 1.0;
393 }
394
395 return FastMath.abs(h) * error1 / FastMath.sqrt(helper.getMainSetDimension() * den);
396
397 }
398
399 }