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.math3.ode.nonstiff;
019
020import org.apache.commons.math3.util.FastMath;
021
022
023/**
024 * This class implements the 8(5,3) Dormand-Prince integrator for Ordinary
025 * Differential Equations.
026 *
027 * <p>This integrator is an embedded Runge-Kutta integrator
028 * of order 8(5,3) used in local extrapolation mode (i.e. the solution
029 * is computed using the high order formula) with stepsize control
030 * (and automatic step initialization) and continuous output. This
031 * method uses 12 functions evaluations per step for integration and 4
032 * evaluations for interpolation. However, since the first
033 * interpolation evaluation is the same as the first integration
034 * evaluation of the next step, we have included it in the integrator
035 * rather than in the interpolator and specified the method was an
036 * <i>fsal</i>. Hence, despite we have 13 stages here, the cost is
037 * really 12 evaluations per step even if no interpolation is done,
038 * and the overcost of interpolation is only 3 evaluations.</p>
039 *
040 * <p>This method is based on an 8(6) method by Dormand and Prince
041 * (i.e. order 8 for the integration and order 6 for error estimation)
042 * modified by Hairer and Wanner to use a 5th order error estimator
043 * with 3rd order correction. This modification was introduced because
044 * the original method failed in some cases (wrong steps can be
045 * accepted when step size is too large, for example in the
046 * Brusselator problem) and also had <i>severe difficulties when
047 * applied to problems with discontinuities</i>. This modification is
048 * explained in the second edition of the first volume (Nonstiff
049 * Problems) of the reference book by Hairer, Norsett and Wanner:
050 * <i>Solving Ordinary Differential Equations</i> (Springer-Verlag,
051 * ISBN 3-540-56670-8).</p>
052 *
053 * @version $Id: DormandPrince853Integrator.java 1416643 2012-12-03 19:37:14Z tn $
054 * @since 1.2
055 */
056
057public class DormandPrince853Integrator extends EmbeddedRungeKuttaIntegrator {
058
059  /** Integrator method name. */
060  private static final String METHOD_NAME = "Dormand-Prince 8 (5, 3)";
061
062  /** Time steps Butcher array. */
063  private static final double[] STATIC_C = {
064    (12.0 - 2.0 * FastMath.sqrt(6.0)) / 135.0, (6.0 - FastMath.sqrt(6.0)) / 45.0, (6.0 - FastMath.sqrt(6.0)) / 30.0,
065    (6.0 + FastMath.sqrt(6.0)) / 30.0, 1.0/3.0, 1.0/4.0, 4.0/13.0, 127.0/195.0, 3.0/5.0,
066    6.0/7.0, 1.0, 1.0
067  };
068
069  /** Internal weights Butcher array. */
070  private static final double[][] STATIC_A = {
071
072    // k2
073    {(12.0 - 2.0 * FastMath.sqrt(6.0)) / 135.0},
074
075    // k3
076    {(6.0 - FastMath.sqrt(6.0)) / 180.0, (6.0 - FastMath.sqrt(6.0)) / 60.0},
077
078    // k4
079    {(6.0 - FastMath.sqrt(6.0)) / 120.0, 0.0, (6.0 - FastMath.sqrt(6.0)) / 40.0},
080
081    // k5
082    {(462.0 + 107.0 * FastMath.sqrt(6.0)) / 3000.0, 0.0,
083     (-402.0 - 197.0 * FastMath.sqrt(6.0)) / 1000.0, (168.0 + 73.0 * FastMath.sqrt(6.0)) / 375.0},
084
085    // k6
086    {1.0 / 27.0, 0.0, 0.0, (16.0 + FastMath.sqrt(6.0)) / 108.0, (16.0 - FastMath.sqrt(6.0)) / 108.0},
087
088    // k7
089    {19.0 / 512.0, 0.0, 0.0, (118.0 + 23.0 * FastMath.sqrt(6.0)) / 1024.0,
090     (118.0 - 23.0 * FastMath.sqrt(6.0)) / 1024.0, -9.0 / 512.0},
091
092    // k8
093    {13772.0 / 371293.0, 0.0, 0.0, (51544.0 + 4784.0 * FastMath.sqrt(6.0)) / 371293.0,
094     (51544.0 - 4784.0 * FastMath.sqrt(6.0)) / 371293.0, -5688.0 / 371293.0, 3072.0 / 371293.0},
095
096    // k9
097    {58656157643.0 / 93983540625.0, 0.0, 0.0,
098     (-1324889724104.0 - 318801444819.0 * FastMath.sqrt(6.0)) / 626556937500.0,
099     (-1324889724104.0 + 318801444819.0 * FastMath.sqrt(6.0)) / 626556937500.0,
100     96044563816.0 / 3480871875.0, 5682451879168.0 / 281950621875.0,
101     -165125654.0 / 3796875.0},
102
103    // k10
104    {8909899.0 / 18653125.0, 0.0, 0.0,
105     (-4521408.0 - 1137963.0 * FastMath.sqrt(6.0)) / 2937500.0,
106     (-4521408.0 + 1137963.0 * FastMath.sqrt(6.0)) / 2937500.0,
107     96663078.0 / 4553125.0, 2107245056.0 / 137915625.0,
108     -4913652016.0 / 147609375.0, -78894270.0 / 3880452869.0},
109
110    // k11
111    {-20401265806.0 / 21769653311.0, 0.0, 0.0,
112     (354216.0 + 94326.0 * FastMath.sqrt(6.0)) / 112847.0,
113     (354216.0 - 94326.0 * FastMath.sqrt(6.0)) / 112847.0,
114     -43306765128.0 / 5313852383.0, -20866708358144.0 / 1126708119789.0,
115     14886003438020.0 / 654632330667.0, 35290686222309375.0 / 14152473387134411.0,
116     -1477884375.0 / 485066827.0},
117
118    // k12
119    {39815761.0 / 17514443.0, 0.0, 0.0,
120     (-3457480.0 - 960905.0 * FastMath.sqrt(6.0)) / 551636.0,
121     (-3457480.0 + 960905.0 * FastMath.sqrt(6.0)) / 551636.0,
122     -844554132.0 / 47026969.0, 8444996352.0 / 302158619.0,
123     -2509602342.0 / 877790785.0, -28388795297996250.0 / 3199510091356783.0,
124     226716250.0 / 18341897.0, 1371316744.0 / 2131383595.0},
125
126    // k13 should be for interpolation only, but since it is the same
127    // stage as the first evaluation of the next step, we perform it
128    // here at no cost by specifying this is an fsal method
129    {104257.0/1920240.0, 0.0, 0.0, 0.0, 0.0, 3399327.0/763840.0,
130     66578432.0/35198415.0, -1674902723.0/288716400.0,
131     54980371265625.0/176692375811392.0, -734375.0/4826304.0,
132     171414593.0/851261400.0, 137909.0/3084480.0}
133
134  };
135
136  /** Propagation weights Butcher array. */
137  private static final double[] STATIC_B = {
138      104257.0/1920240.0,
139      0.0,
140      0.0,
141      0.0,
142      0.0,
143      3399327.0/763840.0,
144      66578432.0/35198415.0,
145      -1674902723.0/288716400.0,
146      54980371265625.0/176692375811392.0,
147      -734375.0/4826304.0,
148      171414593.0/851261400.0,
149      137909.0/3084480.0,
150      0.0
151  };
152
153  /** First error weights array, element 1. */
154  private static final double E1_01 =         116092271.0 / 8848465920.0;
155
156  // elements 2 to 5 are zero, so they are neither stored nor used
157
158  /** First error weights array, element 6. */
159  private static final double E1_06 =          -1871647.0 / 1527680.0;
160
161  /** First error weights array, element 7. */
162  private static final double E1_07 =         -69799717.0 / 140793660.0;
163
164  /** First error weights array, element 8. */
165  private static final double E1_08 =     1230164450203.0 / 739113984000.0;
166
167  /** First error weights array, element 9. */
168  private static final double E1_09 = -1980813971228885.0 / 5654156025964544.0;
169
170  /** First error weights array, element 10. */
171  private static final double E1_10 =         464500805.0 / 1389975552.0;
172
173  /** First error weights array, element 11. */
174  private static final double E1_11 =     1606764981773.0 / 19613062656000.0;
175
176  /** First error weights array, element 12. */
177  private static final double E1_12 =           -137909.0 / 6168960.0;
178
179
180  /** Second error weights array, element 1. */
181  private static final double E2_01 =           -364463.0 / 1920240.0;
182
183  // elements 2 to 5 are zero, so they are neither stored nor used
184
185  /** Second error weights array, element 6. */
186  private static final double E2_06 =           3399327.0 / 763840.0;
187
188  /** Second error weights array, element 7. */
189  private static final double E2_07 =          66578432.0 / 35198415.0;
190
191  /** Second error weights array, element 8. */
192  private static final double E2_08 =       -1674902723.0 / 288716400.0;
193
194  /** Second error weights array, element 9. */
195  private static final double E2_09 =   -74684743568175.0 / 176692375811392.0;
196
197  /** Second error weights array, element 10. */
198  private static final double E2_10 =           -734375.0 / 4826304.0;
199
200  /** Second error weights array, element 11. */
201  private static final double E2_11 =         171414593.0 / 851261400.0;
202
203  /** Second error weights array, element 12. */
204  private static final double E2_12 =             69869.0 / 3084480.0;
205
206  /** Simple constructor.
207   * Build an eighth order Dormand-Prince integrator with the given step bounds
208   * @param minStep minimal step (sign is irrelevant, regardless of
209   * integration direction, forward or backward), the last step can
210   * be smaller than this
211   * @param maxStep maximal step (sign is irrelevant, regardless of
212   * integration direction, forward or backward), the last step can
213   * be smaller than this
214   * @param scalAbsoluteTolerance allowed absolute error
215   * @param scalRelativeTolerance allowed relative error
216   */
217  public DormandPrince853Integrator(final double minStep, final double maxStep,
218                                    final double scalAbsoluteTolerance,
219                                    final double scalRelativeTolerance) {
220    super(METHOD_NAME, true, STATIC_C, STATIC_A, STATIC_B,
221          new DormandPrince853StepInterpolator(),
222          minStep, maxStep, scalAbsoluteTolerance, scalRelativeTolerance);
223  }
224
225  /** Simple constructor.
226   * Build an eighth order Dormand-Prince integrator with the given step bounds
227   * @param minStep minimal step (sign is irrelevant, regardless of
228   * integration direction, forward or backward), the last step can
229   * be smaller than this
230   * @param maxStep maximal step (sign is irrelevant, regardless of
231   * integration direction, forward or backward), the last step can
232   * be smaller than this
233   * @param vecAbsoluteTolerance allowed absolute error
234   * @param vecRelativeTolerance allowed relative error
235   */
236  public DormandPrince853Integrator(final double minStep, final double maxStep,
237                                    final double[] vecAbsoluteTolerance,
238                                    final double[] vecRelativeTolerance) {
239    super(METHOD_NAME, true, STATIC_C, STATIC_A, STATIC_B,
240          new DormandPrince853StepInterpolator(),
241          minStep, maxStep, vecAbsoluteTolerance, vecRelativeTolerance);
242  }
243
244  /** {@inheritDoc} */
245  @Override
246  public int getOrder() {
247    return 8;
248  }
249
250  /** {@inheritDoc} */
251  @Override
252  protected double estimateError(final double[][] yDotK,
253                                 final double[] y0, final double[] y1,
254                                 final double h) {
255    double error1 = 0;
256    double error2 = 0;
257
258    for (int j = 0; j < mainSetDimension; ++j) {
259      final double errSum1 = E1_01 * yDotK[0][j]  + E1_06 * yDotK[5][j] +
260                             E1_07 * yDotK[6][j]  + E1_08 * yDotK[7][j] +
261                             E1_09 * yDotK[8][j]  + E1_10 * yDotK[9][j] +
262                             E1_11 * yDotK[10][j] + E1_12 * yDotK[11][j];
263      final double errSum2 = E2_01 * yDotK[0][j]  + E2_06 * yDotK[5][j] +
264                             E2_07 * yDotK[6][j]  + E2_08 * yDotK[7][j] +
265                             E2_09 * yDotK[8][j]  + E2_10 * yDotK[9][j] +
266                             E2_11 * yDotK[10][j] + E2_12 * yDotK[11][j];
267
268      final double yScale = FastMath.max(FastMath.abs(y0[j]), FastMath.abs(y1[j]));
269      final double tol = (vecAbsoluteTolerance == null) ?
270                         (scalAbsoluteTolerance + scalRelativeTolerance * yScale) :
271                         (vecAbsoluteTolerance[j] + vecRelativeTolerance[j] * yScale);
272      final double ratio1  = errSum1 / tol;
273      error1        += ratio1 * ratio1;
274      final double ratio2  = errSum2 / tol;
275      error2        += ratio2 * ratio2;
276    }
277
278    double den = error1 + 0.01 * error2;
279    if (den <= 0.0) {
280      den = 1.0;
281    }
282
283    return FastMath.abs(h) * error1 / FastMath.sqrt(mainSetDimension * den);
284
285  }
286
287}