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.legacy.core.Field; 021import org.apache.commons.math4.legacy.core.RealFieldElement; 022import org.apache.commons.math4.legacy.ode.FieldEquationsMapper; 023import org.apache.commons.math4.legacy.ode.FieldODEStateAndDerivative; 024import org.apache.commons.math4.legacy.core.MathArrays; 025 026 027/** 028 * This class implements the 8(5,3) Dormand-Prince integrator for Ordinary 029 * Differential Equations. 030 * 031 * <p>This integrator is an embedded Runge-Kutta integrator 032 * of order 8(5,3) used in local extrapolation mode (i.e. the solution 033 * is computed using the high order formula) with stepsize control 034 * (and automatic step initialization) and continuous output. This 035 * method uses 12 functions evaluations per step for integration and 4 036 * evaluations for interpolation. However, since the first 037 * interpolation evaluation is the same as the first integration 038 * evaluation of the next step, we have included it in the integrator 039 * rather than in the interpolator and specified the method was an 040 * <i>fsal</i>. Hence, despite we have 13 stages here, the cost is 041 * really 12 evaluations per step even if no interpolation is done, 042 * and the overcost of interpolation is only 3 evaluations.</p> 043 * 044 * <p>This method is based on an 8(6) method by Dormand and Prince 045 * (i.e. order 8 for the integration and order 6 for error estimation) 046 * modified by Hairer and Wanner to use a 5th order error estimator 047 * with 3rd order correction. This modification was introduced because 048 * the original method failed in some cases (wrong steps can be 049 * accepted when step size is too large, for example in the 050 * Brusselator problem) and also had <i>severe difficulties when 051 * applied to problems with discontinuities</i>. This modification is 052 * explained in the second edition of the first volume (Nonstiff 053 * Problems) of the reference book by Hairer, Norsett and Wanner: 054 * <i>Solving Ordinary Differential Equations</i> (Springer-Verlag, 055 * ISBN 3-540-56670-8).</p> 056 * 057 * @param <T> the type of the field elements 058 * @since 3.6 059 */ 060 061public class DormandPrince853FieldIntegrator<T extends RealFieldElement<T>> 062 extends EmbeddedRungeKuttaFieldIntegrator<T> { 063 064 /** Integrator method name. */ 065 private static final String METHOD_NAME = "Dormand-Prince 8 (5, 3)"; 066 067 /** First error weights array, element 1. */ 068 private final T e1_01; 069 070 // elements 2 to 5 are zero, so they are neither stored nor used 071 072 /** First error weights array, element 6. */ 073 private final T e1_06; 074 075 /** First error weights array, element 7. */ 076 private final T e1_07; 077 078 /** First error weights array, element 8. */ 079 private final T e1_08; 080 081 /** First error weights array, element 9. */ 082 private final T e1_09; 083 084 /** First error weights array, element 10. */ 085 private final T e1_10; 086 087 /** First error weights array, element 11. */ 088 private final T e1_11; 089 090 /** First error weights array, element 12. */ 091 private final T e1_12; 092 093 094 /** Second error weights array, element 1. */ 095 private final T e2_01; 096 097 // elements 2 to 5 are zero, so they are neither stored nor used 098 099 /** Second error weights array, element 6. */ 100 private final T e2_06; 101 102 /** Second error weights array, element 7. */ 103 private final T e2_07; 104 105 /** Second error weights array, element 8. */ 106 private final T e2_08; 107 108 /** Second error weights array, element 9. */ 109 private final T e2_09; 110 111 /** Second error weights array, element 10. */ 112 private final T e2_10; 113 114 /** Second error weights array, element 11. */ 115 private final T e2_11; 116 117 /** Second error weights array, element 12. */ 118 private final T e2_12; 119 120 /** Simple constructor. 121 * Build an eighth order Dormand-Prince integrator with the given step bounds 122 * @param field field to which the time and state vector elements belong 123 * @param minStep minimal step (sign is irrelevant, regardless of 124 * integration direction, forward or backward), the last step can 125 * be smaller than this 126 * @param maxStep maximal step (sign is irrelevant, regardless of 127 * integration direction, forward or backward), the last step can 128 * be smaller than this 129 * @param scalAbsoluteTolerance allowed absolute error 130 * @param scalRelativeTolerance allowed relative error 131 */ 132 public DormandPrince853FieldIntegrator(final Field<T> field, 133 final double minStep, final double maxStep, 134 final double scalAbsoluteTolerance, 135 final double scalRelativeTolerance) { 136 super(field, METHOD_NAME, 12, 137 minStep, maxStep, scalAbsoluteTolerance, scalRelativeTolerance); 138 e1_01 = fraction( 116092271.0, 8848465920.0); 139 e1_06 = fraction( -1871647.0, 1527680.0); 140 e1_07 = fraction( -69799717.0, 140793660.0); 141 e1_08 = fraction( 1230164450203.0, 739113984000.0); 142 e1_09 = fraction(-1980813971228885.0, 5654156025964544.0); 143 e1_10 = fraction( 464500805.0, 1389975552.0); 144 e1_11 = fraction( 1606764981773.0, 19613062656000.0); 145 e1_12 = fraction( -137909.0, 6168960.0); 146 e2_01 = fraction( -364463.0, 1920240.0); 147 e2_06 = fraction( 3399327.0, 763840.0); 148 e2_07 = fraction( 66578432.0, 35198415.0); 149 e2_08 = fraction( -1674902723.0, 288716400.0); 150 e2_09 = fraction( -74684743568175.0, 176692375811392.0); 151 e2_10 = fraction( -734375.0, 4826304.0); 152 e2_11 = fraction( 171414593.0, 851261400.0); 153 e2_12 = fraction( 69869.0, 3084480.0); 154 } 155 156 /** Simple constructor. 157 * Build an eighth order Dormand-Prince integrator with the given step bounds 158 * @param field field to which the time and state vector elements belong 159 * @param minStep minimal step (sign is irrelevant, regardless of 160 * integration direction, forward or backward), the last step can 161 * be smaller than this 162 * @param maxStep maximal step (sign is irrelevant, regardless of 163 * integration direction, forward or backward), the last step can 164 * be smaller than this 165 * @param vecAbsoluteTolerance allowed absolute error 166 * @param vecRelativeTolerance allowed relative error 167 */ 168 public DormandPrince853FieldIntegrator(final Field<T> field, 169 final double minStep, final double maxStep, 170 final double[] vecAbsoluteTolerance, 171 final double[] vecRelativeTolerance) { 172 super(field, METHOD_NAME, 12, 173 minStep, maxStep, vecAbsoluteTolerance, vecRelativeTolerance); 174 e1_01 = fraction( 116092271.0, 8848465920.0); 175 e1_06 = fraction( -1871647.0, 1527680.0); 176 e1_07 = fraction( -69799717.0, 140793660.0); 177 e1_08 = fraction( 1230164450203.0, 739113984000.0); 178 e1_09 = fraction(-1980813971228885.0, 5654156025964544.0); 179 e1_10 = fraction( 464500805.0, 1389975552.0); 180 e1_11 = fraction( 1606764981773.0, 19613062656000.0); 181 e1_12 = fraction( -137909.0, 6168960.0); 182 e2_01 = fraction( -364463.0, 1920240.0); 183 e2_06 = fraction( 3399327.0, 763840.0); 184 e2_07 = fraction( 66578432.0, 35198415.0); 185 e2_08 = fraction( -1674902723.0, 288716400.0); 186 e2_09 = fraction( -74684743568175.0, 176692375811392.0); 187 e2_10 = fraction( -734375.0, 4826304.0); 188 e2_11 = fraction( 171414593.0, 851261400.0); 189 e2_12 = fraction( 69869.0, 3084480.0); 190 } 191 192 /** {@inheritDoc} */ 193 @Override 194 public T[] getC() { 195 196 final T sqrt6 = getField().getOne().multiply(6).sqrt(); 197 198 final T[] c = MathArrays.buildArray(getField(), 15); 199 c[ 0] = sqrt6.add(-6).divide(-67.5); 200 c[ 1] = sqrt6.add(-6).divide(-45.0); 201 c[ 2] = sqrt6.add(-6).divide(-30.0); 202 c[ 3] = sqrt6.add( 6).divide( 30.0); 203 c[ 4] = fraction(1, 3); 204 c[ 5] = fraction(1, 4); 205 c[ 6] = fraction(4, 13); 206 c[ 7] = fraction(127, 195); 207 c[ 8] = fraction(3, 5); 208 c[ 9] = fraction(6, 7); 209 c[10] = getField().getOne(); 210 c[11] = getField().getOne(); 211 c[12] = fraction(1.0, 10.0); 212 c[13] = fraction(1.0, 5.0); 213 c[14] = fraction(7.0, 9.0); 214 215 return c; 216 } 217 218 /** {@inheritDoc} */ 219 @Override 220 public T[][] getA() { 221 222 final T sqrt6 = getField().getOne().multiply(6).sqrt(); 223 224 final T[][] a = MathArrays.buildArray(getField(), 15, -1); 225 for (int i = 0; i < a.length; ++i) { 226 a[i] = MathArrays.buildArray(getField(), i + 1); 227 } 228 229 a[ 0][ 0] = sqrt6.add(-6).divide(-67.5); 230 231 a[ 1][ 0] = sqrt6.add(-6).divide(-180); 232 a[ 1][ 1] = sqrt6.add(-6).divide( -60); 233 234 a[ 2][ 0] = sqrt6.add(-6).divide(-120); 235 a[ 2][ 1] = getField().getZero(); 236 a[ 2][ 2] = sqrt6.add(-6).divide( -40); 237 238 a[ 3][ 0] = sqrt6.multiply(107).add(462).divide( 3000); 239 a[ 3][ 1] = getField().getZero(); 240 a[ 3][ 2] = sqrt6.multiply(197).add(402).divide(-1000); 241 a[ 3][ 3] = sqrt6.multiply( 73).add(168).divide( 375); 242 243 a[ 4][ 0] = fraction(1, 27); 244 a[ 4][ 1] = getField().getZero(); 245 a[ 4][ 2] = getField().getZero(); 246 a[ 4][ 3] = sqrt6.add( 16).divide( 108); 247 a[ 4][ 4] = sqrt6.add(-16).divide(-108); 248 249 a[ 5][ 0] = fraction(19, 512); 250 a[ 5][ 1] = getField().getZero(); 251 a[ 5][ 2] = getField().getZero(); 252 a[ 5][ 3] = sqrt6.multiply( 23).add(118).divide(1024); 253 a[ 5][ 4] = sqrt6.multiply(-23).add(118).divide(1024); 254 a[ 5][ 5] = fraction(-9, 512); 255 256 a[ 6][ 0] = fraction(13772, 371293); 257 a[ 6][ 1] = getField().getZero(); 258 a[ 6][ 2] = getField().getZero(); 259 a[ 6][ 3] = sqrt6.multiply( 4784).add(51544).divide(371293); 260 a[ 6][ 4] = sqrt6.multiply(-4784).add(51544).divide(371293); 261 a[ 6][ 5] = fraction(-5688, 371293); 262 a[ 6][ 6] = fraction( 3072, 371293); 263 264 a[ 7][ 0] = fraction(58656157643.0, 93983540625.0); 265 a[ 7][ 1] = getField().getZero(); 266 a[ 7][ 2] = getField().getZero(); 267 a[ 7][ 3] = sqrt6.multiply(-318801444819.0).add(-1324889724104.0).divide(626556937500.0); 268 a[ 7][ 4] = sqrt6.multiply( 318801444819.0).add(-1324889724104.0).divide(626556937500.0); 269 a[ 7][ 5] = fraction(96044563816.0, 3480871875.0); 270 a[ 7][ 6] = fraction(5682451879168.0, 281950621875.0); 271 a[ 7][ 7] = fraction(-165125654.0, 3796875.0); 272 273 a[ 8][ 0] = fraction(8909899.0, 18653125.0); 274 a[ 8][ 1] = getField().getZero(); 275 a[ 8][ 2] = getField().getZero(); 276 a[ 8][ 3] = sqrt6.multiply(-1137963.0).add(-4521408.0).divide(2937500.0); 277 a[ 8][ 4] = sqrt6.multiply( 1137963.0).add(-4521408.0).divide(2937500.0); 278 a[ 8][ 5] = fraction(96663078.0, 4553125.0); 279 a[ 8][ 6] = fraction(2107245056.0, 137915625.0); 280 a[ 8][ 7] = fraction(-4913652016.0, 147609375.0); 281 a[ 8][ 8] = fraction(-78894270.0, 3880452869.0); 282 283 a[ 9][ 0] = fraction(-20401265806.0, 21769653311.0); 284 a[ 9][ 1] = getField().getZero(); 285 a[ 9][ 2] = getField().getZero(); 286 a[ 9][ 3] = sqrt6.multiply( 94326.0).add(354216.0).divide(112847.0); 287 a[ 9][ 4] = sqrt6.multiply(-94326.0).add(354216.0).divide(112847.0); 288 a[ 9][ 5] = fraction(-43306765128.0, 5313852383.0); 289 a[ 9][ 6] = fraction(-20866708358144.0, 1126708119789.0); 290 a[ 9][ 7] = fraction(14886003438020.0, 654632330667.0); 291 a[ 9][ 8] = fraction(35290686222309375.0, 14152473387134411.0); 292 a[ 9][ 9] = fraction(-1477884375.0, 485066827.0); 293 294 a[10][ 0] = fraction(39815761.0, 17514443.0); 295 a[10][ 1] = getField().getZero(); 296 a[10][ 2] = getField().getZero(); 297 a[10][ 3] = sqrt6.multiply(-960905.0).add(-3457480.0).divide(551636.0); 298 a[10][ 4] = sqrt6.multiply( 960905.0).add(-3457480.0).divide(551636.0); 299 a[10][ 5] = fraction(-844554132.0, 47026969.0); 300 a[10][ 6] = fraction(8444996352.0, 302158619.0); 301 a[10][ 7] = fraction(-2509602342.0, 877790785.0); 302 a[10][ 8] = fraction(-28388795297996250.0, 3199510091356783.0); 303 a[10][ 9] = fraction(226716250.0, 18341897.0); 304 a[10][10] = fraction(1371316744.0, 2131383595.0); 305 306 // the following stage is both for interpolation and the first stage in next step 307 // (the coefficients are identical to the B array) 308 a[11][ 0] = fraction(104257.0, 1920240.0); 309 a[11][ 1] = getField().getZero(); 310 a[11][ 2] = getField().getZero(); 311 a[11][ 3] = getField().getZero(); 312 a[11][ 4] = getField().getZero(); 313 a[11][ 5] = fraction(3399327.0, 763840.0); 314 a[11][ 6] = fraction(66578432.0, 35198415.0); 315 a[11][ 7] = fraction(-1674902723.0, 288716400.0); 316 a[11][ 8] = fraction(54980371265625.0, 176692375811392.0); 317 a[11][ 9] = fraction(-734375.0, 4826304.0); 318 a[11][10] = fraction(171414593.0, 851261400.0); 319 a[11][11] = fraction(137909.0, 3084480.0); 320 321 // the following stages are for interpolation only 322 a[12][ 0] = fraction( 13481885573.0, 240030000000.0); 323 a[12][ 1] = getField().getZero(); 324 a[12][ 2] = getField().getZero(); 325 a[12][ 3] = getField().getZero(); 326 a[12][ 4] = getField().getZero(); 327 a[12][ 5] = getField().getZero(); 328 a[12][ 6] = fraction( 139418837528.0, 549975234375.0); 329 a[12][ 7] = fraction( -11108320068443.0, 45111937500000.0); 330 a[12][ 8] = fraction(-1769651421925959.0, 14249385146080000.0); 331 a[12][ 9] = fraction( 57799439.0, 377055000.0); 332 a[12][10] = fraction( 793322643029.0, 96734250000000.0); 333 a[12][11] = fraction( 1458939311.0, 192780000000.0); 334 a[12][12] = fraction( -4149.0, 500000.0); 335 336 a[13][ 0] = fraction( 1595561272731.0, 50120273500000.0); 337 a[13][ 1] = getField().getZero(); 338 a[13][ 2] = getField().getZero(); 339 a[13][ 3] = getField().getZero(); 340 a[13][ 4] = getField().getZero(); 341 a[13][ 5] = fraction( 975183916491.0, 34457688031250.0); 342 a[13][ 6] = fraction( 38492013932672.0, 718912673015625.0); 343 a[13][ 7] = fraction(-1114881286517557.0, 20298710767500000.0); 344 a[13][ 8] = getField().getZero(); 345 a[13][ 9] = getField().getZero(); 346 a[13][10] = fraction( -2538710946863.0, 23431227861250000.0); 347 a[13][11] = fraction( 8824659001.0, 23066716781250.0); 348 a[13][12] = fraction( -11518334563.0, 33831184612500.0); 349 a[13][13] = fraction( 1912306948.0, 13532473845.0); 350 351 a[14][ 0] = fraction( -13613986967.0, 31741908048.0); 352 a[14][ 1] = getField().getZero(); 353 a[14][ 2] = getField().getZero(); 354 a[14][ 3] = getField().getZero(); 355 a[14][ 4] = getField().getZero(); 356 a[14][ 5] = fraction( -4755612631.0, 1012344804.0); 357 a[14][ 6] = fraction( 42939257944576.0, 5588559685701.0); 358 a[14][ 7] = fraction( 77881972900277.0, 19140370552944.0); 359 a[14][ 8] = fraction( 22719829234375.0, 63689648654052.0); 360 a[14][ 9] = getField().getZero(); 361 a[14][10] = getField().getZero(); 362 a[14][11] = getField().getZero(); 363 a[14][12] = fraction( -1199007803.0, 857031517296.0); 364 a[14][13] = fraction( 157882067000.0, 53564469831.0); 365 a[14][14] = fraction( -290468882375.0, 31741908048.0); 366 367 return a; 368 } 369 370 /** {@inheritDoc} */ 371 @Override 372 public T[] getB() { 373 final T[] b = MathArrays.buildArray(getField(), 16); 374 b[ 0] = fraction(104257, 1920240); 375 b[ 1] = getField().getZero(); 376 b[ 2] = getField().getZero(); 377 b[ 3] = getField().getZero(); 378 b[ 4] = getField().getZero(); 379 b[ 5] = fraction( 3399327.0, 763840.0); 380 b[ 6] = fraction( 66578432.0, 35198415.0); 381 b[ 7] = fraction( -1674902723.0, 288716400.0); 382 b[ 8] = fraction( 54980371265625.0, 176692375811392.0); 383 b[ 9] = fraction( -734375.0, 4826304.0); 384 b[10] = fraction( 171414593.0, 851261400.0); 385 b[11] = fraction( 137909.0, 3084480.0); 386 b[12] = getField().getZero(); 387 b[13] = getField().getZero(); 388 b[14] = getField().getZero(); 389 b[15] = getField().getZero(); 390 return b; 391 } 392 393 /** {@inheritDoc} */ 394 @Override 395 protected DormandPrince853FieldStepInterpolator<T> 396 createInterpolator(final boolean forward, T[][] yDotK, 397 final FieldODEStateAndDerivative<T> globalPreviousState, 398 final FieldODEStateAndDerivative<T> globalCurrentState, final FieldEquationsMapper<T> mapper) { 399 return new DormandPrince853FieldStepInterpolator<>(getField(), forward, yDotK, 400 globalPreviousState, globalCurrentState, 401 globalPreviousState, globalCurrentState, 402 mapper); 403 } 404 405 /** {@inheritDoc} */ 406 @Override 407 public int getOrder() { 408 return 8; 409 } 410 411 /** {@inheritDoc} */ 412 @Override 413 protected T estimateError(final T[][] yDotK, final T[] y0, final T[] y1, final T h) { 414 T error1 = h.getField().getZero(); 415 T error2 = h.getField().getZero(); 416 417 for (int j = 0; j < mainSetDimension; ++j) { 418 final T errSum1 = yDotK[ 0][j].multiply(e1_01). 419 add(yDotK[ 5][j].multiply(e1_06)). 420 add(yDotK[ 6][j].multiply(e1_07)). 421 add(yDotK[ 7][j].multiply(e1_08)). 422 add(yDotK[ 8][j].multiply(e1_09)). 423 add(yDotK[ 9][j].multiply(e1_10)). 424 add(yDotK[10][j].multiply(e1_11)). 425 add(yDotK[11][j].multiply(e1_12)); 426 final T errSum2 = yDotK[ 0][j].multiply(e2_01). 427 add(yDotK[ 5][j].multiply(e2_06)). 428 add(yDotK[ 6][j].multiply(e2_07)). 429 add(yDotK[ 7][j].multiply(e2_08)). 430 add(yDotK[ 8][j].multiply(e2_09)). 431 add(yDotK[ 9][j].multiply(e2_10)). 432 add(yDotK[10][j].multiply(e2_11)). 433 add(yDotK[11][j].multiply(e2_12)); 434 435 final T yScale = RealFieldElement.max(y0[j].abs(), y1[j].abs()); 436 final T tol = vecAbsoluteTolerance == null ? 437 yScale.multiply(scalRelativeTolerance).add(scalAbsoluteTolerance) : 438 yScale.multiply(vecRelativeTolerance[j]).add(vecAbsoluteTolerance[j]); 439 final T ratio1 = errSum1.divide(tol); 440 error1 = error1.add(ratio1.multiply(ratio1)); 441 final T ratio2 = errSum2.divide(tol); 442 error2 = error2.add(ratio2.multiply(ratio2)); 443 } 444 445 T den = error1.add(error2.multiply(0.01)); 446 if (den.getReal() <= 0.0) { 447 den = h.getField().getOne(); 448 } 449 450 return h.abs().multiply(error1).divide(den.multiply(mainSetDimension).sqrt()); 451 } 452}