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: 785473 $ $Date: 2009-06-17 00:02:35 -0400 (Wed, 17 Jun 2009) $ 052 * @since 1.2 053 */ 054 055 public abstract class RungeKuttaIntegrator extends AbstractIntegrator { 056 057 /** Simple constructor. 058 * Build a Runge-Kutta integrator with the given 059 * step. The default step handler does nothing. 060 * @param name name of the method 061 * @param c time steps from Butcher array (without the first zero) 062 * @param a internal weights from Butcher array (without the first empty row) 063 * @param b propagation weights for the high order method from Butcher array 064 * @param prototype prototype of the step interpolator to use 065 * @param step integration step 066 */ 067 protected RungeKuttaIntegrator(final String name, 068 final double[] c, final double[][] a, final double[] b, 069 final RungeKuttaStepInterpolator prototype, 070 final double step) { 071 super(name); 072 this.c = c; 073 this.a = a; 074 this.b = b; 075 this.prototype = prototype; 076 this.step = Math.abs(step); 077 } 078 079 /** {@inheritDoc} */ 080 public double integrate(final FirstOrderDifferentialEquations equations, 081 final double t0, final double[] y0, 082 final double t, final double[] y) 083 throws DerivativeException, IntegratorException { 084 085 sanityChecks(equations, t0, y0, t, y); 086 setEquations(equations); 087 resetEvaluations(); 088 final boolean forward = (t > t0); 089 090 // create some internal working arrays 091 final int stages = c.length + 1; 092 if (y != y0) { 093 System.arraycopy(y0, 0, y, 0, y0.length); 094 } 095 final double[][] yDotK = new double[stages][]; 096 for (int i = 0; i < stages; ++i) { 097 yDotK [i] = new double[y0.length]; 098 } 099 final double[] yTmp = new double[y0.length]; 100 101 // set up an interpolator sharing the integrator arrays 102 AbstractStepInterpolator interpolator; 103 if (requiresDenseOutput() || (! eventsHandlersManager.isEmpty())) { 104 final RungeKuttaStepInterpolator rki = (RungeKuttaStepInterpolator) prototype.copy(); 105 rki.reinitialize(this, yTmp, yDotK, forward); 106 interpolator = rki; 107 } else { 108 interpolator = new DummyStepInterpolator(yTmp, forward); 109 } 110 interpolator.storeTime(t0); 111 112 // set up integration control objects 113 stepStart = t0; 114 stepSize = forward ? step : -step; 115 for (StepHandler handler : stepHandlers) { 116 handler.reset(); 117 } 118 CombinedEventsManager manager = addEndTimeChecker(t0, t, eventsHandlersManager); 119 boolean lastStep = false; 120 121 // main integration loop 122 while (!lastStep) { 123 124 interpolator.shift(); 125 126 for (boolean loop = true; loop;) { 127 128 // first stage 129 computeDerivatives(stepStart, y, yDotK[0]); 130 131 // next stages 132 for (int k = 1; k < stages; ++k) { 133 134 for (int j = 0; j < y0.length; ++j) { 135 double sum = a[k-1][0] * yDotK[0][j]; 136 for (int l = 1; l < k; ++l) { 137 sum += a[k-1][l] * yDotK[l][j]; 138 } 139 yTmp[j] = y[j] + stepSize * sum; 140 } 141 142 computeDerivatives(stepStart + c[k-1] * stepSize, yTmp, yDotK[k]); 143 144 } 145 146 // estimate the state at the end of the step 147 for (int j = 0; j < y0.length; ++j) { 148 double sum = b[0] * yDotK[0][j]; 149 for (int l = 1; l < stages; ++l) { 150 sum += b[l] * yDotK[l][j]; 151 } 152 yTmp[j] = y[j] + stepSize * sum; 153 } 154 155 // discrete events handling 156 interpolator.storeTime(stepStart + stepSize); 157 if (manager.evaluateStep(interpolator)) { 158 final double dt = manager.getEventTime() - stepStart; 159 if (Math.abs(dt) <= Math.ulp(stepStart)) { 160 // rejecting the step would lead to a too small next step, we accept it 161 loop = false; 162 } else { 163 // reject the step to match exactly the next switch time 164 stepSize = dt; 165 } 166 } else { 167 loop = false; 168 } 169 170 } 171 172 // the step has been accepted 173 final double nextStep = stepStart + stepSize; 174 System.arraycopy(yTmp, 0, y, 0, y0.length); 175 manager.stepAccepted(nextStep, y); 176 lastStep = manager.stop(); 177 178 // provide the step data to the step handler 179 interpolator.storeTime(nextStep); 180 for (StepHandler handler : stepHandlers) { 181 handler.handleStep(interpolator, lastStep); 182 } 183 stepStart = nextStep; 184 185 if (manager.reset(stepStart, y) && ! lastStep) { 186 // some events handler has triggered changes that 187 // invalidate the derivatives, we need to recompute them 188 computeDerivatives(stepStart, y, yDotK[0]); 189 } 190 191 // make sure step size is set to default before next step 192 stepSize = forward ? step : -step; 193 194 } 195 196 final double stopTime = stepStart; 197 stepStart = Double.NaN; 198 stepSize = Double.NaN; 199 return stopTime; 200 201 } 202 203 /** Time steps from Butcher array (without the first zero). */ 204 private double[] c; 205 206 /** Internal weights from Butcher array (without the first empty row). */ 207 private double[][] a; 208 209 /** External weights for the high order method from Butcher array. */ 210 private double[] b; 211 212 /** Prototype of the step interpolator. */ 213 private RungeKuttaStepInterpolator prototype; 214 215 /** Integration step. */ 216 private double step; 217 218 }