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 20dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 21dee0849a9704d532af0b550146cbafbaa6ee1d19Raymondimport org.apache.commons.math.ode.AbstractIntegrator; 22dee0849a9704d532af0b550146cbafbaa6ee1d19Raymondimport org.apache.commons.math.ode.DerivativeException; 23dee0849a9704d532af0b550146cbafbaa6ee1d19Raymondimport org.apache.commons.math.ode.FirstOrderDifferentialEquations; 24dee0849a9704d532af0b550146cbafbaa6ee1d19Raymondimport org.apache.commons.math.ode.IntegratorException; 25dee0849a9704d532af0b550146cbafbaa6ee1d19Raymondimport org.apache.commons.math.ode.sampling.AbstractStepInterpolator; 26dee0849a9704d532af0b550146cbafbaa6ee1d19Raymondimport org.apache.commons.math.ode.sampling.DummyStepInterpolator; 27dee0849a9704d532af0b550146cbafbaa6ee1d19Raymondimport org.apache.commons.math.ode.sampling.StepHandler; 28dee0849a9704d532af0b550146cbafbaa6ee1d19Raymondimport org.apache.commons.math.util.FastMath; 29dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 30dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond/** 31dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * This class implements the common part of all fixed step Runge-Kutta 32dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * integrators for Ordinary Differential Equations. 33dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * 34dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * <p>These methods are explicit Runge-Kutta methods, their Butcher 35dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * arrays are as follows : 36dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * <pre> 37dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * 0 | 38dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * c2 | a21 39dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * c3 | a31 a32 40dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * ... | ... 41dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * cs | as1 as2 ... ass-1 42dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * |-------------------------- 43dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * | b1 b2 ... bs-1 bs 44dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * </pre> 45dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * </p> 46dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * 47dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * @see EulerIntegrator 48dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * @see ClassicalRungeKuttaIntegrator 49dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * @see GillIntegrator 50dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * @see MidpointIntegrator 51dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * @version $Revision: 1073158 $ $Date: 2011-02-21 22:46:52 +0100 (lun. 21 févr. 2011) $ 52dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * @since 1.2 53dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond */ 54dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 55dee0849a9704d532af0b550146cbafbaa6ee1d19Raymondpublic abstract class RungeKuttaIntegrator extends AbstractIntegrator { 56dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 57dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond /** Time steps from Butcher array (without the first zero). */ 58dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond private final double[] c; 59dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 60dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond /** Internal weights from Butcher array (without the first empty row). */ 61dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond private final double[][] a; 62dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 63dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond /** External weights for the high order method from Butcher array. */ 64dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond private final double[] b; 65dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 66dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond /** Prototype of the step interpolator. */ 67dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond private final RungeKuttaStepInterpolator prototype; 68dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 69dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond /** Integration step. */ 70dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond private final double step; 71dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 72dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond /** Simple constructor. 73dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * Build a Runge-Kutta integrator with the given 74dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * step. The default step handler does nothing. 75dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * @param name name of the method 76dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * @param c time steps from Butcher array (without the first zero) 77dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * @param a internal weights from Butcher array (without the first empty row) 78dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * @param b propagation weights for the high order method from Butcher array 79dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * @param prototype prototype of the step interpolator to use 80dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * @param step integration step 81dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond */ 82dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond protected RungeKuttaIntegrator(final String name, 83dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond final double[] c, final double[][] a, final double[] b, 84dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond final RungeKuttaStepInterpolator prototype, 85dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond final double step) { 86dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond super(name); 87dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond this.c = c; 88dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond this.a = a; 89dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond this.b = b; 90dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond this.prototype = prototype; 91dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond this.step = FastMath.abs(step); 92dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond } 93dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 94dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond /** {@inheritDoc} */ 95dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond public double integrate(final FirstOrderDifferentialEquations equations, 96dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond final double t0, final double[] y0, 97dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond final double t, final double[] y) 98dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond throws DerivativeException, IntegratorException { 99dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 100dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond sanityChecks(equations, t0, y0, t, y); 101dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond setEquations(equations); 102dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond resetEvaluations(); 103dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond final boolean forward = t > t0; 104dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 105dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond // create some internal working arrays 106dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond final int stages = c.length + 1; 107dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond if (y != y0) { 108dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond System.arraycopy(y0, 0, y, 0, y0.length); 109dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond } 110dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond final double[][] yDotK = new double[stages][]; 111dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond for (int i = 0; i < stages; ++i) { 112dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond yDotK [i] = new double[y0.length]; 113dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond } 114dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond final double[] yTmp = new double[y0.length]; 115dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond final double[] yDotTmp = new double[y0.length]; 116dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 117dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond // set up an interpolator sharing the integrator arrays 118dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond AbstractStepInterpolator interpolator; 119dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond if (requiresDenseOutput()) { 120dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond final RungeKuttaStepInterpolator rki = (RungeKuttaStepInterpolator) prototype.copy(); 121dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond rki.reinitialize(this, yTmp, yDotK, forward); 122dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond interpolator = rki; 123dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond } else { 124dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond interpolator = new DummyStepInterpolator(yTmp, yDotK[stages - 1], forward); 125dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond } 126dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond interpolator.storeTime(t0); 127dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 128dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond // set up integration control objects 129dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond stepStart = t0; 130dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond stepSize = forward ? step : -step; 131dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond for (StepHandler handler : stepHandlers) { 132dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond handler.reset(); 133dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond } 134dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond setStateInitialized(false); 135dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 136dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond // main integration loop 137dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond isLastStep = false; 138dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond do { 139dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 140dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond interpolator.shift(); 141dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 142dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond // first stage 143dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond computeDerivatives(stepStart, y, yDotK[0]); 144dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 145dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond // next stages 146dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond for (int k = 1; k < stages; ++k) { 147dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 148dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond for (int j = 0; j < y0.length; ++j) { 149dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond double sum = a[k-1][0] * yDotK[0][j]; 150dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond for (int l = 1; l < k; ++l) { 151dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond sum += a[k-1][l] * yDotK[l][j]; 152dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond } 153dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond yTmp[j] = y[j] + stepSize * sum; 154dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond } 155dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 156dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond computeDerivatives(stepStart + c[k-1] * stepSize, yTmp, yDotK[k]); 157dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 158dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond } 159dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 160dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond // estimate the state at the end of the step 161dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond for (int j = 0; j < y0.length; ++j) { 162dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond double sum = b[0] * yDotK[0][j]; 163dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond for (int l = 1; l < stages; ++l) { 164dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond sum += b[l] * yDotK[l][j]; 165dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond } 166dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond yTmp[j] = y[j] + stepSize * sum; 167dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond } 168dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 169dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond // discrete events handling 170dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond interpolator.storeTime(stepStart + stepSize); 171dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond System.arraycopy(yTmp, 0, y, 0, y0.length); 172dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond System.arraycopy(yDotK[stages - 1], 0, yDotTmp, 0, y0.length); 173dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond stepStart = acceptStep(interpolator, y, yDotTmp, t); 174dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 175dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond if (!isLastStep) { 176dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 177dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond // prepare next step 178dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond interpolator.storeTime(stepStart); 179dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 180dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond // stepsize control for next step 181dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond final double nextT = stepStart + stepSize; 182dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond final boolean nextIsLast = forward ? (nextT >= t) : (nextT <= t); 183dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond if (nextIsLast) { 184dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond stepSize = t - stepStart; 185dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond } 186dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond } 187dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 188dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond } while (!isLastStep); 189dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 190dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond final double stopTime = stepStart; 191dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond stepStart = Double.NaN; 192dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond stepSize = Double.NaN; 193dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond return stopTime; 194dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 195dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond } 196dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 197dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond} 198