1 files changed, 269 insertions, 0 deletions
diff --git a/src/main/java/org/apache/commons/math3/ode/nonstiff/RungeKuttaIntegrator.java b/src/main/java/org/apache/commons/math3/ode/nonstiff/RungeKuttaIntegrator.java
new file mode 100644
index 0000000..5f7d5d8
--- /dev/null
+++ b/src/main/java/org/apache/commons/math3/ode/nonstiff/RungeKuttaIntegrator.java
@@ -0,0 +1,269 @@
+/*
+ * Licensed to the Apache Software Foundation (ASF) under one or more
+ * contributor license agreements.  See the NOTICE file distributed with
+ * this work for additional information regarding copyright ownership.
+ * The ASF licenses this file to You under the Apache License, Version 2.0
+ * (the "License"); you may not use this file except in compliance with
+ * the License.  You may obtain a copy of the License at
+ *
+ *      http://www.apache.org/licenses/LICENSE-2.0
+ *
+ * Unless required by applicable law or agreed to in writing, software
+ * distributed under the License is distributed on an "AS IS" BASIS,
+ * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+ * See the License for the specific language governing permissions and
+ * limitations under the License.
+ */
+
+package org.apache.commons.math3.ode.nonstiff;
+
+
+import org.apache.commons.math3.exception.DimensionMismatchException;
+import org.apache.commons.math3.exception.MaxCountExceededException;
+import org.apache.commons.math3.exception.NoBracketingException;
+import org.apache.commons.math3.exception.NumberIsTooSmallException;
+import org.apache.commons.math3.ode.AbstractIntegrator;
+import org.apache.commons.math3.ode.ExpandableStatefulODE;
+import org.apache.commons.math3.ode.FirstOrderDifferentialEquations;
+import org.apache.commons.math3.util.FastMath;
+
+/**
+ * This class implements the common part of all fixed step Runge-Kutta
+ * integrators for Ordinary Differential Equations.
+ *
+ * <p>These methods are explicit Runge-Kutta methods, their Butcher
+ * arrays are as follows :
+ * <pre>
+ *    0  |
+ *   c2  | a21
+ *   c3  | a31  a32
+ *   ... |        ...
+ *   cs  | as1  as2  ...  ass-1
+ *       |--------------------------
+ *       |  b1   b2  ...   bs-1  bs
+ * </pre>
+ * </p>
+ *
+ * @see EulerIntegrator
+ * @see ClassicalRungeKuttaIntegrator
+ * @see GillIntegrator
+ * @see MidpointIntegrator
+ * @since 1.2
+ */
+
+public abstract class RungeKuttaIntegrator extends AbstractIntegrator {
+
+    /** Time steps from Butcher array (without the first zero). */
+    private final double[] c;
+
+    /** Internal weights from Butcher array (without the first empty row). */
+    private final double[][] a;
+
+    /** External weights for the high order method from Butcher array. */
+    private final double[] b;
+
+    /** Prototype of the step interpolator. */
+    private final RungeKuttaStepInterpolator prototype;
+
+    /** Integration step. */
+    private final double step;
+
+  /** Simple constructor.
+   * Build a Runge-Kutta integrator with the given
+   * step. The default step handler does nothing.
+   * @param name name of the method
+   * @param c time steps from Butcher array (without the first zero)
+   * @param a internal weights from Butcher array (without the first empty row)
+   * @param b propagation weights for the high order method from Butcher array
+   * @param prototype prototype of the step interpolator to use
+   * @param step integration step
+   */
+  protected RungeKuttaIntegrator(final String name,
+                                 final double[] c, final double[][] a, final double[] b,
+                                 final RungeKuttaStepInterpolator prototype,
+                                 final double step) {
+    super(name);
+    this.c          = c;
+    this.a          = a;
+    this.b          = b;
+    this.prototype  = prototype;
+    this.step       = FastMath.abs(step);
+  }
+
+  /** {@inheritDoc} */
+  @Override
+  public void integrate(final ExpandableStatefulODE equations, final double t)
+      throws NumberIsTooSmallException, DimensionMismatchException,
+             MaxCountExceededException, NoBracketingException {
+
+    sanityChecks(equations, t);
+    setEquations(equations);
+    final boolean forward = t > equations.getTime();
+
+    // create some internal working arrays
+    final double[] y0      = equations.getCompleteState();
+    final double[] y       = y0.clone();
+    final int stages       = c.length + 1;
+    final double[][] yDotK = new double[stages][];
+    for (int i = 0; i < stages; ++i) {
+      yDotK [i] = new double[y0.length];
+    }
+    final double[] yTmp    = y0.clone();
+    final double[] yDotTmp = new double[y0.length];
+
+    // set up an interpolator sharing the integrator arrays
+    final RungeKuttaStepInterpolator interpolator = (RungeKuttaStepInterpolator) prototype.copy();
+    interpolator.reinitialize(this, yTmp, yDotK, forward,
+                              equations.getPrimaryMapper(), equations.getSecondaryMappers());
+    interpolator.storeTime(equations.getTime());
+
+    // set up integration control objects
+    stepStart = equations.getTime();
+    if (forward) {
+        if (stepStart + step >= t) {
+            stepSize = t - stepStart;
+        } else {
+            stepSize = step;
+        }
+    } else {
+        if (stepStart - step <= t) {
+            stepSize = t - stepStart;
+        } else {
+            stepSize = -step;
+        }
+    }
+    initIntegration(equations.getTime(), y0, t);
+
+    // main integration loop
+    isLastStep = false;
+    do {
+
+      interpolator.shift();
+
+      // first stage
+      computeDerivatives(stepStart, y, yDotK[0]);
+
+      // next stages
+      for (int k = 1; k < stages; ++k) {
+
+          for (int j = 0; j < y0.length; ++j) {
+              double sum = a[k-1][0] * yDotK[0][j];
+              for (int l = 1; l < k; ++l) {
+                  sum += a[k-1][l] * yDotK[l][j];
+              }
+              yTmp[j] = y[j] + stepSize * sum;
+          }
+
+          computeDerivatives(stepStart + c[k-1] * stepSize, yTmp, yDotK[k]);
+
+      }
+
+      // estimate the state at the end of the step
+      for (int j = 0; j < y0.length; ++j) {
+          double sum    = b[0] * yDotK[0][j];
+          for (int l = 1; l < stages; ++l) {
+              sum    += b[l] * yDotK[l][j];
+          }
+          yTmp[j] = y[j] + stepSize * sum;
+      }
+
+      // discrete events handling
+      interpolator.storeTime(stepStart + stepSize);
+      System.arraycopy(yTmp, 0, y, 0, y0.length);
+      System.arraycopy(yDotK[stages - 1], 0, yDotTmp, 0, y0.length);
+      stepStart = acceptStep(interpolator, y, yDotTmp, t);
+
+      if (!isLastStep) {
+
+          // prepare next step
+          interpolator.storeTime(stepStart);
+
+          // stepsize control for next step
+          final double  nextT      = stepStart + stepSize;
+          final boolean nextIsLast = forward ? (nextT >= t) : (nextT <= t);
+          if (nextIsLast) {
+              stepSize = t - stepStart;
+          }
+      }
+
+    } while (!isLastStep);
+
+    // dispatch results
+    equations.setTime(stepStart);
+    equations.setCompleteState(y);
+
+    stepStart = Double.NaN;
+    stepSize  = Double.NaN;
+
+  }
+
+  /** Fast computation of a single step of ODE integration.
+   * <p>This method is intended for the limited use case of
+   * very fast computation of only one step without using any of the
+   * rich features of general integrators that may take some time
+   * to set up (i.e. no step handlers, no events handlers, no additional
+   * states, no interpolators, no error control, no evaluations count,
+   * no sanity checks ...). It handles the strict minimum of computation,
+   * so it can be embedded in outer loops.</p>
+   * <p>
+   * This method is <em>not</em> used at all by the {@link #integrate(ExpandableStatefulODE, double)}
+   * method. It also completely ignores the step set at construction time, and
+   * uses only a single step to go from {@code t0} to {@code t}.
+   * </p>
+   * <p>
+   * As this method does not use any of the state-dependent features of the integrator,
+   * it should be reasonably thread-safe <em>if and only if</em> the provided differential
+   * equations are themselves thread-safe.
+   * </p>
+   * @param equations differential equations to integrate
+   * @param t0 initial time
+   * @param y0 initial value of the state vector at t0
+   * @param t target time for the integration
+   * (can be set to a value smaller than {@code t0} for backward integration)
+   * @return state vector at {@code t}
+   */
+  public double[] singleStep(final FirstOrderDifferentialEquations equations,
+                             final double t0, final double[] y0, final double t) {
+
+      // create some internal working arrays
+      final double[] y       = y0.clone();
+      final int stages       = c.length + 1;
+      final double[][] yDotK = new double[stages][];
+      for (int i = 0; i < stages; ++i) {
+          yDotK [i] = new double[y0.length];
+      }
+      final double[] yTmp    = y0.clone();
+
+      // first stage
+      final double h = t - t0;
+      equations.computeDerivatives(t0, y, yDotK[0]);
+
+      // next stages
+      for (int k = 1; k < stages; ++k) {
+
+          for (int j = 0; j < y0.length; ++j) {
+              double sum = a[k-1][0] * yDotK[0][j];
+              for (int l = 1; l < k; ++l) {
+                  sum += a[k-1][l] * yDotK[l][j];
+              }
+              yTmp[j] = y[j] + h * sum;
+          }
+
+          equations.computeDerivatives(t0 + c[k-1] * h, yTmp, yDotK[k]);
+
+      }
+
+      // estimate the state at the end of the step
+      for (int j = 0; j < y0.length; ++j) {
+          double sum = b[0] * yDotK[0][j];
+          for (int l = 1; l < stages; ++l) {
+              sum += b[l] * yDotK[l][j];
+          }
+          y[j] += h * sum;
+      }
+
+      return y;
+
+  }
+
+}