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