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 018 package org.apache.commons.math.ode.nonstiff; 019 020 021 import org.apache.commons.math.ode.AbstractIntegrator; 022 import org.apache.commons.math.ode.DerivativeException; 023 import org.apache.commons.math.ode.FirstOrderDifferentialEquations; 024 import org.apache.commons.math.ode.IntegratorException; 025 import org.apache.commons.math.ode.events.CombinedEventsManager; 026 import org.apache.commons.math.ode.sampling.AbstractStepInterpolator; 027 import org.apache.commons.math.ode.sampling.DummyStepInterpolator; 028 import org.apache.commons.math.ode.sampling.StepHandler; 029 030 /** 031 * This class implements the common part of all fixed step Runge-Kutta 032 * integrators for Ordinary Differential Equations. 033 * 034 * <p>These methods are explicit Runge-Kutta methods, their Butcher 035 * arrays are as follows : 036 * <pre> 037 * 0 | 038 * c2 | a21 039 * c3 | a31 a32 040 * ... | ... 041 * cs | as1 as2 ... ass-1 042 * |-------------------------- 043 * | b1 b2 ... bs-1 bs 044 * </pre> 045 * </p> 046 * 047 * @see EulerIntegrator 048 * @see ClassicalRungeKuttaIntegrator 049 * @see GillIntegrator 050 * @see MidpointIntegrator 051 * @version $Revision: 927202 $ $Date: 2010-03-24 18:11:51 -0400 (Wed, 24 Mar 2010) $ 052 * @since 1.2 053 */ 054 055 public abstract class RungeKuttaIntegrator extends AbstractIntegrator { 056 057 /** Time steps from Butcher array (without the first zero). */ 058 private final double[] c; 059 060 /** Internal weights from Butcher array (without the first empty row). */ 061 private final double[][] a; 062 063 /** External weights for the high order method from Butcher array. */ 064 private final double[] b; 065 066 /** Prototype of the step interpolator. */ 067 private final RungeKuttaStepInterpolator prototype; 068 069 /** Integration step. */ 070 private final double step; 071 072 /** Simple constructor. 073 * Build a Runge-Kutta integrator with the given 074 * step. The default step handler does nothing. 075 * @param name name of the method 076 * @param c time steps from Butcher array (without the first zero) 077 * @param a internal weights from Butcher array (without the first empty row) 078 * @param b propagation weights for the high order method from Butcher array 079 * @param prototype prototype of the step interpolator to use 080 * @param step integration step 081 */ 082 protected RungeKuttaIntegrator(final String name, 083 final double[] c, final double[][] a, final double[] b, 084 final RungeKuttaStepInterpolator prototype, 085 final double step) { 086 super(name); 087 this.c = c; 088 this.a = a; 089 this.b = b; 090 this.prototype = prototype; 091 this.step = Math.abs(step); 092 } 093 094 /** {@inheritDoc} */ 095 public double integrate(final FirstOrderDifferentialEquations equations, 096 final double t0, final double[] y0, 097 final double t, final double[] y) 098 throws DerivativeException, IntegratorException { 099 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 116 // set up an interpolator sharing the integrator arrays 117 AbstractStepInterpolator interpolator; 118 if (requiresDenseOutput() || (! eventsHandlersManager.isEmpty())) { 119 final RungeKuttaStepInterpolator rki = (RungeKuttaStepInterpolator) prototype.copy(); 120 rki.reinitialize(this, yTmp, yDotK, forward); 121 interpolator = rki; 122 } else { 123 interpolator = new DummyStepInterpolator(yTmp, yDotK[stages - 1], forward); 124 } 125 interpolator.storeTime(t0); 126 127 // set up integration control objects 128 stepStart = t0; 129 stepSize = forward ? step : -step; 130 for (StepHandler handler : stepHandlers) { 131 handler.reset(); 132 } 133 CombinedEventsManager manager = addEndTimeChecker(t0, t, eventsHandlersManager); 134 boolean lastStep = false; 135 136 // main integration loop 137 while (!lastStep) { 138 139 interpolator.shift(); 140 141 for (boolean loop = true; loop;) { 142 143 // first stage 144 computeDerivatives(stepStart, y, yDotK[0]); 145 146 // next stages 147 for (int k = 1; k < stages; ++k) { 148 149 for (int j = 0; j < y0.length; ++j) { 150 double sum = a[k-1][0] * yDotK[0][j]; 151 for (int l = 1; l < k; ++l) { 152 sum += a[k-1][l] * yDotK[l][j]; 153 } 154 yTmp[j] = y[j] + stepSize * sum; 155 } 156 157 computeDerivatives(stepStart + c[k-1] * stepSize, yTmp, yDotK[k]); 158 159 } 160 161 // estimate the state at the end of the step 162 for (int j = 0; j < y0.length; ++j) { 163 double sum = b[0] * yDotK[0][j]; 164 for (int l = 1; l < stages; ++l) { 165 sum += b[l] * yDotK[l][j]; 166 } 167 yTmp[j] = y[j] + stepSize * sum; 168 } 169 170 // discrete events handling 171 interpolator.storeTime(stepStart + stepSize); 172 if (manager.evaluateStep(interpolator)) { 173 final double dt = manager.getEventTime() - stepStart; 174 if (Math.abs(dt) <= Math.ulp(stepStart)) { 175 // we cannot simply truncate the step, reject the current computation 176 // and let the loop compute another state with the truncated step. 177 // it is so small (much probably exactly 0 due to limited accuracy) 178 // that the code above would fail handling it. 179 // So we set up an artificial 0 size step by copying states 180 interpolator.storeTime(stepStart); 181 System.arraycopy(y, 0, yTmp, 0, y0.length); 182 stepSize = 0; 183 loop = false; 184 } else { 185 // reject the step to match exactly the next switch time 186 stepSize = dt; 187 } 188 } else { 189 loop = false; 190 } 191 192 } 193 194 // the step has been accepted 195 final double nextStep = stepStart + stepSize; 196 System.arraycopy(yTmp, 0, y, 0, y0.length); 197 manager.stepAccepted(nextStep, y); 198 lastStep = manager.stop(); 199 200 // provide the step data to the step handler 201 interpolator.storeTime(nextStep); 202 for (StepHandler handler : stepHandlers) { 203 handler.handleStep(interpolator, lastStep); 204 } 205 stepStart = nextStep; 206 207 if (manager.reset(stepStart, y) && ! lastStep) { 208 // some events handler has triggered changes that 209 // invalidate the derivatives, we need to recompute them 210 computeDerivatives(stepStart, y, yDotK[0]); 211 } 212 213 // make sure step size is set to default before next step 214 stepSize = forward ? step : -step; 215 216 } 217 218 final double stopTime = stepStart; 219 stepStart = Double.NaN; 220 stepSize = Double.NaN; 221 return stopTime; 222 223 } 224 225 }