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