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