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.legacy.ode.nonstiff; 019 020import org.apache.commons.math4.core.jdkmath.JdkMath; 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 * @since 1.2 054 */ 055 056public class DormandPrince853Integrator extends EmbeddedRungeKuttaIntegrator { 057 058 /** Integrator method name. */ 059 private static final String METHOD_NAME = "Dormand-Prince 8 (5, 3)"; 060 061 /** Time steps Butcher array. */ 062 private static final double[] STATIC_C = { 063 (12.0 - 2.0 * JdkMath.sqrt(6.0)) / 135.0, (6.0 - JdkMath.sqrt(6.0)) / 45.0, (6.0 - JdkMath.sqrt(6.0)) / 30.0, 064 (6.0 + JdkMath.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, 065 6.0/7.0, 1.0, 1.0 066 }; 067 068 /** Internal weights Butcher array. */ 069 private static final double[][] STATIC_A = { 070 071 // k2 072 {(12.0 - 2.0 * JdkMath.sqrt(6.0)) / 135.0}, 073 074 // k3 075 {(6.0 - JdkMath.sqrt(6.0)) / 180.0, (6.0 - JdkMath.sqrt(6.0)) / 60.0}, 076 077 // k4 078 {(6.0 - JdkMath.sqrt(6.0)) / 120.0, 0.0, (6.0 - JdkMath.sqrt(6.0)) / 40.0}, 079 080 // k5 081 {(462.0 + 107.0 * JdkMath.sqrt(6.0)) / 3000.0, 0.0, 082 (-402.0 - 197.0 * JdkMath.sqrt(6.0)) / 1000.0, (168.0 + 73.0 * JdkMath.sqrt(6.0)) / 375.0}, 083 084 // k6 085 {1.0 / 27.0, 0.0, 0.0, (16.0 + JdkMath.sqrt(6.0)) / 108.0, (16.0 - JdkMath.sqrt(6.0)) / 108.0}, 086 087 // k7 088 {19.0 / 512.0, 0.0, 0.0, (118.0 + 23.0 * JdkMath.sqrt(6.0)) / 1024.0, 089 (118.0 - 23.0 * JdkMath.sqrt(6.0)) / 1024.0, -9.0 / 512.0}, 090 091 // k8 092 {13772.0 / 371293.0, 0.0, 0.0, (51544.0 + 4784.0 * JdkMath.sqrt(6.0)) / 371293.0, 093 (51544.0 - 4784.0 * JdkMath.sqrt(6.0)) / 371293.0, -5688.0 / 371293.0, 3072.0 / 371293.0}, 094 095 // k9 096 {58656157643.0 / 93983540625.0, 0.0, 0.0, 097 (-1324889724104.0 - 318801444819.0 * JdkMath.sqrt(6.0)) / 626556937500.0, 098 (-1324889724104.0 + 318801444819.0 * JdkMath.sqrt(6.0)) / 626556937500.0, 099 96044563816.0 / 3480871875.0, 5682451879168.0 / 281950621875.0, 100 -165125654.0 / 3796875.0}, 101 102 // k10 103 {8909899.0 / 18653125.0, 0.0, 0.0, 104 (-4521408.0 - 1137963.0 * JdkMath.sqrt(6.0)) / 2937500.0, 105 (-4521408.0 + 1137963.0 * JdkMath.sqrt(6.0)) / 2937500.0, 106 96663078.0 / 4553125.0, 2107245056.0 / 137915625.0, 107 -4913652016.0 / 147609375.0, -78894270.0 / 3880452869.0}, 108 109 // k11 110 {-20401265806.0 / 21769653311.0, 0.0, 0.0, 111 (354216.0 + 94326.0 * JdkMath.sqrt(6.0)) / 112847.0, 112 (354216.0 - 94326.0 * JdkMath.sqrt(6.0)) / 112847.0, 113 -43306765128.0 / 5313852383.0, -20866708358144.0 / 1126708119789.0, 114 14886003438020.0 / 654632330667.0, 35290686222309375.0 / 14152473387134411.0, 115 -1477884375.0 / 485066827.0}, 116 117 // k12 118 {39815761.0 / 17514443.0, 0.0, 0.0, 119 (-3457480.0 - 960905.0 * JdkMath.sqrt(6.0)) / 551636.0, 120 (-3457480.0 + 960905.0 * JdkMath.sqrt(6.0)) / 551636.0, 121 -844554132.0 / 47026969.0, 8444996352.0 / 302158619.0, 122 -2509602342.0 / 877790785.0, -28388795297996250.0 / 3199510091356783.0, 123 226716250.0 / 18341897.0, 1371316744.0 / 2131383595.0}, 124 125 // k13 should be for interpolation only, but since it is the same 126 // stage as the first evaluation of the next step, we perform it 127 // here at no cost by specifying this is an fsal method 128 {104257.0/1920240.0, 0.0, 0.0, 0.0, 0.0, 3399327.0/763840.0, 129 66578432.0/35198415.0, -1674902723.0/288716400.0, 130 54980371265625.0/176692375811392.0, -734375.0/4826304.0, 131 171414593.0/851261400.0, 137909.0/3084480.0} 132 }; 133 134 /** Propagation weights Butcher array. */ 135 private static final double[] STATIC_B = { 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 (sign is irrelevant, regardless of 207 * integration direction, forward or backward), the last step can 208 * be smaller than this 209 * @param maxStep maximal step (sign is irrelevant, regardless of 210 * integration direction, forward or backward), the last step can 211 * be smaller than this 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 (sign is irrelevant, regardless of 226 * integration direction, forward or backward), the last step can 227 * be smaller than this 228 * @param maxStep maximal step (sign is irrelevant, regardless of 229 * integration direction, forward or backward), the last step can 230 * be smaller than this 231 * @param vecAbsoluteTolerance allowed absolute error 232 * @param vecRelativeTolerance allowed relative error 233 */ 234 public DormandPrince853Integrator(final double minStep, final double maxStep, 235 final double[] vecAbsoluteTolerance, 236 final double[] vecRelativeTolerance) { 237 super(METHOD_NAME, true, STATIC_C, STATIC_A, STATIC_B, 238 new DormandPrince853StepInterpolator(), 239 minStep, maxStep, vecAbsoluteTolerance, vecRelativeTolerance); 240 } 241 242 /** {@inheritDoc} */ 243 @Override 244 public int getOrder() { 245 return 8; 246 } 247 248 /** {@inheritDoc} */ 249 @Override 250 protected double estimateError(final double[][] yDotK, 251 final double[] y0, final double[] y1, 252 final double h) { 253 double error1 = 0; 254 double error2 = 0; 255 256 for (int j = 0; j < mainSetDimension; ++j) { 257 final double errSum1 = E1_01 * yDotK[0][j] + E1_06 * yDotK[5][j] + 258 E1_07 * yDotK[6][j] + E1_08 * yDotK[7][j] + 259 E1_09 * yDotK[8][j] + E1_10 * yDotK[9][j] + 260 E1_11 * yDotK[10][j] + E1_12 * yDotK[11][j]; 261 final double errSum2 = E2_01 * yDotK[0][j] + E2_06 * yDotK[5][j] + 262 E2_07 * yDotK[6][j] + E2_08 * yDotK[7][j] + 263 E2_09 * yDotK[8][j] + E2_10 * yDotK[9][j] + 264 E2_11 * yDotK[10][j] + E2_12 * yDotK[11][j]; 265 266 final double yScale = JdkMath.max(JdkMath.abs(y0[j]), JdkMath.abs(y1[j])); 267 final double tol = (vecAbsoluteTolerance == null) ? 268 (scalAbsoluteTolerance + scalRelativeTolerance * yScale) : 269 (vecAbsoluteTolerance[j] + vecRelativeTolerance[j] * yScale); 270 final double ratio1 = errSum1 / tol; 271 error1 += ratio1 * ratio1; 272 final double ratio2 = errSum2 / tol; 273 error2 += ratio2 * ratio2; 274 } 275 276 double den = error1 + 0.01 * error2; 277 if (den <= 0.0) { 278 den = 1.0; 279 } 280 281 return JdkMath.abs(h) * error1 / JdkMath.sqrt(mainSetDimension * den); 282 } 283}