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.math3.ode.nonstiff;
019
020
021 import org.apache.commons.math3.exception.DimensionMismatchException;
022 import org.apache.commons.math3.exception.MaxCountExceededException;
023 import org.apache.commons.math3.exception.NoBracketingException;
024 import org.apache.commons.math3.exception.NumberIsTooSmallException;
025 import org.apache.commons.math3.ode.AbstractIntegrator;
026 import org.apache.commons.math3.ode.ExpandableStatefulODE;
027 import org.apache.commons.math3.util.FastMath;
028
029 /**
030 * This class implements the common part of all fixed step Runge-Kutta
031 * integrators for Ordinary Differential Equations.
032 *
033 * <p>These methods are explicit Runge-Kutta methods, their Butcher
034 * arrays are as follows :
035 * <pre>
036 * 0 |
037 * c2 | a21
038 * c3 | a31 a32
039 * ... | ...
040 * cs | as1 as2 ... ass-1
041 * |--------------------------
042 * | b1 b2 ... bs-1 bs
043 * </pre>
044 * </p>
045 *
046 * @see EulerIntegrator
047 * @see ClassicalRungeKuttaIntegrator
048 * @see GillIntegrator
049 * @see MidpointIntegrator
050 * @version $Id: RungeKuttaIntegrator.java 1416643 2012-12-03 19:37:14Z tn $
051 * @since 1.2
052 */
053
054 public abstract class RungeKuttaIntegrator extends AbstractIntegrator {
055
056 /** Time steps from Butcher array (without the first zero). */
057 private final double[] c;
058
059 /** Internal weights from Butcher array (without the first empty row). */
060 private final double[][] a;
061
062 /** External weights for the high order method from Butcher array. */
063 private final double[] b;
064
065 /** Prototype of the step interpolator. */
066 private final RungeKuttaStepInterpolator prototype;
067
068 /** Integration step. */
069 private final double step;
070
071 /** Simple constructor.
072 * Build a Runge-Kutta integrator with the given
073 * step. The default step handler does nothing.
074 * @param name name of the method
075 * @param c time steps from Butcher array (without the first zero)
076 * @param a internal weights from Butcher array (without the first empty row)
077 * @param b propagation weights for the high order method from Butcher array
078 * @param prototype prototype of the step interpolator to use
079 * @param step integration step
080 */
081 protected RungeKuttaIntegrator(final String name,
082 final double[] c, final double[][] a, final double[] b,
083 final RungeKuttaStepInterpolator prototype,
084 final double step) {
085 super(name);
086 this.c = c;
087 this.a = a;
088 this.b = b;
089 this.prototype = prototype;
090 this.step = FastMath.abs(step);
091 }
092
093 /** {@inheritDoc} */
094 @Override
095 public void integrate(final ExpandableStatefulODE equations, final double t)
096 throws NumberIsTooSmallException, DimensionMismatchException,
097 MaxCountExceededException, NoBracketingException {
098
099 sanityChecks(equations, t);
100 setEquations(equations);
101 final boolean forward = t > equations.getTime();
102
103 // create some internal working arrays
104 final double[] y0 = equations.getCompleteState();
105 final double[] y = y0.clone();
106 final int stages = c.length + 1;
107 final double[][] yDotK = new double[stages][];
108 for (int i = 0; i < stages; ++i) {
109 yDotK [i] = new double[y0.length];
110 }
111 final double[] yTmp = y0.clone();
112 final double[] yDotTmp = new double[y0.length];
113
114 // set up an interpolator sharing the integrator arrays
115 final RungeKuttaStepInterpolator interpolator = (RungeKuttaStepInterpolator) prototype.copy();
116 interpolator.reinitialize(this, yTmp, yDotK, forward,
117 equations.getPrimaryMapper(), equations.getSecondaryMappers());
118 interpolator.storeTime(equations.getTime());
119
120 // set up integration control objects
121 stepStart = equations.getTime();
122 stepSize = forward ? step : -step;
123 initIntegration(equations.getTime(), y0, t);
124
125 // main integration loop
126 isLastStep = false;
127 do {
128
129 interpolator.shift();
130
131 // first stage
132 computeDerivatives(stepStart, y, yDotK[0]);
133
134 // next stages
135 for (int k = 1; k < stages; ++k) {
136
137 for (int j = 0; j < y0.length; ++j) {
138 double sum = a[k-1][0] * yDotK[0][j];
139 for (int l = 1; l < k; ++l) {
140 sum += a[k-1][l] * yDotK[l][j];
141 }
142 yTmp[j] = y[j] + stepSize * sum;
143 }
144
145 computeDerivatives(stepStart + c[k-1] * stepSize, yTmp, yDotK[k]);
146
147 }
148
149 // estimate the state at the end of the step
150 for (int j = 0; j < y0.length; ++j) {
151 double sum = b[0] * yDotK[0][j];
152 for (int l = 1; l < stages; ++l) {
153 sum += b[l] * yDotK[l][j];
154 }
155 yTmp[j] = y[j] + stepSize * sum;
156 }
157
158 // discrete events handling
159 interpolator.storeTime(stepStart + stepSize);
160 System.arraycopy(yTmp, 0, y, 0, y0.length);
161 System.arraycopy(yDotK[stages - 1], 0, yDotTmp, 0, y0.length);
162 stepStart = acceptStep(interpolator, y, yDotTmp, t);
163
164 if (!isLastStep) {
165
166 // prepare next step
167 interpolator.storeTime(stepStart);
168
169 // stepsize control for next step
170 final double nextT = stepStart + stepSize;
171 final boolean nextIsLast = forward ? (nextT >= t) : (nextT <= t);
172 if (nextIsLast) {
173 stepSize = t - stepStart;
174 }
175 }
176
177 } while (!isLastStep);
178
179 // dispatch results
180 equations.setTime(stepStart);
181 equations.setCompleteState(y);
182
183 stepStart = Double.NaN;
184 stepSize = Double.NaN;
185
186 }
187
188 }