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.analysis.solvers.UnivariateSolver; 021import org.apache.commons.math4.legacy.exception.DimensionMismatchException; 022import org.apache.commons.math4.legacy.exception.MaxCountExceededException; 023import org.apache.commons.math4.legacy.exception.NoBracketingException; 024import org.apache.commons.math4.legacy.exception.NumberIsTooSmallException; 025import org.apache.commons.math4.legacy.ode.ExpandableStatefulODE; 026import org.apache.commons.math4.legacy.ode.events.EventHandler; 027import org.apache.commons.math4.legacy.ode.sampling.AbstractStepInterpolator; 028import org.apache.commons.math4.legacy.ode.sampling.StepHandler; 029import org.apache.commons.math4.core.jdkmath.JdkMath; 030 031/** 032 * This class implements a Gragg-Bulirsch-Stoer integrator for 033 * Ordinary Differential Equations. 034 * 035 * <p>The Gragg-Bulirsch-Stoer algorithm is one of the most efficient 036 * ones currently available for smooth problems. It uses Richardson 037 * extrapolation to estimate what would be the solution if the step 038 * size could be decreased down to zero.</p> 039 * 040 * <p> 041 * This method changes both the step size and the order during 042 * integration, in order to minimize computation cost. It is 043 * particularly well suited when a very high precision is needed. The 044 * limit where this method becomes more efficient than high-order 045 * embedded Runge-Kutta methods like {@link DormandPrince853Integrator 046 * Dormand-Prince 8(5,3)} depends on the problem. Results given in the 047 * Hairer, Norsett and Wanner book show for example that this limit 048 * occurs for accuracy around 1e-6 when integrating Saltzam-Lorenz 049 * equations (the authors note this problem is <i>extremely sensitive 050 * to the errors in the first integration steps</i>), and around 1e-11 051 * for a two dimensional celestial mechanics problems with seven 052 * bodies (pleiades problem, involving quasi-collisions for which 053 * <i>automatic step size control is essential</i>). 054 * </p> 055 * 056 * <p> 057 * This implementation is basically a reimplementation in Java of the 058 * <a 059 * href="http://www.unige.ch/math/folks/hairer/prog/nonstiff/odex.f">odex</a> 060 * fortran code by E. Hairer and G. Wanner. The redistribution policy 061 * for this code is available <a 062 * href="http://www.unige.ch/~hairer/prog/licence.txt">here</a>, for 063 * convenience, it is reproduced below.</p> 064 * 065 * <table border="" style="text-align: center; background-color: #E0E0E0"> 066 * <caption>odex redistribution policy</caption> 067 * <tr><td>Copyright (c) 2004, Ernst Hairer</td></tr> 068 * 069 * <tr><td>Redistribution and use in source and binary forms, with or 070 * without modification, are permitted provided that the following 071 * conditions are met: 072 * <ul> 073 * <li>Redistributions of source code must retain the above copyright 074 * notice, this list of conditions and the following disclaimer.</li> 075 * <li>Redistributions in binary form must reproduce the above copyright 076 * notice, this list of conditions and the following disclaimer in the 077 * documentation and/or other materials provided with the distribution.</li> 078 * </ul></td></tr> 079 * 080 * <tr><td><strong>THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND 081 * CONTRIBUTORS "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, 082 * BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS 083 * FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE REGENTS OR 084 * CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, 085 * EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, 086 * PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR 087 * PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF 088 * LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING 089 * NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS 090 * SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.</strong></td></tr> 091 * </table> 092 * 093 * @since 1.2 094 */ 095 096public class GraggBulirschStoerIntegrator extends AdaptiveStepsizeIntegrator { 097 098 /** Integrator method name. */ 099 private static final String METHOD_NAME = "Gragg-Bulirsch-Stoer"; 100 101 /** maximal order. */ 102 private int maxOrder; 103 104 /** step size sequence. */ 105 private int[] sequence; 106 107 /** overall cost of applying step reduction up to iteration k+1, in number of calls. */ 108 private int[] costPerStep; 109 110 /** cost per unit step. */ 111 private double[] costPerTimeUnit; 112 113 /** optimal steps for each order. */ 114 private double[] optimalStep; 115 116 /** extrapolation coefficients. */ 117 private double[][] coeff; 118 119 /** stability check enabling parameter. */ 120 private boolean performTest; 121 122 /** maximal number of checks for each iteration. */ 123 private int maxChecks; 124 125 /** maximal number of iterations for which checks are performed. */ 126 private int maxIter; 127 128 /** stepsize reduction factor in case of stability check failure. */ 129 private double stabilityReduction; 130 131 /** first stepsize control factor. */ 132 private double stepControl1; 133 134 /** second stepsize control factor. */ 135 private double stepControl2; 136 137 /** third stepsize control factor. */ 138 private double stepControl3; 139 140 /** fourth stepsize control factor. */ 141 private double stepControl4; 142 143 /** first order control factor. */ 144 private double orderControl1; 145 146 /** second order control factor. */ 147 private double orderControl2; 148 149 /** use interpolation error in stepsize control. */ 150 private boolean useInterpolationError; 151 152 /** interpolation order control parameter. */ 153 private int mudif; 154 155 /** Simple constructor. 156 * Build a Gragg-Bulirsch-Stoer integrator with the given step 157 * bounds. All tuning parameters are set to their default 158 * values. The default step handler does nothing. 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 scalAbsoluteTolerance allowed absolute error 166 * @param scalRelativeTolerance allowed relative error 167 */ 168 public GraggBulirschStoerIntegrator(final double minStep, final double maxStep, 169 final double scalAbsoluteTolerance, 170 final double scalRelativeTolerance) { 171 super(METHOD_NAME, minStep, maxStep, 172 scalAbsoluteTolerance, scalRelativeTolerance); 173 setStabilityCheck(true, -1, -1, -1); 174 setControlFactors(-1, -1, -1, -1); 175 setOrderControl(-1, -1, -1); 176 setInterpolationControl(true, -1); 177 } 178 179 /** Simple constructor. 180 * Build a Gragg-Bulirsch-Stoer integrator with the given step 181 * bounds. All tuning parameters are set to their default 182 * values. The default step handler does nothing. 183 * @param minStep minimal step (must be positive even for backward 184 * integration), the last step can be smaller than this 185 * @param maxStep maximal step (must be positive even for backward 186 * integration) 187 * @param vecAbsoluteTolerance allowed absolute error 188 * @param vecRelativeTolerance allowed relative error 189 */ 190 public GraggBulirschStoerIntegrator(final double minStep, final double maxStep, 191 final double[] vecAbsoluteTolerance, 192 final double[] vecRelativeTolerance) { 193 super(METHOD_NAME, minStep, maxStep, 194 vecAbsoluteTolerance, vecRelativeTolerance); 195 setStabilityCheck(true, -1, -1, -1); 196 setControlFactors(-1, -1, -1, -1); 197 setOrderControl(-1, -1, -1); 198 setInterpolationControl(true, -1); 199 } 200 201 /** Set the stability check controls. 202 * <p>The stability check is performed on the first few iterations of 203 * the extrapolation scheme. If this test fails, the step is rejected 204 * and the stepsize is reduced.</p> 205 * <p>By default, the test is performed, at most during two 206 * iterations at each step, and at most once for each of these 207 * iterations. The default stepsize reduction factor is 0.5.</p> 208 * @param performStabilityCheck if true, stability check will be performed, 209 if false, the check will be skipped 210 * @param maxNumIter maximal number of iterations for which checks are 211 * performed (the number of iterations is reset to default if negative 212 * or null) 213 * @param maxNumChecks maximal number of checks for each iteration 214 * (the number of checks is reset to default if negative or null) 215 * @param stepsizeReductionFactor stepsize reduction factor in case of 216 * failure (the factor is reset to default if lower than 0.0001 or 217 * greater than 0.9999) 218 */ 219 public void setStabilityCheck(final boolean performStabilityCheck, 220 final int maxNumIter, final int maxNumChecks, 221 final double stepsizeReductionFactor) { 222 223 this.performTest = performStabilityCheck; 224 this.maxIter = (maxNumIter <= 0) ? 2 : maxNumIter; 225 this.maxChecks = (maxNumChecks <= 0) ? 1 : maxNumChecks; 226 227 if (stepsizeReductionFactor < 0.0001 || stepsizeReductionFactor > 0.9999) { 228 this.stabilityReduction = 0.5; 229 } else { 230 this.stabilityReduction = stepsizeReductionFactor; 231 } 232 } 233 234 /** Set the step size control factors. 235 236 * <p>The new step size hNew is computed from the old one h by: 237 * <pre> 238 * hNew = h * stepControl2 / (err/stepControl1)^(1/(2k+1)) 239 * </pre> 240 * where err is the scaled error and k the iteration number of the 241 * extrapolation scheme (counting from 0). The default values are 242 * 0.65 for stepControl1 and 0.94 for stepControl2. 243 * <p>The step size is subject to the restriction: 244 * <pre> 245 * stepControl3^(1/(2k+1))/stepControl4 <= hNew/h <= 1/stepControl3^(1/(2k+1)) 246 * </pre> 247 * The default values are 0.02 for stepControl3 and 4.0 for 248 * stepControl4. 249 * @param control1 first stepsize control factor (the factor is 250 * reset to default if lower than 0.0001 or greater than 0.9999) 251 * @param control2 second stepsize control factor (the factor 252 * is reset to default if lower than 0.0001 or greater than 0.9999) 253 * @param control3 third stepsize control factor (the factor is 254 * reset to default if lower than 0.0001 or greater than 0.9999) 255 * @param control4 fourth stepsize control factor (the factor 256 * is reset to default if lower than 1.0001 or greater than 999.9) 257 */ 258 public void setControlFactors(final double control1, final double control2, 259 final double control3, final double control4) { 260 261 if (control1 < 0.0001 || control1 > 0.9999) { 262 this.stepControl1 = 0.65; 263 } else { 264 this.stepControl1 = control1; 265 } 266 267 if (control2 < 0.0001 || control2 > 0.9999) { 268 this.stepControl2 = 0.94; 269 } else { 270 this.stepControl2 = control2; 271 } 272 273 if (control3 < 0.0001 || control3 > 0.9999) { 274 this.stepControl3 = 0.02; 275 } else { 276 this.stepControl3 = control3; 277 } 278 279 if (control4 < 1.0001 || control4 > 999.9) { 280 this.stepControl4 = 4.0; 281 } else { 282 this.stepControl4 = control4; 283 } 284 } 285 286 /** Set the order control parameters. 287 * <p>The Gragg-Bulirsch-Stoer method changes both the step size and 288 * the order during integration, in order to minimize computation 289 * cost. Each extrapolation step increases the order by 2, so the 290 * maximal order that will be used is always even, it is twice the 291 * maximal number of columns in the extrapolation table.</p> 292 * <pre> 293 * order is decreased if w(k-1) <= w(k) * orderControl1 294 * order is increased if w(k) <= w(k-1) * orderControl2 295 * </pre> 296 * <p>where w is the table of work per unit step for each order 297 * (number of function calls divided by the step length), and k is 298 * the current order.</p> 299 * <p>The default maximal order after construction is 18 (i.e. the 300 * maximal number of columns is 9). The default values are 0.8 for 301 * orderControl1 and 0.9 for orderControl2.</p> 302 * @param maximalOrder maximal order in the extrapolation table (the 303 * maximal order is reset to default if order <= 6 or odd) 304 * @param control1 first order control factor (the factor is 305 * reset to default if lower than 0.0001 or greater than 0.9999) 306 * @param control2 second order control factor (the factor 307 * is reset to default if lower than 0.0001 or greater than 0.9999) 308 */ 309 public void setOrderControl(final int maximalOrder, 310 final double control1, final double control2) { 311 312 if (maximalOrder <= 6 || (maximalOrder & 1) != 0) { 313 this.maxOrder = 18; 314 } 315 316 if (control1 < 0.0001 || control1 > 0.9999) { 317 this.orderControl1 = 0.8; 318 } else { 319 this.orderControl1 = control1; 320 } 321 322 if (control2 < 0.0001 || control2 > 0.9999) { 323 this.orderControl2 = 0.9; 324 } else { 325 this.orderControl2 = control2; 326 } 327 328 // reinitialize the arrays 329 initializeArrays(); 330 } 331 332 /** {@inheritDoc} */ 333 @Override 334 public void addStepHandler (final StepHandler handler) { 335 336 super.addStepHandler(handler); 337 338 // reinitialize the arrays 339 initializeArrays(); 340 } 341 342 /** {@inheritDoc} */ 343 @Override 344 public void addEventHandler(final EventHandler function, 345 final double maxCheckInterval, 346 final double convergence, 347 final int maxIterationCount, 348 final UnivariateSolver solver) { 349 super.addEventHandler(function, maxCheckInterval, convergence, 350 maxIterationCount, solver); 351 352 // reinitialize the arrays 353 initializeArrays(); 354 } 355 356 /** Initialize the integrator internal arrays. */ 357 private void initializeArrays() { 358 359 final int size = maxOrder / 2; 360 361 if (sequence == null || sequence.length != size) { 362 // all arrays should be reallocated with the right size 363 sequence = new int[size]; 364 costPerStep = new int[size]; 365 coeff = new double[size][]; 366 costPerTimeUnit = new double[size]; 367 optimalStep = new double[size]; 368 } 369 370 // step size sequence: 2, 6, 10, 14, ... 371 for (int k = 0; k < size; ++k) { 372 sequence[k] = 4 * k + 2; 373 } 374 375 // initialize the order selection cost array 376 // (number of function calls for each column of the extrapolation table) 377 costPerStep[0] = sequence[0] + 1; 378 for (int k = 1; k < size; ++k) { 379 costPerStep[k] = costPerStep[k-1] + sequence[k]; 380 } 381 382 // initialize the extrapolation tables 383 for (int k = 0; k < size; ++k) { 384 coeff[k] = (k > 0) ? new double[k] : null; 385 for (int l = 0; l < k; ++l) { 386 final double ratio = ((double) sequence[k]) / sequence[k-l-1]; 387 coeff[k][l] = 1.0 / (ratio * ratio - 1.0); 388 } 389 } 390 } 391 392 /** Set the interpolation order control parameter. 393 * The interpolation order for dense output is 2k - mudif + 1. The 394 * default value for mudif is 4 and the interpolation error is used 395 * in stepsize control by default. 396 397 * @param useInterpolationErrorForControl if true, interpolation error is used 398 * for stepsize control 399 * @param mudifControlParameter interpolation order control parameter (the parameter 400 * is reset to default if <= 0 or >= 7) 401 */ 402 public void setInterpolationControl(final boolean useInterpolationErrorForControl, 403 final int mudifControlParameter) { 404 405 this.useInterpolationError = useInterpolationErrorForControl; 406 407 if (mudifControlParameter <= 0 || mudifControlParameter >= 7) { 408 this.mudif = 4; 409 } else { 410 this.mudif = mudifControlParameter; 411 } 412 } 413 414 /** Update scaling array. 415 * @param y1 first state vector to use for scaling 416 * @param y2 second state vector to use for scaling 417 * @param scale scaling array to update (can be shorter than state) 418 */ 419 private void rescale(final double[] y1, final double[] y2, final double[] scale) { 420 if (vecAbsoluteTolerance == null) { 421 for (int i = 0; i < scale.length; ++i) { 422 final double yi = JdkMath.max(JdkMath.abs(y1[i]), JdkMath.abs(y2[i])); 423 scale[i] = scalAbsoluteTolerance + scalRelativeTolerance * yi; 424 } 425 } else { 426 for (int i = 0; i < scale.length; ++i) { 427 final double yi = JdkMath.max(JdkMath.abs(y1[i]), JdkMath.abs(y2[i])); 428 scale[i] = vecAbsoluteTolerance[i] + vecRelativeTolerance[i] * yi; 429 } 430 } 431 } 432 433 /** Perform integration over one step using substeps of a modified 434 * midpoint method. 435 * @param t0 initial time 436 * @param y0 initial value of the state vector at t0 437 * @param step global step 438 * @param k iteration number (from 0 to sequence.length - 1) 439 * @param scale scaling array (can be shorter than state) 440 * @param f placeholder where to put the state vector derivatives at each substep 441 * (element 0 already contains initial derivative) 442 * @param yMiddle placeholder where to put the state vector at the middle of the step 443 * @param yEnd placeholder where to put the state vector at the end 444 * @param yTmp placeholder for one state vector 445 * @return true if computation was done properly, 446 * false if stability check failed before end of computation 447 * @exception MaxCountExceededException if the number of functions evaluations is exceeded 448 * @exception DimensionMismatchException if arrays dimensions do not match equations settings 449 */ 450 private boolean tryStep(final double t0, final double[] y0, final double step, final int k, 451 final double[] scale, final double[][] f, 452 final double[] yMiddle, final double[] yEnd, 453 final double[] yTmp) 454 throws MaxCountExceededException, DimensionMismatchException { 455 456 final int n = sequence[k]; 457 final double subStep = step / n; 458 final double subStep2 = 2 * subStep; 459 460 // first substep 461 double t = t0 + subStep; 462 for (int i = 0; i < y0.length; ++i) { 463 yTmp[i] = y0[i]; 464 yEnd[i] = y0[i] + subStep * f[0][i]; 465 } 466 computeDerivatives(t, yEnd, f[1]); 467 468 // other substeps 469 for (int j = 1; j < n; ++j) { 470 471 if (2 * j == n) { 472 // save the point at the middle of the step 473 System.arraycopy(yEnd, 0, yMiddle, 0, y0.length); 474 } 475 476 t += subStep; 477 for (int i = 0; i < y0.length; ++i) { 478 final double middle = yEnd[i]; 479 yEnd[i] = yTmp[i] + subStep2 * f[j][i]; 480 yTmp[i] = middle; 481 } 482 483 computeDerivatives(t, yEnd, f[j+1]); 484 485 // stability check 486 if (performTest && j <= maxChecks && k < maxIter) { 487 double initialNorm = 0.0; 488 for (int l = 0; l < scale.length; ++l) { 489 final double ratio = f[0][l] / scale[l]; 490 initialNorm += ratio * ratio; 491 } 492 double deltaNorm = 0.0; 493 for (int l = 0; l < scale.length; ++l) { 494 final double ratio = (f[j+1][l] - f[0][l]) / scale[l]; 495 deltaNorm += ratio * ratio; 496 } 497 if (deltaNorm > 4 * JdkMath.max(1.0e-15, initialNorm)) { 498 return false; 499 } 500 } 501 } 502 503 // correction of the last substep (at t0 + step) 504 for (int i = 0; i < y0.length; ++i) { 505 yEnd[i] = 0.5 * (yTmp[i] + yEnd[i] + subStep * f[n][i]); 506 } 507 508 return true; 509 } 510 511 /** Extrapolate a vector. 512 * @param offset offset to use in the coefficients table 513 * @param k index of the last updated point 514 * @param diag working diagonal of the Aitken-Neville's 515 * triangle, without the last element 516 * @param last last element 517 */ 518 private void extrapolate(final int offset, final int k, 519 final double[][] diag, final double[] last) { 520 521 // update the diagonal 522 for (int j = 1; j < k; ++j) { 523 for (int i = 0; i < last.length; ++i) { 524 // Aitken-Neville's recursive formula 525 diag[k-j-1][i] = diag[k-j][i] + 526 coeff[k+offset][j-1] * (diag[k-j][i] - diag[k-j-1][i]); 527 } 528 } 529 530 // update the last element 531 for (int i = 0; i < last.length; ++i) { 532 // Aitken-Neville's recursive formula 533 last[i] = diag[0][i] + coeff[k+offset][k-1] * (diag[0][i] - last[i]); 534 } 535 } 536 537 /** {@inheritDoc} */ 538 @Override 539 public void integrate(final ExpandableStatefulODE equations, final double t) 540 throws NumberIsTooSmallException, DimensionMismatchException, 541 MaxCountExceededException, NoBracketingException { 542 543 sanityChecks(equations, t); 544 setEquations(equations); 545 final boolean forward = t > equations.getTime(); 546 547 // create some internal working arrays 548 final double[] y0 = equations.getCompleteState(); 549 final double[] y = y0.clone(); 550 final double[] yDot0 = new double[y.length]; 551 final double[] y1 = new double[y.length]; 552 final double[] yTmp = new double[y.length]; 553 final double[] yTmpDot = new double[y.length]; 554 555 final double[][] diagonal = new double[sequence.length-1][]; 556 final double[][] y1Diag = new double[sequence.length-1][]; 557 for (int k = 0; k < sequence.length-1; ++k) { 558 diagonal[k] = new double[y.length]; 559 y1Diag[k] = new double[y.length]; 560 } 561 562 final double[][][] fk = new double[sequence.length][][]; 563 for (int k = 0; k < sequence.length; ++k) { 564 565 fk[k] = new double[sequence[k] + 1][]; 566 567 // all substeps start at the same point, so share the first array 568 fk[k][0] = yDot0; 569 570 for (int l = 0; l < sequence[k]; ++l) { 571 fk[k][l+1] = new double[y0.length]; 572 } 573 } 574 575 if (y != y0) { 576 System.arraycopy(y0, 0, y, 0, y0.length); 577 } 578 579 final double[] yDot1 = new double[y0.length]; 580 final double[][] yMidDots = new double[1 + 2 * sequence.length][y0.length]; 581 582 // initial scaling 583 final double[] scale = new double[mainSetDimension]; 584 rescale(y, y, scale); 585 586 // initial order selection 587 final double tol = 588 (vecRelativeTolerance == null) ? scalRelativeTolerance : vecRelativeTolerance[0]; 589 final double log10R = JdkMath.log10(JdkMath.max(1.0e-10, tol)); 590 int targetIter = JdkMath.max(1, 591 JdkMath.min(sequence.length - 2, 592 (int) JdkMath.floor(0.5 - 0.6 * log10R))); 593 594 // set up an interpolator sharing the integrator arrays 595 final AbstractStepInterpolator interpolator = 596 new GraggBulirschStoerStepInterpolator(y, yDot0, 597 y1, yDot1, 598 yMidDots, forward, 599 equations.getPrimaryMapper(), 600 equations.getSecondaryMappers()); 601 interpolator.storeTime(equations.getTime()); 602 603 stepStart = equations.getTime(); 604 double hNew = 0; 605 double maxError = Double.MAX_VALUE; 606 boolean previousRejected = false; 607 boolean firstTime = true; 608 boolean newStep = true; 609 boolean firstStepAlreadyComputed = false; 610 initIntegration(equations.getTime(), y0, t); 611 costPerTimeUnit[0] = 0; 612 isLastStep = false; 613 do { 614 615 double error; 616 boolean reject = false; 617 618 if (newStep) { 619 620 interpolator.shift(); 621 622 // first evaluation, at the beginning of the step 623 if (! firstStepAlreadyComputed) { 624 computeDerivatives(stepStart, y, yDot0); 625 } 626 627 if (firstTime) { 628 hNew = initializeStep(forward, 2 * targetIter + 1, scale, 629 stepStart, y, yDot0, yTmp, yTmpDot); 630 } 631 632 newStep = false; 633 } 634 635 stepSize = hNew; 636 637 // step adjustment near bounds 638 if ((forward && (stepStart + stepSize > t)) || 639 ((! forward) && (stepStart + stepSize < t))) { 640 stepSize = t - stepStart; 641 } 642 final double nextT = stepStart + stepSize; 643 isLastStep = forward ? (nextT >= t) : (nextT <= t); 644 645 // iterate over several substep sizes 646 int k = -1; 647 for (boolean loop = true; loop;) { 648 649 ++k; 650 651 // modified midpoint integration with the current substep 652 if ( ! tryStep(stepStart, y, stepSize, k, scale, fk[k], 653 (k == 0) ? yMidDots[0] : diagonal[k-1], 654 (k == 0) ? y1 : y1Diag[k-1], 655 yTmp)) { 656 657 // the stability check failed, we reduce the global step 658 hNew = JdkMath.abs(filterStep(stepSize * stabilityReduction, forward, false)); 659 reject = true; 660 loop = false; 661 } else { 662 663 // the substep was computed successfully 664 if (k > 0) { 665 666 // extrapolate the state at the end of the step 667 // using last iteration data 668 extrapolate(0, k, y1Diag, y1); 669 rescale(y, y1, scale); 670 671 // estimate the error at the end of the step. 672 error = 0; 673 for (int j = 0; j < mainSetDimension; ++j) { 674 final double e = JdkMath.abs(y1[j] - y1Diag[0][j]) / scale[j]; 675 error += e * e; 676 } 677 error = JdkMath.sqrt(error / mainSetDimension); 678 679 if (error > 1.0e15 || (k > 1 && error > maxError)) { 680 // error is too big, we reduce the global step 681 hNew = JdkMath.abs(filterStep(stepSize * stabilityReduction, forward, false)); 682 reject = true; 683 loop = false; 684 } else { 685 686 maxError = JdkMath.max(4 * error, 1.0); 687 688 // compute optimal stepsize for this order 689 final double exp = 1.0 / (2 * k + 1); 690 double fac = stepControl2 / JdkMath.pow(error / stepControl1, exp); 691 final double pow = JdkMath.pow(stepControl3, exp); 692 fac = JdkMath.max(pow / stepControl4, JdkMath.min(1 / pow, fac)); 693 optimalStep[k] = JdkMath.abs(filterStep(stepSize * fac, forward, true)); 694 costPerTimeUnit[k] = costPerStep[k] / optimalStep[k]; 695 696 // check convergence 697 switch (k - targetIter) { 698 699 case -1 : 700 if (targetIter > 1 && !previousRejected) { 701 702 // check if we can stop iterations now 703 if (error <= 1.0) { 704 // convergence have been reached just before targetIter 705 loop = false; 706 } else { 707 // estimate if there is a chance convergence will 708 // be reached on next iteration, using the 709 // asymptotic evolution of error 710 final double ratio = ((double) sequence [targetIter] * sequence[targetIter + 1]) / 711 (sequence[0] * sequence[0]); 712 if (error > ratio * ratio) { 713 // we don't expect to converge on next iteration 714 // we reject the step immediately and reduce order 715 reject = true; 716 loop = false; 717 targetIter = k; 718 if (targetIter > 1 && 719 costPerTimeUnit[targetIter-1] < 720 orderControl1 * costPerTimeUnit[targetIter]) { 721 --targetIter; 722 } 723 hNew = optimalStep[targetIter]; 724 } 725 } 726 } 727 break; 728 729 case 0: 730 if (error <= 1.0) { 731 // convergence has been reached exactly at targetIter 732 loop = false; 733 } else { 734 // estimate if there is a chance convergence will 735 // be reached on next iteration, using the 736 // asymptotic evolution of error 737 final double ratio = ((double) sequence[k+1]) / sequence[0]; 738 if (error > ratio * ratio) { 739 // we don't expect to converge on next iteration 740 // we reject the step immediately 741 reject = true; 742 loop = false; 743 if (targetIter > 1 && 744 costPerTimeUnit[targetIter-1] < 745 orderControl1 * costPerTimeUnit[targetIter]) { 746 --targetIter; 747 } 748 hNew = optimalStep[targetIter]; 749 } 750 } 751 break; 752 753 case 1 : 754 if (error > 1.0) { 755 reject = true; 756 if (targetIter > 1 && 757 costPerTimeUnit[targetIter-1] < 758 orderControl1 * costPerTimeUnit[targetIter]) { 759 --targetIter; 760 } 761 hNew = optimalStep[targetIter]; 762 } 763 loop = false; 764 break; 765 766 default : 767 if ((firstTime || isLastStep) && error <= 1.0) { 768 loop = false; 769 } 770 break; 771 } 772 } 773 } 774 } 775 } 776 777 if (! reject) { 778 // derivatives at end of step 779 computeDerivatives(stepStart + stepSize, y1, yDot1); 780 } 781 782 // dense output handling 783 double hInt = getMaxStep(); 784 if (! reject) { 785 786 // extrapolate state at middle point of the step 787 for (int j = 1; j <= k; ++j) { 788 extrapolate(0, j, diagonal, yMidDots[0]); 789 } 790 791 final int mu = 2 * k - mudif + 3; 792 793 for (int l = 0; l < mu; ++l) { 794 795 // derivative at middle point of the step 796 final int l2 = l / 2; 797 double factor = JdkMath.pow(0.5 * sequence[l2], l); 798 int middleIndex = fk[l2].length / 2; 799 for (int i = 0; i < y0.length; ++i) { 800 yMidDots[l+1][i] = factor * fk[l2][middleIndex + l][i]; 801 } 802 for (int j = 1; j <= k - l2; ++j) { 803 factor = JdkMath.pow(0.5 * sequence[j + l2], l); 804 middleIndex = fk[l2+j].length / 2; 805 for (int i = 0; i < y0.length; ++i) { 806 diagonal[j-1][i] = factor * fk[l2+j][middleIndex+l][i]; 807 } 808 extrapolate(l2, j, diagonal, yMidDots[l+1]); 809 } 810 for (int i = 0; i < y0.length; ++i) { 811 yMidDots[l+1][i] *= stepSize; 812 } 813 814 // compute centered differences to evaluate next derivatives 815 for (int j = (l + 1) / 2; j <= k; ++j) { 816 for (int m = fk[j].length - 1; m >= 2 * (l + 1); --m) { 817 for (int i = 0; i < y0.length; ++i) { 818 fk[j][m][i] -= fk[j][m-2][i]; 819 } 820 } 821 } 822 } 823 824 if (mu >= 0) { 825 826 // estimate the dense output coefficients 827 final GraggBulirschStoerStepInterpolator gbsInterpolator 828 = (GraggBulirschStoerStepInterpolator) interpolator; 829 gbsInterpolator.computeCoefficients(mu, stepSize); 830 831 if (useInterpolationError) { 832 // use the interpolation error to limit stepsize 833 final double interpError = gbsInterpolator.estimateError(scale); 834 hInt = JdkMath.abs(stepSize / JdkMath.max(JdkMath.pow(interpError, 1.0 / (mu+4)), 835 0.01)); 836 if (interpError > 10.0) { 837 hNew = hInt; 838 reject = true; 839 } 840 } 841 } 842 } 843 844 if (! reject) { 845 846 // Discrete events handling 847 interpolator.storeTime(stepStart + stepSize); 848 stepStart = acceptStep(interpolator, y1, yDot1, t); 849 850 // prepare next step 851 interpolator.storeTime(stepStart); 852 System.arraycopy(y1, 0, y, 0, y0.length); 853 System.arraycopy(yDot1, 0, yDot0, 0, y0.length); 854 firstStepAlreadyComputed = true; 855 856 int optimalIter; 857 if (k == 1) { 858 optimalIter = 2; 859 if (previousRejected) { 860 optimalIter = 1; 861 } 862 } else if (k <= targetIter) { 863 optimalIter = k; 864 if (costPerTimeUnit[k-1] < orderControl1 * costPerTimeUnit[k]) { 865 optimalIter = k-1; 866 } else if (costPerTimeUnit[k] < orderControl2 * costPerTimeUnit[k-1]) { 867 optimalIter = JdkMath.min(k+1, sequence.length - 2); 868 } 869 } else { 870 optimalIter = k - 1; 871 if (k > 2 && 872 costPerTimeUnit[k-2] < orderControl1 * costPerTimeUnit[k-1]) { 873 optimalIter = k - 2; 874 } 875 if (costPerTimeUnit[k] < orderControl2 * costPerTimeUnit[optimalIter]) { 876 optimalIter = JdkMath.min(k, sequence.length - 2); 877 } 878 } 879 880 if (previousRejected) { 881 // after a rejected step neither order nor stepsize 882 // should increase 883 targetIter = JdkMath.min(optimalIter, k); 884 hNew = JdkMath.min(JdkMath.abs(stepSize), optimalStep[targetIter]); 885 } else { 886 // stepsize control 887 if (optimalIter <= k) { 888 hNew = optimalStep[optimalIter]; 889 } else { 890 if (k < targetIter && 891 costPerTimeUnit[k] < orderControl2 * costPerTimeUnit[k-1]) { 892 hNew = filterStep(optimalStep[k] * costPerStep[optimalIter+1] / costPerStep[k], 893 forward, false); 894 } else { 895 hNew = filterStep(optimalStep[k] * costPerStep[optimalIter] / costPerStep[k], 896 forward, false); 897 } 898 } 899 900 targetIter = optimalIter; 901 } 902 903 newStep = true; 904 } 905 906 hNew = JdkMath.min(hNew, hInt); 907 if (! forward) { 908 hNew = -hNew; 909 } 910 911 firstTime = false; 912 913 if (reject) { 914 isLastStep = false; 915 previousRejected = true; 916 } else { 917 previousRejected = false; 918 } 919 } while (!isLastStep); 920 921 // dispatch results 922 equations.setTime(stepStart); 923 equations.setCompleteState(y); 924 925 resetInternalState(); 926 } 927}