Home | History | Annotate | Download | only in nonstiff
      1 /*
      2  * Licensed to the Apache Software Foundation (ASF) 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 ASF 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  *      http://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.apache.commons.math.ode.nonstiff;
     19 
     20 import org.apache.commons.math.util.FastMath;
     21 
     22 
     23 /**
     24  * This class implements the 8(5,3) Dormand-Prince integrator for Ordinary
     25  * Differential Equations.
     26  *
     27  * <p>This integrator is an embedded Runge-Kutta integrator
     28  * of order 8(5,3) used in local extrapolation mode (i.e. the solution
     29  * is computed using the high order formula) with stepsize control
     30  * (and automatic step initialization) and continuous output. This
     31  * method uses 12 functions evaluations per step for integration and 4
     32  * evaluations for interpolation. However, since the first
     33  * interpolation evaluation is the same as the first integration
     34  * evaluation of the next step, we have included it in the integrator
     35  * rather than in the interpolator and specified the method was an
     36  * <i>fsal</i>. Hence, despite we have 13 stages here, the cost is
     37  * really 12 evaluations per step even if no interpolation is done,
     38  * and the overcost of interpolation is only 3 evaluations.</p>
     39  *
     40  * <p>This method is based on an 8(6) method by Dormand and Prince
     41  * (i.e. order 8 for the integration and order 6 for error estimation)
     42  * modified by Hairer and Wanner to use a 5th order error estimator
     43  * with 3rd order correction. This modification was introduced because
     44  * the original method failed in some cases (wrong steps can be
     45  * accepted when step size is too large, for example in the
     46  * Brusselator problem) and also had <i>severe difficulties when
     47  * applied to problems with discontinuities</i>. This modification is
     48  * explained in the second edition of the first volume (Nonstiff
     49  * Problems) of the reference book by Hairer, Norsett and Wanner:
     50  * <i>Solving Ordinary Differential Equations</i> (Springer-Verlag,
     51  * ISBN 3-540-56670-8).</p>
     52  *
     53  * @version $Revision: 990655 $ $Date: 2010-08-29 23:49:40 +0200 (dim. 29 aot 2010) $
     54  * @since 1.2
     55  */
     56 
     57 public class DormandPrince853Integrator extends EmbeddedRungeKuttaIntegrator {
     58 
     59   /** Integrator method name. */
     60   private static final String METHOD_NAME = "Dormand-Prince 8 (5, 3)";
     61 
     62   /** Time steps Butcher array. */
     63   private static final double[] STATIC_C = {
     64     (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,
     65     (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,
     66     6.0/7.0, 1.0, 1.0
     67   };
     68 
     69   /** Internal weights Butcher array. */
     70   private static final double[][] STATIC_A = {
     71 
     72     // k2
     73     {(12.0 - 2.0 * FastMath.sqrt(6.0)) / 135.0},
     74 
     75     // k3
     76     {(6.0 - FastMath.sqrt(6.0)) / 180.0, (6.0 - FastMath.sqrt(6.0)) / 60.0},
     77 
     78     // k4
     79     {(6.0 - FastMath.sqrt(6.0)) / 120.0, 0.0, (6.0 - FastMath.sqrt(6.0)) / 40.0},
     80 
     81     // k5
     82     {(462.0 + 107.0 * FastMath.sqrt(6.0)) / 3000.0, 0.0,
     83      (-402.0 - 197.0 * FastMath.sqrt(6.0)) / 1000.0, (168.0 + 73.0 * FastMath.sqrt(6.0)) / 375.0},
     84 
     85     // k6
     86     {1.0 / 27.0, 0.0, 0.0, (16.0 + FastMath.sqrt(6.0)) / 108.0, (16.0 - FastMath.sqrt(6.0)) / 108.0},
     87 
     88     // k7
     89     {19.0 / 512.0, 0.0, 0.0, (118.0 + 23.0 * FastMath.sqrt(6.0)) / 1024.0,
     90      (118.0 - 23.0 * FastMath.sqrt(6.0)) / 1024.0, -9.0 / 512.0},
     91 
     92     // k8
     93     {13772.0 / 371293.0, 0.0, 0.0, (51544.0 + 4784.0 * FastMath.sqrt(6.0)) / 371293.0,
     94      (51544.0 - 4784.0 * FastMath.sqrt(6.0)) / 371293.0, -5688.0 / 371293.0, 3072.0 / 371293.0},
     95 
     96     // k9
     97     {58656157643.0 / 93983540625.0, 0.0, 0.0,
     98      (-1324889724104.0 - 318801444819.0 * FastMath.sqrt(6.0)) / 626556937500.0,
     99      (-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 (must be positive even for backward
    209    * integration), the last step can be smaller than this
    210    * @param maxStep maximal step (must be positive even for backward
    211    * integration)
    212    * @param scalAbsoluteTolerance allowed absolute error
    213    * @param scalRelativeTolerance allowed relative error
    214    */
    215   public DormandPrince853Integrator(final double minStep, final double maxStep,
    216                                     final double scalAbsoluteTolerance,
    217                                     final double scalRelativeTolerance) {
    218     super(METHOD_NAME, true, STATIC_C, STATIC_A, STATIC_B,
    219           new DormandPrince853StepInterpolator(),
    220           minStep, maxStep, scalAbsoluteTolerance, scalRelativeTolerance);
    221   }
    222 
    223   /** Simple constructor.
    224    * Build an eighth order Dormand-Prince integrator with the given step bounds
    225    * @param minStep minimal step (must be positive even for backward
    226    * integration), the last step can be smaller than this
    227    * @param maxStep maximal step (must be positive even for backward
    228    * integration)
    229    * @param vecAbsoluteTolerance allowed absolute error
    230    * @param vecRelativeTolerance allowed relative error
    231    */
    232   public DormandPrince853Integrator(final double minStep, final double maxStep,
    233                                     final double[] vecAbsoluteTolerance,
    234                                     final double[] vecRelativeTolerance) {
    235     super(METHOD_NAME, true, STATIC_C, STATIC_A, STATIC_B,
    236           new DormandPrince853StepInterpolator(),
    237           minStep, maxStep, vecAbsoluteTolerance, vecRelativeTolerance);
    238   }
    239 
    240   /** {@inheritDoc} */
    241   @Override
    242   public int getOrder() {
    243     return 8;
    244   }
    245 
    246   /** {@inheritDoc} */
    247   @Override
    248   protected double estimateError(final double[][] yDotK,
    249                                  final double[] y0, final double[] y1,
    250                                  final double h) {
    251     double error1 = 0;
    252     double error2 = 0;
    253 
    254     for (int j = 0; j < mainSetDimension; ++j) {
    255       final double errSum1 = E1_01 * yDotK[0][j]  + E1_06 * yDotK[5][j] +
    256                              E1_07 * yDotK[6][j]  + E1_08 * yDotK[7][j] +
    257                              E1_09 * yDotK[8][j]  + E1_10 * yDotK[9][j] +
    258                              E1_11 * yDotK[10][j] + E1_12 * yDotK[11][j];
    259       final double errSum2 = E2_01 * yDotK[0][j]  + E2_06 * yDotK[5][j] +
    260                              E2_07 * yDotK[6][j]  + E2_08 * yDotK[7][j] +
    261                              E2_09 * yDotK[8][j]  + E2_10 * yDotK[9][j] +
    262                              E2_11 * yDotK[10][j] + E2_12 * yDotK[11][j];
    263 
    264       final double yScale = FastMath.max(FastMath.abs(y0[j]), FastMath.abs(y1[j]));
    265       final double tol = (vecAbsoluteTolerance == null) ?
    266                          (scalAbsoluteTolerance + scalRelativeTolerance * yScale) :
    267                          (vecAbsoluteTolerance[j] + vecRelativeTolerance[j] * yScale);
    268       final double ratio1  = errSum1 / tol;
    269       error1        += ratio1 * ratio1;
    270       final double ratio2  = errSum2 / tol;
    271       error2        += ratio2 * ratio2;
    272     }
    273 
    274     double den = error1 + 0.01 * error2;
    275     if (den <= 0.0) {
    276       den = 1.0;
    277     }
    278 
    279     return FastMath.abs(h) * error1 / FastMath.sqrt(mainSetDimension * den);
    280 
    281   }
    282 
    283 }
    284