001/*
002 * Licensed to the Apache Software Foundation (ASF) under one or more
003 * contributor license agreements.  See the NOTICE file distributed with
004 * this work for additional information regarding copyright ownership.
005 * The ASF licenses this file to You under the Apache License, Version 2.0
006 * (the "License"); you may not use this file except in compliance with
007 * the License.  You may obtain a copy of the License at
008 *
009 *      http://www.apache.org/licenses/LICENSE-2.0
010 *
011 * Unless required by applicable law or agreed to in writing, software
012 * distributed under the License is distributed on an "AS IS" BASIS,
013 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
014 * See the License for the specific language governing permissions and
015 * limitations under the License.
016 */
017
018package org.apache.commons.math4.ode.nonstiff;
019
020import org.apache.commons.math4.Field;
021import org.apache.commons.math4.RealFieldElement;
022import org.apache.commons.math4.ode.FieldEquationsMapper;
023import org.apache.commons.math4.ode.FieldODEStateAndDerivative;
024import org.apache.commons.math4.util.MathArrays;
025import org.apache.commons.math4.util.MathUtils;
026
027
028/**
029 * This class implements the 8(5,3) Dormand-Prince integrator for Ordinary
030 * Differential Equations.
031 *
032 * <p>This integrator is an embedded Runge-Kutta integrator
033 * of order 8(5,3) used in local extrapolation mode (i.e. the solution
034 * is computed using the high order formula) with stepsize control
035 * (and automatic step initialization) and continuous output. This
036 * method uses 12 functions evaluations per step for integration and 4
037 * evaluations for interpolation. However, since the first
038 * interpolation evaluation is the same as the first integration
039 * evaluation of the next step, we have included it in the integrator
040 * rather than in the interpolator and specified the method was an
041 * <i>fsal</i>. Hence, despite we have 13 stages here, the cost is
042 * really 12 evaluations per step even if no interpolation is done,
043 * and the overcost of interpolation is only 3 evaluations.</p>
044 *
045 * <p>This method is based on an 8(6) method by Dormand and Prince
046 * (i.e. order 8 for the integration and order 6 for error estimation)
047 * modified by Hairer and Wanner to use a 5th order error estimator
048 * with 3rd order correction. This modification was introduced because
049 * the original method failed in some cases (wrong steps can be
050 * accepted when step size is too large, for example in the
051 * Brusselator problem) and also had <i>severe difficulties when
052 * applied to problems with discontinuities</i>. This modification is
053 * explained in the second edition of the first volume (Nonstiff
054 * Problems) of the reference book by Hairer, Norsett and Wanner:
055 * <i>Solving Ordinary Differential Equations</i> (Springer-Verlag,
056 * ISBN 3-540-56670-8).</p>
057 *
058 * @param <T> the type of the field elements
059 * @since 3.6
060 */
061
062public class DormandPrince853FieldIntegrator<T extends RealFieldElement<T>>
063    extends EmbeddedRungeKuttaFieldIntegrator<T> {
064
065    /** Integrator method name. */
066    private static final String METHOD_NAME = "Dormand-Prince 8 (5, 3)";
067
068    /** First error weights array, element 1. */
069    private final T e1_01;
070
071    // elements 2 to 5 are zero, so they are neither stored nor used
072
073    /** First error weights array, element 6. */
074    private final T e1_06;
075
076    /** First error weights array, element 7. */
077    private final T e1_07;
078
079    /** First error weights array, element 8. */
080    private final T e1_08;
081
082    /** First error weights array, element 9. */
083    private final T e1_09;
084
085    /** First error weights array, element 10. */
086    private final T e1_10;
087
088    /** First error weights array, element 11. */
089    private final T e1_11;
090
091    /** First error weights array, element 12. */
092    private final T e1_12;
093
094
095    /** Second error weights array, element 1. */
096    private final T e2_01;
097
098    // elements 2 to 5 are zero, so they are neither stored nor used
099
100    /** Second error weights array, element 6. */
101    private final T e2_06;
102
103    /** Second error weights array, element 7. */
104    private final T e2_07;
105
106    /** Second error weights array, element 8. */
107    private final T e2_08;
108
109    /** Second error weights array, element 9. */
110    private final T e2_09;
111
112    /** Second error weights array, element 10. */
113    private final T e2_10;
114
115    /** Second error weights array, element 11. */
116    private final T e2_11;
117
118    /** Second error weights array, element 12. */
119    private final T e2_12;
120
121    /** Simple constructor.
122     * Build an eighth order Dormand-Prince integrator with the given step bounds
123     * @param field field to which the time and state vector elements belong
124     * @param minStep minimal step (sign is irrelevant, regardless of
125     * integration direction, forward or backward), the last step can
126     * be smaller than this
127     * @param maxStep maximal step (sign is irrelevant, regardless of
128     * integration direction, forward or backward), the last step can
129     * be smaller than this
130     * @param scalAbsoluteTolerance allowed absolute error
131     * @param scalRelativeTolerance allowed relative error
132     */
133    public DormandPrince853FieldIntegrator(final Field<T> field,
134                                           final double minStep, final double maxStep,
135                                           final double scalAbsoluteTolerance,
136                                           final double scalRelativeTolerance) {
137        super(field, METHOD_NAME, 12,
138              minStep, maxStep, scalAbsoluteTolerance, scalRelativeTolerance);
139        e1_01 = fraction(        116092271.0,       8848465920.0);
140        e1_06 = fraction(         -1871647.0,          1527680.0);
141        e1_07 = fraction(        -69799717.0,        140793660.0);
142        e1_08 = fraction(    1230164450203.0,     739113984000.0);
143        e1_09 = fraction(-1980813971228885.0, 5654156025964544.0);
144        e1_10 = fraction(        464500805.0,       1389975552.0);
145        e1_11 = fraction(    1606764981773.0,   19613062656000.0);
146        e1_12 = fraction(          -137909.0,          6168960.0);
147        e2_01 = fraction(          -364463.0,          1920240.0);
148        e2_06 = fraction(          3399327.0,           763840.0);
149        e2_07 = fraction(         66578432.0,         35198415.0);
150        e2_08 = fraction(      -1674902723.0,        288716400.0);
151        e2_09 = fraction(  -74684743568175.0,  176692375811392.0);
152        e2_10 = fraction(          -734375.0,          4826304.0);
153        e2_11 = fraction(        171414593.0,        851261400.0);
154        e2_12 = fraction(            69869.0,          3084480.0);
155    }
156
157    /** Simple constructor.
158     * Build an eighth order Dormand-Prince integrator with the given step bounds
159     * @param field field to which the time and state vector elements belong
160     * @param minStep minimal step (sign is irrelevant, regardless of
161     * integration direction, forward or backward), the last step can
162     * be smaller than this
163     * @param maxStep maximal step (sign is irrelevant, regardless of
164     * integration direction, forward or backward), the last step can
165     * be smaller than this
166     * @param vecAbsoluteTolerance allowed absolute error
167     * @param vecRelativeTolerance allowed relative error
168     */
169    public DormandPrince853FieldIntegrator(final Field<T> field,
170                                           final double minStep, final double maxStep,
171                                           final double[] vecAbsoluteTolerance,
172                                           final double[] vecRelativeTolerance) {
173        super(field, METHOD_NAME, 12,
174              minStep, maxStep, vecAbsoluteTolerance, vecRelativeTolerance);
175        e1_01 = fraction(        116092271.0,       8848465920.0);
176        e1_06 = fraction(         -1871647.0,          1527680.0);
177        e1_07 = fraction(        -69799717.0,        140793660.0);
178        e1_08 = fraction(    1230164450203.0,     739113984000.0);
179        e1_09 = fraction(-1980813971228885.0, 5654156025964544.0);
180        e1_10 = fraction(        464500805.0,       1389975552.0);
181        e1_11 = fraction(    1606764981773.0,   19613062656000.0);
182        e1_12 = fraction(          -137909.0,          6168960.0);
183        e2_01 = fraction(          -364463.0,          1920240.0);
184        e2_06 = fraction(          3399327.0,           763840.0);
185        e2_07 = fraction(         66578432.0,         35198415.0);
186        e2_08 = fraction(      -1674902723.0,        288716400.0);
187        e2_09 = fraction(  -74684743568175.0,  176692375811392.0);
188        e2_10 = fraction(          -734375.0,          4826304.0);
189        e2_11 = fraction(        171414593.0,        851261400.0);
190        e2_12 = fraction(            69869.0,          3084480.0);
191    }
192
193    /** {@inheritDoc} */
194    @Override
195    public T[] getC() {
196
197        final T sqrt6 = getField().getOne().multiply(6).sqrt();
198
199        final T[] c = MathArrays.buildArray(getField(), 15);
200        c[ 0] = sqrt6.add(-6).divide(-67.5);
201        c[ 1] = sqrt6.add(-6).divide(-45.0);
202        c[ 2] = sqrt6.add(-6).divide(-30.0);
203        c[ 3] = sqrt6.add( 6).divide( 30.0);
204        c[ 4] = fraction(1, 3);
205        c[ 5] = fraction(1, 4);
206        c[ 6] = fraction(4, 13);
207        c[ 7] = fraction(127, 195);
208        c[ 8] = fraction(3, 5);
209        c[ 9] = fraction(6, 7);
210        c[10] = getField().getOne();
211        c[11] = getField().getOne();
212        c[12] = fraction(1.0, 10.0);
213        c[13] = fraction(1.0, 5.0);
214        c[14] = fraction(7.0, 9.0);
215
216        return c;
217
218    }
219
220    /** {@inheritDoc} */
221    @Override
222    public T[][] getA() {
223
224        final T sqrt6 = getField().getOne().multiply(6).sqrt();
225
226        final T[][] a = MathArrays.buildArray(getField(), 15, -1);
227        for (int i = 0; i < a.length; ++i) {
228            a[i] = MathArrays.buildArray(getField(), i + 1);
229        }
230
231        a[ 0][ 0] = sqrt6.add(-6).divide(-67.5);
232
233        a[ 1][ 0] = sqrt6.add(-6).divide(-180);
234        a[ 1][ 1] = sqrt6.add(-6).divide( -60);
235
236        a[ 2][ 0] = sqrt6.add(-6).divide(-120);
237        a[ 2][ 1] = getField().getZero();
238        a[ 2][ 2] = sqrt6.add(-6).divide( -40);
239
240        a[ 3][ 0] = sqrt6.multiply(107).add(462).divide( 3000);
241        a[ 3][ 1] = getField().getZero();
242        a[ 3][ 2] = sqrt6.multiply(197).add(402).divide(-1000);
243        a[ 3][ 3] = sqrt6.multiply( 73).add(168).divide(  375);
244
245        a[ 4][ 0] = fraction(1, 27);
246        a[ 4][ 1] = getField().getZero();
247        a[ 4][ 2] = getField().getZero();
248        a[ 4][ 3] = sqrt6.add( 16).divide( 108);
249        a[ 4][ 4] = sqrt6.add(-16).divide(-108);
250
251        a[ 5][ 0] = fraction(19, 512);
252        a[ 5][ 1] = getField().getZero();
253        a[ 5][ 2] = getField().getZero();
254        a[ 5][ 3] = sqrt6.multiply( 23).add(118).divide(1024);
255        a[ 5][ 4] = sqrt6.multiply(-23).add(118).divide(1024);
256        a[ 5][ 5] = fraction(-9, 512);
257
258        a[ 6][ 0] = fraction(13772, 371293);
259        a[ 6][ 1] = getField().getZero();
260        a[ 6][ 2] = getField().getZero();
261        a[ 6][ 3] = sqrt6.multiply( 4784).add(51544).divide(371293);
262        a[ 6][ 4] = sqrt6.multiply(-4784).add(51544).divide(371293);
263        a[ 6][ 5] = fraction(-5688, 371293);
264        a[ 6][ 6] = fraction( 3072, 371293);
265
266        a[ 7][ 0] = fraction(58656157643.0, 93983540625.0);
267        a[ 7][ 1] = getField().getZero();
268        a[ 7][ 2] = getField().getZero();
269        a[ 7][ 3] = sqrt6.multiply(-318801444819.0).add(-1324889724104.0).divide(626556937500.0);
270        a[ 7][ 4] = sqrt6.multiply( 318801444819.0).add(-1324889724104.0).divide(626556937500.0);
271        a[ 7][ 5] = fraction(96044563816.0, 3480871875.0);
272        a[ 7][ 6] = fraction(5682451879168.0, 281950621875.0);
273        a[ 7][ 7] = fraction(-165125654.0, 3796875.0);
274
275        a[ 8][ 0] = fraction(8909899.0, 18653125.0);
276        a[ 8][ 1] = getField().getZero();
277        a[ 8][ 2] = getField().getZero();
278        a[ 8][ 3] = sqrt6.multiply(-1137963.0).add(-4521408.0).divide(2937500.0);
279        a[ 8][ 4] = sqrt6.multiply( 1137963.0).add(-4521408.0).divide(2937500.0);
280        a[ 8][ 5] = fraction(96663078.0, 4553125.0);
281        a[ 8][ 6] = fraction(2107245056.0, 137915625.0);
282        a[ 8][ 7] = fraction(-4913652016.0, 147609375.0);
283        a[ 8][ 8] = fraction(-78894270.0, 3880452869.0);
284
285        a[ 9][ 0] = fraction(-20401265806.0, 21769653311.0);
286        a[ 9][ 1] = getField().getZero();
287        a[ 9][ 2] = getField().getZero();
288        a[ 9][ 3] = sqrt6.multiply( 94326.0).add(354216.0).divide(112847.0);
289        a[ 9][ 4] = sqrt6.multiply(-94326.0).add(354216.0).divide(112847.0);
290        a[ 9][ 5] = fraction(-43306765128.0, 5313852383.0);
291        a[ 9][ 6] = fraction(-20866708358144.0, 1126708119789.0);
292        a[ 9][ 7] = fraction(14886003438020.0, 654632330667.0);
293        a[ 9][ 8] = fraction(35290686222309375.0, 14152473387134411.0);
294        a[ 9][ 9] = fraction(-1477884375.0, 485066827.0);
295
296        a[10][ 0] = fraction(39815761.0, 17514443.0);
297        a[10][ 1] = getField().getZero();
298        a[10][ 2] = getField().getZero();
299        a[10][ 3] = sqrt6.multiply(-960905.0).add(-3457480.0).divide(551636.0);
300        a[10][ 4] = sqrt6.multiply( 960905.0).add(-3457480.0).divide(551636.0);
301        a[10][ 5] = fraction(-844554132.0, 47026969.0);
302        a[10][ 6] = fraction(8444996352.0, 302158619.0);
303        a[10][ 7] = fraction(-2509602342.0, 877790785.0);
304        a[10][ 8] = fraction(-28388795297996250.0, 3199510091356783.0);
305        a[10][ 9] = fraction(226716250.0, 18341897.0);
306        a[10][10] = fraction(1371316744.0, 2131383595.0);
307
308        // the following stage is both for interpolation and the first stage in next step
309        // (the coefficients are identical to the B array)
310        a[11][ 0] = fraction(104257.0, 1920240.0);
311        a[11][ 1] = getField().getZero();
312        a[11][ 2] = getField().getZero();
313        a[11][ 3] = getField().getZero();
314        a[11][ 4] = getField().getZero();
315        a[11][ 5] = fraction(3399327.0, 763840.0);
316        a[11][ 6] = fraction(66578432.0, 35198415.0);
317        a[11][ 7] = fraction(-1674902723.0, 288716400.0);
318        a[11][ 8] = fraction(54980371265625.0, 176692375811392.0);
319        a[11][ 9] = fraction(-734375.0, 4826304.0);
320        a[11][10] = fraction(171414593.0, 851261400.0);
321        a[11][11] = fraction(137909.0, 3084480.0);
322
323        // the following stages are for interpolation only
324        a[12][ 0] = fraction(      13481885573.0, 240030000000.0);
325        a[12][ 1] = getField().getZero();
326        a[12][ 2] = getField().getZero();
327        a[12][ 3] = getField().getZero();
328        a[12][ 4] = getField().getZero();
329        a[12][ 5] = getField().getZero();
330        a[12][ 6] = fraction(     139418837528.0, 549975234375.0);
331        a[12][ 7] = fraction(  -11108320068443.0, 45111937500000.0);
332        a[12][ 8] = fraction(-1769651421925959.0, 14249385146080000.0);
333        a[12][ 9] = fraction(         57799439.0, 377055000.0);
334        a[12][10] = fraction(     793322643029.0, 96734250000000.0);
335        a[12][11] = fraction(       1458939311.0, 192780000000.0);
336        a[12][12]  = fraction(            -4149.0, 500000.0);
337
338        a[13][ 0] = fraction(    1595561272731.0, 50120273500000.0);
339        a[13][ 1] = getField().getZero();
340        a[13][ 2] = getField().getZero();
341        a[13][ 3] = getField().getZero();
342        a[13][ 4] = getField().getZero();
343        a[13][ 5] = fraction(     975183916491.0, 34457688031250.0);
344        a[13][ 6] = fraction(   38492013932672.0, 718912673015625.0);
345        a[13][ 7] = fraction(-1114881286517557.0, 20298710767500000.0);
346        a[13][ 8] = getField().getZero();
347        a[13][ 9] = getField().getZero();
348        a[13][10] = fraction(   -2538710946863.0, 23431227861250000.0);
349        a[13][11] = fraction(       8824659001.0, 23066716781250.0);
350        a[13][12] = fraction(     -11518334563.0, 33831184612500.0);
351        a[13][13] = fraction(       1912306948.0, 13532473845.0);
352
353        a[14][ 0] = fraction(     -13613986967.0, 31741908048.0);
354        a[14][ 1] = getField().getZero();
355        a[14][ 2] = getField().getZero();
356        a[14][ 3] = getField().getZero();
357        a[14][ 4] = getField().getZero();
358        a[14][ 5] = fraction(      -4755612631.0, 1012344804.0);
359        a[14][ 6] = fraction(   42939257944576.0, 5588559685701.0);
360        a[14][ 7] = fraction(   77881972900277.0, 19140370552944.0);
361        a[14][ 8] = fraction(   22719829234375.0, 63689648654052.0);
362        a[14][ 9] = getField().getZero();
363        a[14][10] = getField().getZero();
364        a[14][11] = getField().getZero();
365        a[14][12] = fraction(      -1199007803.0, 857031517296.0);
366        a[14][13] = fraction(     157882067000.0, 53564469831.0);
367        a[14][14] = fraction(    -290468882375.0, 31741908048.0);
368
369        return a;
370
371    }
372
373    /** {@inheritDoc} */
374    @Override
375    public T[] getB() {
376        final T[] b = MathArrays.buildArray(getField(), 16);
377        b[ 0] = fraction(104257, 1920240);
378        b[ 1] = getField().getZero();
379        b[ 2] = getField().getZero();
380        b[ 3] = getField().getZero();
381        b[ 4] = getField().getZero();
382        b[ 5] = fraction(        3399327.0,          763840.0);
383        b[ 6] = fraction(       66578432.0,        35198415.0);
384        b[ 7] = fraction(    -1674902723.0,       288716400.0);
385        b[ 8] = fraction( 54980371265625.0, 176692375811392.0);
386        b[ 9] = fraction(        -734375.0,         4826304.0);
387        b[10] = fraction(      171414593.0,       851261400.0);
388        b[11] = fraction(         137909.0,         3084480.0);
389        b[12] = getField().getZero();
390        b[13] = getField().getZero();
391        b[14] = getField().getZero();
392        b[15] = getField().getZero();
393        return b;
394    }
395
396    /** {@inheritDoc} */
397    @Override
398    protected DormandPrince853FieldStepInterpolator<T>
399        createInterpolator(final boolean forward, T[][] yDotK,
400                           final FieldODEStateAndDerivative<T> globalPreviousState,
401                           final FieldODEStateAndDerivative<T> globalCurrentState, final FieldEquationsMapper<T> mapper) {
402        return new DormandPrince853FieldStepInterpolator<>(getField(), forward, yDotK,
403                                                            globalPreviousState, globalCurrentState,
404                                                            globalPreviousState, globalCurrentState,
405                                                            mapper);
406    }
407
408    /** {@inheritDoc} */
409    @Override
410    public int getOrder() {
411        return 8;
412    }
413
414    /** {@inheritDoc} */
415    @Override
416    protected T estimateError(final T[][] yDotK, final T[] y0, final T[] y1, final T h) {
417        T error1 = h.getField().getZero();
418        T error2 = h.getField().getZero();
419
420        for (int j = 0; j < mainSetDimension; ++j) {
421            final T errSum1 =      yDotK[ 0][j].multiply(e1_01).
422                               add(yDotK[ 5][j].multiply(e1_06)).
423                               add(yDotK[ 6][j].multiply(e1_07)).
424                               add(yDotK[ 7][j].multiply(e1_08)).
425                               add(yDotK[ 8][j].multiply(e1_09)).
426                               add(yDotK[ 9][j].multiply(e1_10)).
427                               add(yDotK[10][j].multiply(e1_11)).
428                               add(yDotK[11][j].multiply(e1_12));
429            final T errSum2 =      yDotK[ 0][j].multiply(e2_01).
430                               add(yDotK[ 5][j].multiply(e2_06)).
431                               add(yDotK[ 6][j].multiply(e2_07)).
432                               add(yDotK[ 7][j].multiply(e2_08)).
433                               add(yDotK[ 8][j].multiply(e2_09)).
434                               add(yDotK[ 9][j].multiply(e2_10)).
435                               add(yDotK[10][j].multiply(e2_11)).
436                               add(yDotK[11][j].multiply(e2_12));
437
438            final T yScale = MathUtils.max(y0[j].abs(), y1[j].abs());
439            final T tol = vecAbsoluteTolerance == null ?
440                          yScale.multiply(scalRelativeTolerance).add(scalAbsoluteTolerance) :
441                          yScale.multiply(vecRelativeTolerance[j]).add(vecAbsoluteTolerance[j]);
442            final T ratio1  = errSum1.divide(tol);
443            error1        = error1.add(ratio1.multiply(ratio1));
444            final T ratio2  = errSum2.divide(tol);
445            error2        = error2.add(ratio2.multiply(ratio2));
446        }
447
448        T den = error1.add(error2.multiply(0.01));
449        if (den.getReal() <= 0.0) {
450            den = h.getField().getOne();
451        }
452
453        return h.abs().multiply(error1).divide(den.multiply(mainSetDimension).sqrt());
454
455    }
456
457}