Raymond | dee0849 | 2015-04-02 10:43:13 -0700 | [diff] [blame] | 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 août 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 | } |