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 }