1dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond/* 2dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * Licensed to the Apache Software Foundation (ASF) under one or more 3dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * contributor license agreements. See the NOTICE file distributed with 4dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * this work for additional information regarding copyright ownership. 5dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * The ASF licenses this file to You under the Apache License, Version 2.0 6dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * (the "License"); you may not use this file except in compliance with 7dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * the License. You may obtain a copy of the License at 8dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * 9dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * http://www.apache.org/licenses/LICENSE-2.0 10dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * 11dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * Unless required by applicable law or agreed to in writing, software 12dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * distributed under the License is distributed on an "AS IS" BASIS, 13dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. 14dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * See the License for the specific language governing permissions and 15dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * limitations under the License. 16dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond */ 17dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 18dee0849a9704d532af0b550146cbafbaa6ee1d19Raymondpackage org.apache.commons.math.ode.nonstiff; 19dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 20dee0849a9704d532af0b550146cbafbaa6ee1d19Raymondimport org.apache.commons.math.ode.DerivativeException; 21dee0849a9704d532af0b550146cbafbaa6ee1d19Raymondimport org.apache.commons.math.ode.FirstOrderDifferentialEquations; 22dee0849a9704d532af0b550146cbafbaa6ee1d19Raymondimport org.apache.commons.math.ode.IntegratorException; 23dee0849a9704d532af0b550146cbafbaa6ee1d19Raymondimport org.apache.commons.math.ode.sampling.AbstractStepInterpolator; 24dee0849a9704d532af0b550146cbafbaa6ee1d19Raymondimport org.apache.commons.math.ode.sampling.DummyStepInterpolator; 25dee0849a9704d532af0b550146cbafbaa6ee1d19Raymondimport org.apache.commons.math.ode.sampling.StepHandler; 26dee0849a9704d532af0b550146cbafbaa6ee1d19Raymondimport org.apache.commons.math.util.FastMath; 27dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 28dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond/** 29dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * This class implements the common part of all embedded Runge-Kutta 30dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * integrators for Ordinary Differential Equations. 31dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * 32dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * <p>These methods are embedded explicit Runge-Kutta methods with two 33dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * sets of coefficients allowing to estimate the error, their Butcher 34dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * arrays are as follows : 35dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * <pre> 36dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * 0 | 37dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * c2 | a21 38dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * c3 | a31 a32 39dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * ... | ... 40dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * cs | as1 as2 ... ass-1 41dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * |-------------------------- 42dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * | b1 b2 ... bs-1 bs 43dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * | b'1 b'2 ... b's-1 b's 44dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * </pre> 45dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * </p> 46dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * 47dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * <p>In fact, we rather use the array defined by ej = bj - b'j to 48dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * compute directly the error rather than computing two estimates and 49dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * then comparing them.</p> 50dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * 51dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * <p>Some methods are qualified as <i>fsal</i> (first same as last) 52dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * methods. This means the last evaluation of the derivatives in one 53dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * step is the same as the first in the next step. Then, this 54dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * evaluation can be reused from one step to the next one and the cost 55dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * of such a method is really s-1 evaluations despite the method still 56dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * has s stages. This behaviour is true only for successful steps, if 57dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * the step is rejected after the error estimation phase, no 58dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * evaluation is saved. For an <i>fsal</i> method, we have cs = 1 and 59dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * asi = bi for all i.</p> 60dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * 61dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * @version $Revision: 1073158 $ $Date: 2011-02-21 22:46:52 +0100 (lun. 21 févr. 2011) $ 62dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * @since 1.2 63dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond */ 64dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 65dee0849a9704d532af0b550146cbafbaa6ee1d19Raymondpublic abstract class EmbeddedRungeKuttaIntegrator 66dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond extends AdaptiveStepsizeIntegrator { 67dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 68dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond /** Indicator for <i>fsal</i> methods. */ 69dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond private final boolean fsal; 70dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 71dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond /** Time steps from Butcher array (without the first zero). */ 72dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond private final double[] c; 73dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 74dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond /** Internal weights from Butcher array (without the first empty row). */ 75dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond private final double[][] a; 76dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 77dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond /** External weights for the high order method from Butcher array. */ 78dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond private final double[] b; 79dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 80dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond /** Prototype of the step interpolator. */ 81dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond private final RungeKuttaStepInterpolator prototype; 82dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 83dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond /** Stepsize control exponent. */ 84dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond private final double exp; 85dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 86dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond /** Safety factor for stepsize control. */ 87dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond private double safety; 88dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 89dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond /** Minimal reduction factor for stepsize control. */ 90dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond private double minReduction; 91dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 92dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond /** Maximal growth factor for stepsize control. */ 93dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond private double maxGrowth; 94dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 95dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond /** Build a Runge-Kutta integrator with the given Butcher array. 96dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * @param name name of the method 97dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * @param fsal indicate that the method is an <i>fsal</i> 98dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * @param c time steps from Butcher array (without the first zero) 99dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * @param a internal weights from Butcher array (without the first empty row) 100dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * @param b propagation weights for the high order method from Butcher array 101dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * @param prototype prototype of the step interpolator to use 102dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * @param minStep minimal step (must be positive even for backward 103dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * integration), the last step can be smaller than this 104dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * @param maxStep maximal step (must be positive even for backward 105dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * integration) 106dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * @param scalAbsoluteTolerance allowed absolute error 107dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * @param scalRelativeTolerance allowed relative error 108dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond */ 109dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond protected EmbeddedRungeKuttaIntegrator(final String name, final boolean fsal, 110dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond final double[] c, final double[][] a, final double[] b, 111dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond final RungeKuttaStepInterpolator prototype, 112dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond final double minStep, final double maxStep, 113dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond final double scalAbsoluteTolerance, 114dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond final double scalRelativeTolerance) { 115dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 116dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond super(name, minStep, maxStep, scalAbsoluteTolerance, scalRelativeTolerance); 117dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 118dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond this.fsal = fsal; 119dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond this.c = c; 120dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond this.a = a; 121dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond this.b = b; 122dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond this.prototype = prototype; 123dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 124dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond exp = -1.0 / getOrder(); 125dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 126dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond // set the default values of the algorithm control parameters 127dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond setSafety(0.9); 128dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond setMinReduction(0.2); 129dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond setMaxGrowth(10.0); 130dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 131dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond } 132dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 133dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond /** Build a Runge-Kutta integrator with the given Butcher array. 134dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * @param name name of the method 135dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * @param fsal indicate that the method is an <i>fsal</i> 136dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * @param c time steps from Butcher array (without the first zero) 137dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * @param a internal weights from Butcher array (without the first empty row) 138dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * @param b propagation weights for the high order method from Butcher array 139dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * @param prototype prototype of the step interpolator to use 140dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * @param minStep minimal step (must be positive even for backward 141dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * integration), the last step can be smaller than this 142dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * @param maxStep maximal step (must be positive even for backward 143dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * integration) 144dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * @param vecAbsoluteTolerance allowed absolute error 145dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * @param vecRelativeTolerance allowed relative error 146dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond */ 147dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond protected EmbeddedRungeKuttaIntegrator(final String name, final boolean fsal, 148dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond final double[] c, final double[][] a, final double[] b, 149dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond final RungeKuttaStepInterpolator prototype, 150dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond final double minStep, final double maxStep, 151dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond final double[] vecAbsoluteTolerance, 152dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond final double[] vecRelativeTolerance) { 153dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 154dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond super(name, minStep, maxStep, vecAbsoluteTolerance, vecRelativeTolerance); 155dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 156dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond this.fsal = fsal; 157dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond this.c = c; 158dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond this.a = a; 159dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond this.b = b; 160dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond this.prototype = prototype; 161dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 162dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond exp = -1.0 / getOrder(); 163dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 164dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond // set the default values of the algorithm control parameters 165dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond setSafety(0.9); 166dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond setMinReduction(0.2); 167dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond setMaxGrowth(10.0); 168dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 169dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond } 170dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 171dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond /** Get the order of the method. 172dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * @return order of the method 173dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond */ 174dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond public abstract int getOrder(); 175dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 176dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond /** Get the safety factor for stepsize control. 177dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * @return safety factor 178dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond */ 179dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond public double getSafety() { 180dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond return safety; 181dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond } 182dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 183dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond /** Set the safety factor for stepsize control. 184dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * @param safety safety factor 185dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond */ 186dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond public void setSafety(final double safety) { 187dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond this.safety = safety; 188dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond } 189dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 190dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond /** {@inheritDoc} */ 191dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond @Override 192dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond public double integrate(final FirstOrderDifferentialEquations equations, 193dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond final double t0, final double[] y0, 194dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond final double t, final double[] y) 195dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond throws DerivativeException, IntegratorException { 196dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 197dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond sanityChecks(equations, t0, y0, t, y); 198dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond setEquations(equations); 199dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond resetEvaluations(); 200dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond final boolean forward = t > t0; 201dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 202dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond // create some internal working arrays 203dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond final int stages = c.length + 1; 204dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond if (y != y0) { 205dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond System.arraycopy(y0, 0, y, 0, y0.length); 206dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond } 207dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond final double[][] yDotK = new double[stages][y0.length]; 208dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond final double[] yTmp = new double[y0.length]; 209dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond final double[] yDotTmp = new double[y0.length]; 210dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 211dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond // set up an interpolator sharing the integrator arrays 212dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond AbstractStepInterpolator interpolator; 213dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond if (requiresDenseOutput()) { 214dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond final RungeKuttaStepInterpolator rki = (RungeKuttaStepInterpolator) prototype.copy(); 215dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond rki.reinitialize(this, yTmp, yDotK, forward); 216dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond interpolator = rki; 217dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond } else { 218dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond interpolator = new DummyStepInterpolator(yTmp, yDotK[stages - 1], forward); 219dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond } 220dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond interpolator.storeTime(t0); 221dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 222dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond // set up integration control objects 223dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond stepStart = t0; 224dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond double hNew = 0; 225dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond boolean firstTime = true; 226dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond for (StepHandler handler : stepHandlers) { 227dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond handler.reset(); 228dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond } 229dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond setStateInitialized(false); 230dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 231dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond // main integration loop 232dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond isLastStep = false; 233dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond do { 234dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 235dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond interpolator.shift(); 236dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 237dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond // iterate over step size, ensuring local normalized error is smaller than 1 238dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond double error = 10; 239dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond while (error >= 1.0) { 240dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 241dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond if (firstTime || !fsal) { 242dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond // first stage 243dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond computeDerivatives(stepStart, y, yDotK[0]); 244dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond } 245dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 246dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond if (firstTime) { 247dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond final double[] scale = new double[mainSetDimension]; 248dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond if (vecAbsoluteTolerance == null) { 249dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond for (int i = 0; i < scale.length; ++i) { 250dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond scale[i] = scalAbsoluteTolerance + scalRelativeTolerance * FastMath.abs(y[i]); 251dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond } 252dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond } else { 253dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond for (int i = 0; i < scale.length; ++i) { 254dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond scale[i] = vecAbsoluteTolerance[i] + vecRelativeTolerance[i] * FastMath.abs(y[i]); 255dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond } 256dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond } 257dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond hNew = initializeStep(equations, forward, getOrder(), scale, 258dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond stepStart, y, yDotK[0], yTmp, yDotK[1]); 259dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond firstTime = false; 260dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond } 261dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 262dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond stepSize = hNew; 263dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 264dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond // next stages 265dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond for (int k = 1; k < stages; ++k) { 266dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 267dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond for (int j = 0; j < y0.length; ++j) { 268dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond double sum = a[k-1][0] * yDotK[0][j]; 269dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond for (int l = 1; l < k; ++l) { 270dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond sum += a[k-1][l] * yDotK[l][j]; 271dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond } 272dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond yTmp[j] = y[j] + stepSize * sum; 273dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond } 274dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 275dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond computeDerivatives(stepStart + c[k-1] * stepSize, yTmp, yDotK[k]); 276dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 277dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond } 278dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 279dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond // estimate the state at the end of the step 280dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond for (int j = 0; j < y0.length; ++j) { 281dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond double sum = b[0] * yDotK[0][j]; 282dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond for (int l = 1; l < stages; ++l) { 283dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond sum += b[l] * yDotK[l][j]; 284dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond } 285dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond yTmp[j] = y[j] + stepSize * sum; 286dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond } 287dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 288dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond // estimate the error at the end of the step 289dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond error = estimateError(yDotK, y, yTmp, stepSize); 290dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond if (error >= 1.0) { 291dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond // reject the step and attempt to reduce error by stepsize control 292dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond final double factor = 293dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond FastMath.min(maxGrowth, 294dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond FastMath.max(minReduction, safety * FastMath.pow(error, exp))); 295dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond hNew = filterStep(stepSize * factor, forward, false); 296dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond } 297dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 298dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond } 299dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 300dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond // local error is small enough: accept the step, trigger events and step handlers 301dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond interpolator.storeTime(stepStart + stepSize); 302dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond System.arraycopy(yTmp, 0, y, 0, y0.length); 303dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond System.arraycopy(yDotK[stages - 1], 0, yDotTmp, 0, y0.length); 304dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond stepStart = acceptStep(interpolator, y, yDotTmp, t); 305dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 306dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond if (!isLastStep) { 307dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 308dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond // prepare next step 309dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond interpolator.storeTime(stepStart); 310dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 311dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond if (fsal) { 312dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond // save the last evaluation for the next step 313dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond System.arraycopy(yDotTmp, 0, yDotK[0], 0, y0.length); 314dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond } 315dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 316dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond // stepsize control for next step 317dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond final double factor = 318dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond FastMath.min(maxGrowth, FastMath.max(minReduction, safety * FastMath.pow(error, exp))); 319dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond final double scaledH = stepSize * factor; 320dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond final double nextT = stepStart + scaledH; 321dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond final boolean nextIsLast = forward ? (nextT >= t) : (nextT <= t); 322dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond hNew = filterStep(scaledH, forward, nextIsLast); 323dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 324dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond final double filteredNextT = stepStart + hNew; 325dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond final boolean filteredNextIsLast = forward ? (filteredNextT >= t) : (filteredNextT <= t); 326dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond if (filteredNextIsLast) { 327dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond hNew = t - stepStart; 328dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond } 329dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 330dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond } 331dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 332dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond } while (!isLastStep); 333dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 334dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond final double stopTime = stepStart; 335dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond resetInternalState(); 336dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond return stopTime; 337dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 338dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond } 339dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 340dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond /** Get the minimal reduction factor for stepsize control. 341dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * @return minimal reduction factor 342dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond */ 343dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond public double getMinReduction() { 344dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond return minReduction; 345dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond } 346dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 347dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond /** Set the minimal reduction factor for stepsize control. 348dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * @param minReduction minimal reduction factor 349dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond */ 350dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond public void setMinReduction(final double minReduction) { 351dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond this.minReduction = minReduction; 352dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond } 353dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 354dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond /** Get the maximal growth factor for stepsize control. 355dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * @return maximal growth factor 356dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond */ 357dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond public double getMaxGrowth() { 358dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond return maxGrowth; 359dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond } 360dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 361dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond /** Set the maximal growth factor for stepsize control. 362dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * @param maxGrowth maximal growth factor 363dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond */ 364dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond public void setMaxGrowth(final double maxGrowth) { 365dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond this.maxGrowth = maxGrowth; 366dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond } 367dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 368dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond /** Compute the error ratio. 369dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * @param yDotK derivatives computed during the first stages 370dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * @param y0 estimate of the step at the start of the step 371dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * @param y1 estimate of the step at the end of the step 372dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * @param h current step 373dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * @return error ratio, greater than 1 if step should be rejected 374dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond */ 375dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond protected abstract double estimateError(double[][] yDotK, 376dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond double[] y0, double[] y1, 377dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond double h); 378dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 379dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond} 380