commons-dev mailing list archives

Site index · List index
Message view « Date » · « Thread »
Top « Date » · « Thread »
From l..@apache.org
Subject svn commit: r512066 [3/5] - in /jakarta/commons/proper/math/trunk/src: java/org/apache/commons/math/ode/ mantissa/src/org/spaceroots/mantissa/ode/ mantissa/src/org/spaceroots/mantissa/ode/doc-files/ mantissa/tests-src/org/spaceroots/mantissa/ode/ test/...
Date Mon, 26 Feb 2007 23:12:45 GMT
Added: jakarta/commons/proper/math/trunk/src/java/org/apache/commons/math/ode/GraggBulirschStoerIntegrator.java
URL: http://svn.apache.org/viewvc/jakarta/commons/proper/math/trunk/src/java/org/apache/commons/math/ode/GraggBulirschStoerIntegrator.java?view=auto&rev=512066
==============================================================================
--- jakarta/commons/proper/math/trunk/src/java/org/apache/commons/math/ode/GraggBulirschStoerIntegrator.java (added)
+++ jakarta/commons/proper/math/trunk/src/java/org/apache/commons/math/ode/GraggBulirschStoerIntegrator.java Mon Feb 26 15:12:40 2007
@@ -0,0 +1,1020 @@
+/*
+ * 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.math.ode;
+
+/**
+ * This class implements a Gragg-Bulirsch-Stoer integrator for
+ * Ordinary Differential Equations.
+
+ * <p>The Gragg-Bulirsch-Stoer algorithm is one of the most efficient
+ * ones currently available for smooth problems. It uses Richardson
+ * extrapolation to estimate what would be the solution if the step
+ * size could be decreased down to zero.</p>
+
+ * <p>
+ * This method changes both the step size and the order during
+ * integration, in order to minimize computation cost. It is
+ * particularly well suited when a very high precision is needed. The
+ * limit where this method becomes more efficient than high-order
+ * Runge-Kutta-Fehlberg methods like {@link DormandPrince853Integrator
+ * Dormand-Prince 8(5,3)} depends on the problem. Results given in the
+ * Hairer, Norsett and Wanner book show for example that this limit
+ * occurs for accuracy around 1e-6 when integrating Saltzam-Lorenz
+ * equations (the authors note this problem is <i>extremely sensitive
+ * to the errors in the first integration steps</i>), and around 1e-11
+ * for a two dimensional celestial mechanics problems with seven
+ * bodies (pleiades problem, involving quasi-collisions for which
+ * <i>automatic step size control is essential</i>).
+ * </p>
+
+ * <p>
+ * This implementation is basically a reimplementation in Java of the
+ * <a
+ * href="http://www.unige.ch/math/folks/hairer/prog/nonstiff/odex.f">odex</a>
+ * fortran code by E. Hairer and G. Wanner. The redistribution policy
+ * for this code is available <a
+ * href="http://www.unige.ch/~hairer/prog/licence.txt">here</a>, for
+ * convenience, it is reproduced below.</p>
+ * </p>
+
+ * <table border="0" width="80%" cellpadding="10" align="center" bgcolor="#E0E0E0">
+ * <tr><td>Copyright (c) 2004, Ernst Hairer</td></tr>
+
+ * <tr><td>Redistribution and use in source and binary forms, with or
+ * without modification, are permitted provided that the following
+ * conditions are met:
+ * <ul>
+ *  <li>Redistributions of source code must retain the above copyright
+ *      notice, this list of conditions and the following disclaimer.</li>
+ *  <li>Redistributions in binary form must reproduce the above copyright
+ *      notice, this list of conditions and the following disclaimer in the
+ *      documentation and/or other materials provided with the distribution.</li>
+ * </ul></td></tr>
+
+ * <tr><td><strong>THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND
+ * CONTRIBUTORS "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING,
+ * BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS
+ * FOR A  PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE REGENTS OR
+ * CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
+ * EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
+ * PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
+ * PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
+ * LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
+ * NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
+ * SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.</strong></td></tr>
+ * </table>
+
+ * @author E. Hairer and G. Wanner (fortran version)
+ * @author L. Maisonobe (Java port)
+ * @version $Id: GraggBulirschStoerIntegrator.java 1705 2006-09-17 19:57:39Z luc $
+
+ */
+
+public class GraggBulirschStoerIntegrator
+  extends AdaptiveStepsizeIntegrator {
+
+  private static final String methodName = "Gragg-Bulirsch-Stoer";
+
+  /** Simple constructor.
+   * Build a Gragg-Bulirsch-Stoer integrator with the given step
+   * bounds. All tuning parameters are set to their default
+   * values. The default step handler does nothing.
+   * @param minStep minimal step (must be positive even for backward
+   * integration), the last step can be smaller than this
+   * @param maxStep maximal step (must be positive even for backward
+   * integration)
+   * @param scalAbsoluteTolerance allowed absolute error
+   * @param scalRelativeTolerance allowed relative error
+   */
+  public GraggBulirschStoerIntegrator(double minStep, double maxStep,
+                                      double scalAbsoluteTolerance,
+                                      double scalRelativeTolerance) {
+    super(minStep, maxStep, scalAbsoluteTolerance, scalRelativeTolerance);
+    denseOutput = (handler.requiresDenseOutput()
+                   || (! switchesHandler.isEmpty()));
+    setStabilityCheck(true, -1, -1, -1);
+    setStepsizeControl(-1, -1, -1, -1);
+    setOrderControl(-1, -1, -1);
+    setInterpolationControl(true, -1);
+  }
+
+  /** Simple constructor.
+   * Build a Gragg-Bulirsch-Stoer integrator with the given step
+   * bounds. All tuning parameters are set to their default
+   * values. The default step handler does nothing.
+   * @param minStep minimal step (must be positive even for backward
+   * integration), the last step can be smaller than this
+   * @param maxStep maximal step (must be positive even for backward
+   * integration)
+   * @param vecAbsoluteTolerance allowed absolute error
+   * @param vecRelativeTolerance allowed relative error
+   */
+  public GraggBulirschStoerIntegrator(double minStep, double maxStep,
+                                      double[] vecAbsoluteTolerance,
+                                      double[] vecRelativeTolerance) {
+    super(minStep, maxStep, vecAbsoluteTolerance, vecRelativeTolerance);
+    denseOutput = (handler.requiresDenseOutput()
+                   || (! switchesHandler.isEmpty()));
+    setStabilityCheck(true, -1, -1, -1);
+    setStepsizeControl(-1, -1, -1, -1);
+    setOrderControl(-1, -1, -1);
+    setInterpolationControl(true, -1);
+  }
+
+  /** Set the stability check controls.
+   * <p>The stability check is performed on the first few iterations of
+   * the extrapolation scheme. If this test fails, the step is rejected
+   * and the stepsize is reduced.</p>
+   * <p>By default, the test is performed, at most during two
+   * iterations at each step, and at most once for each of these
+   * iterations. The default stepsize reduction factor is 0.5.</p>
+   * @param performTest if true, stability check will be performed,
+     if false, the check will be skipped
+   * @param maxIter maximal number of iterations for which checks are
+   * performed (the number of iterations is reset to default if negative
+   * or null)
+   * @param maxChecks maximal number of checks for each iteration
+   * (the number of checks is reset to default if negative or null)
+   * @param stabilityReduction stepsize reduction factor in case of
+   * failure (the factor is reset to default if lower than 0.0001 or
+   * greater than 0.9999)
+   */
+  public void setStabilityCheck(boolean performTest,
+                                int maxIter, int maxChecks,
+                                double stabilityReduction) {
+
+    this.performTest = performTest;
+    this.maxIter     = (maxIter   <= 0) ? 2 : maxIter;
+    this.maxChecks   = (maxChecks <= 0) ? 1 : maxChecks;
+
+    if ((stabilityReduction < 0.0001) || (stabilityReduction > 0.9999)) {
+      this.stabilityReduction = 0.5;
+    } else {
+      this.stabilityReduction = stabilityReduction;
+    }
+
+  }
+
+  /** Set the step size control factors.
+
+   * <p>The new step size hNew is computed from the old one h by:
+   * <pre>
+   * hNew = h * stepControl2 / (err/stepControl1)^(1/(2k+1))
+   * </pre>
+   * where err is the scaled error and k the iteration number of the
+   * extrapolation scheme (counting from 0). The default values are
+   * 0.65 for stepControl1 and 0.94 for stepControl2.</p>
+   * <p>The step size is subject to the restriction:
+   * <pre>
+   * stepControl3^(1/(2k+1))/stepControl4 <= hNew/h <= 1/stepControl3^(1/(2k+1))
+   * </pre>
+   * The default values are 0.02 for stepControl3 and 4.0 for
+   * stepControl4.</p>
+   * @param stepControl1 first stepsize control factor (the factor is
+   * reset to default if lower than 0.0001 or greater than 0.9999)
+   * @param stepControl2 second stepsize control factor (the factor
+   * is reset to default if lower than 0.0001 or greater than 0.9999)
+   * @param stepControl3 third stepsize control factor (the factor is
+   * reset to default if lower than 0.0001 or greater than 0.9999)
+   * @param stepControl4 fourth stepsize control factor (the factor
+   * is reset to default if lower than 1.0001 or greater than 999.9)
+   */
+  public void setStepsizeControl(double stepControl1, double stepControl2,
+                                 double stepControl3, double stepControl4) {
+
+    if ((stepControl1 < 0.0001) || (stepControl1 > 0.9999)) {
+      this.stepControl1 = 0.65;
+    } else {
+      this.stepControl1 = stepControl1;
+    }
+
+    if ((stepControl2 < 0.0001) || (stepControl2 > 0.9999)) {
+      this.stepControl2 = 0.94;
+    } else {
+      this.stepControl2 = stepControl2;
+    }
+
+    if ((stepControl3 < 0.0001) || (stepControl3 > 0.9999)) {
+      this.stepControl3 = 0.02;
+    } else {
+      this.stepControl3 = stepControl3;
+    }
+
+    if ((stepControl4 < 1.0001) || (stepControl4 > 999.9)) {
+      this.stepControl4 = 4.0;
+    } else {
+      this.stepControl4 = stepControl4;
+    }
+
+  }
+
+  /** Set the order control parameters.
+   * <p>The Gragg-Bulirsch-Stoer method changes both the step size and
+   * the order during integration, in order to minimize computation
+   * cost. Each extrapolation step increases the order by 2, so the
+   * maximal order that will be used is always even, it is twice the
+   * maximal number of columns in the extrapolation table.</p>
+   * <pre>
+   * order is decreased if w(k-1) <= w(k)   * orderControl1
+   * order is increased if w(k)   <= w(k-1) * orderControl2
+   * </pre>
+   * <p>where w is the table of work per unit step for each order
+   * (number of function calls divided by the step length), and k is
+   * the current order.</p>
+   * <p>The default maximal order after construction is 18 (i.e. the
+   * maximal number of columns is 9). The default values are 0.8 for
+   * orderControl1 and 0.9 for orderControl2.</p>
+   * @param maxOrder maximal order in the extrapolation table (the
+   * maximal order is reset to default if order <= 6 or odd)
+   * @param orderControl1 first order control factor (the factor is
+   * reset to default if lower than 0.0001 or greater than 0.9999)
+   * @param orderControl2 second order control factor (the factor
+   * is reset to default if lower than 0.0001 or greater than 0.9999)
+   */
+  public void setOrderControl(int maxOrder,
+                              double orderControl1, double orderControl2) {
+
+    if ((maxOrder <= 6) || (maxOrder % 2 != 0)) {
+      this.maxOrder = 18;
+    }
+
+    if ((orderControl1 < 0.0001) || (orderControl1 > 0.9999)) {
+      this.orderControl1 = 0.8;
+    } else {
+      this.orderControl1 = orderControl1;
+    }
+
+    if ((orderControl2 < 0.0001) || (orderControl2 > 0.9999)) {
+      this.orderControl2 = 0.9;
+    } else {
+      this.orderControl2 = orderControl2;
+    }
+
+    // reinitialize the arrays
+    initializeArrays();
+
+  }
+
+  /** Set the step handler for this integrator.
+   * The handler will be called by the integrator for each accepted
+   * step.
+   * @param handler handler for the accepted steps
+   */
+  public void setStepHandler (StepHandler handler) {
+
+    super.setStepHandler(handler);
+    denseOutput = (handler.requiresDenseOutput()
+                   || (! switchesHandler.isEmpty()));
+
+    // reinitialize the arrays
+    initializeArrays();
+
+  }
+
+  /** Add a switching function to the integrator.
+   * @param function switching function
+   * @param maxCheckInterval maximal time interval between switching
+   * function checks (this interval prevents missing sign changes in
+   * case the integration steps becomes very large)
+   * @param convergence convergence threshold in the event time search
+   */
+  public void addSwitchingFunction(SwitchingFunction function,
+                                   double maxCheckInterval,
+                                   double convergence) {
+    super.addSwitchingFunction(function, maxCheckInterval, convergence);
+    denseOutput = (handler.requiresDenseOutput()
+                   || (! switchesHandler.isEmpty()));
+
+    // reinitialize the arrays
+    initializeArrays();
+
+  }
+
+  /** Initialize the integrator internal arrays. */
+  private void initializeArrays() {
+
+    int size = maxOrder / 2;
+
+    if ((sequence == null) || (sequence.length != size)) {
+      // all arrays should be reallocated with the right size
+      sequence        = new int[size];
+      costPerStep     = new int[size];
+      coeff           = new double[size][];
+      costPerTimeUnit = new double[size];
+      optimalStep     = new double[size];
+    }
+
+    if (denseOutput) {
+      // step size sequence: 2, 6, 10, 14, ...
+      for (int k = 0; k < size; ++k) {
+        sequence[k] = 4 * k + 2;
+      }
+    } else {
+      // step size sequence: 2, 4, 6, 8, ...
+      for (int k = 0; k < size; ++k) {
+        sequence[k] = 2 * (k + 1); 
+      }
+    }
+
+    // initialize the order selection cost array
+    // (number of function calls for each column of the extrapolation table)
+    costPerStep[0] = sequence[0] + 1;
+    for (int k = 1; k < size; ++k) {
+      costPerStep[k] = costPerStep[k-1] + sequence[k];
+    }
+
+    // initialize the extrapolation tables
+    for (int k = 0; k < size; ++k) {
+      coeff[k] = (k > 0) ? new double[k] : null;
+      for (int l = 0; l < k; ++l) {
+        double ratio = ((double) sequence[k]) / sequence[k-l-1];
+        coeff[k][l] = 1.0 / (ratio * ratio - 1.0);
+      }
+    }
+
+  }
+
+  /** Set the interpolation order control parameter.
+   * The interpolation order for dense output is 2k - mudif + 1. The
+   * default value for mudif is 4 and the interpolation error is used
+   * in stepsize control by default.
+
+   * @param useInterpolationError if true, interpolation error is used
+   * for stepsize control
+   * @param mudif interpolation order control parameter (the parameter
+   * is reset to default if <= 0 or >= 7)
+   */
+  public void setInterpolationControl(boolean useInterpolationError,
+                                      int mudif) {
+
+    this.useInterpolationError = useInterpolationError;
+
+    if ((mudif <= 0) || (mudif >= 7)) {
+      this.mudif = 4;
+    } else {
+      this.mudif = mudif;
+    }
+
+  }
+
+  /** Get the name of the method.
+   * @return name of the method
+   */
+  public String getName() {
+    return methodName;
+  }
+
+  /** Update scaling array.
+   * @param y1 first state vector to use for scaling
+   * @param y2 second state vector to use for scaling
+   * @param scale scaling array to update
+   */
+  private void rescale(double[] y1, double[] y2, double[] scale) {
+    if (vecAbsoluteTolerance == null) {
+      for (int i = 0; i < scale.length; ++i) {
+        double yi = Math.max(Math.abs(y1[i]), Math.abs(y2[i]));
+        scale[i] = scalAbsoluteTolerance + scalRelativeTolerance * yi;
+      }
+    } else {
+      for (int i = 0; i < scale.length; ++i) {
+        double yi = Math.max(Math.abs(y1[i]), Math.abs(y2[i]));
+        scale[i] = vecAbsoluteTolerance[i] + vecRelativeTolerance[i] * yi;
+      }
+    }
+  }
+
+  /** Perform integration over one step using substeps of a modified
+   * midpoint method.
+   * @param equations differential equations to integrate
+   * @param t0 initial time
+   * @param y0 initial value of the state vector at t0
+   * @param step global step
+   * @param k iteration number (from 0 to sequence.length - 1)
+   * @param scale scaling array
+   * @param f placeholder where to put the state vector derivatives at each substep
+   *          (element 0 already contains initial derivative)
+   * @param yMiddle placeholder where to put the state vector at the middle of the step
+   * @param yEnd placeholder where to put the state vector at the end
+   * @param yTmp placeholder for one state vector
+   * @return true if computation was done properly,
+   *         false if stability check failed before end of computation
+   * @throws DerivativeException this exception is propagated to the caller if the
+   * underlying user function triggers one
+   */
+  private boolean tryStep(FirstOrderDifferentialEquations equations,
+                          double t0, double[] y0, double step, int k,
+                          double[] scale, double[][] f,
+                          double[] yMiddle, double[] yEnd,
+                          double[] yTmp)
+    throws DerivativeException {
+
+    int    n        = sequence[k];
+    double subStep  = step / n;
+    double subStep2 = 2 * subStep;
+
+    // first substep
+    double t = t0 + subStep;
+    for (int i = 0; i < y0.length; ++i) {
+      yTmp[i] = y0[i];
+      yEnd[i] = y0[i] + subStep * f[0][i];
+    }
+    equations.computeDerivatives(t, yEnd, f[1]);
+
+    // other substeps
+    for (int j = 1; j < n; ++j) {
+
+      if (2 * j == n) {
+        // save the point at the middle of the step
+        System.arraycopy(yEnd, 0, yMiddle, 0, y0.length);
+      }
+
+      t += subStep;
+      for (int i = 0; i < y0.length; ++i) {
+        double middle = yEnd[i];
+        yEnd[i]       = yTmp[i] + subStep2 * f[j][i];
+        yTmp[i]       = middle;
+      }
+
+      equations.computeDerivatives(t, yEnd, f[j+1]);
+
+      // stability check
+      if (performTest && (j <= maxChecks) && (k < maxIter)) {
+        double initialNorm = 0.0;
+        for (int l = 0; l < y0.length; ++l) {
+          double ratio = f[0][l] / scale[l];
+          initialNorm += ratio * ratio;
+        }
+        double deltaNorm = 0.0;
+        for (int l = 0; l < y0.length; ++l) {
+          double ratio = (f[j+1][l] - f[0][l]) / scale[l];
+          deltaNorm += ratio * ratio;
+        }
+        if (deltaNorm > 4 * Math.max(1.0e-15, initialNorm)) {
+          return false;
+        }
+      }
+
+    }
+
+    // correction of the last substep (at t0 + step)
+    for (int i = 0; i < y0.length; ++i) {
+      yEnd[i] = 0.5 * (yTmp[i] + yEnd[i] + subStep * f[n][i]);
+    }
+
+    return true;
+
+  }
+
+  /** Extrapolate a vector.
+   * @param offset offset to use in the coefficients table
+   * @param k index of the last updated point
+   * @param diag working diagonal of the Aitken-Neville's
+   * triangle, without the last element
+   * @param last last element
+   */
+  private void extrapolate(int offset, int k, double[][] diag, double[] last) {
+
+    // update the diagonal
+    for (int j = 1; j < k; ++j) {
+      for (int i = 0; i < last.length; ++i) {
+        // Aitken-Neville's recursive formula
+        diag[k-j-1][i] = diag[k-j][i]
+                       + coeff[k+offset][j-1] * (diag[k-j][i] - diag[k-j-1][i]);
+      }
+    }
+
+    // update the last element
+    for (int i = 0; i < last.length; ++i) {
+      // Aitken-Neville's recursive formula
+      last[i] = diag[0][i] + coeff[k+offset][k-1] * (diag[0][i] - last[i]);
+    }
+  }
+
+  public void integrate(FirstOrderDifferentialEquations equations,
+                        double t0, double[] y0, double t, double[] y)
+  throws DerivativeException, IntegratorException {
+
+    // sanity check
+    if (equations.getDimension() != y0.length) {
+      throw new IntegratorException("dimensions mismatch: "
+                                    + "ODE problem has dimension {0}"
+                                    + ", state vector has dimension {1}",
+                                    new String[] {
+                                      Integer.toString(equations.getDimension()),
+                                      Integer.toString(y0.length)
+                                    });
+    }
+    if (Math.abs(t - t0) <= 1.0e-12 * Math.max(Math.abs(t0), Math.abs(t))) {
+      throw new IntegratorException("too small integration interval: length = {0}",
+                                    new String[] {
+                                      Double.toString(Math.abs(t - t0))
+                                    });
+    }
+
+    boolean forward = (t > t0);
+
+    // create some internal working arrays
+    double[] yDot0   = new double[y0.length];
+    double[] y1      = new double[y0.length];
+    double[] yTmp    = new double[y0.length];
+    double[] yTmpDot = new double[y0.length];
+
+    double[][] diagonal = new double[sequence.length-1][];
+    double[][] y1Diag = new double[sequence.length-1][];
+    for (int k = 0; k < sequence.length-1; ++k) {
+      diagonal[k] = new double[y0.length];
+      y1Diag[k] = new double[y0.length];
+    }
+
+    double[][][] fk  = new double[sequence.length][][];
+    for (int k = 0; k < sequence.length; ++k) {
+
+      fk[k]    = new double[sequence[k] + 1][];
+
+      // all substeps start at the same point, so share the first array
+      fk[k][0] = yDot0;
+
+      for (int l = 0; l < sequence[k]; ++l) {
+        fk[k][l+1] = new double[y0.length];
+      }
+
+    }
+
+    if (y != y0) {
+      System.arraycopy(y0, 0, y, 0, y0.length);
+    }
+
+    double[] yDot1      = null;
+    double[][] yMidDots = null;
+    if (denseOutput) {
+      yDot1    = new double[y0.length];
+      yMidDots = new double[1 + 2 * sequence.length][];
+      for (int j = 0; j < yMidDots.length; ++j) {
+        yMidDots[j] = new double[y0.length];
+      }
+    } else {
+      yMidDots    = new double[1][];
+      yMidDots[0] = new double[y0.length];
+    }
+
+    // initial scaling
+    double[] scale = new double[y0.length];
+    rescale(y, y, scale);
+
+    // initial order selection
+    double log10R = Math.log(Math.max(1.0e-10,
+                                      (vecRelativeTolerance == null)
+                                      ? scalRelativeTolerance
+                                      : vecRelativeTolerance[0]))
+                  / Math.log(10.0);
+    int targetIter = Math.max(1,
+                              Math.min(sequence.length - 2,
+                                       (int) Math.floor(0.5 - 0.6 * log10R)));
+    // set up an interpolator sharing the integrator arrays
+    AbstractStepInterpolator interpolator = null;
+    if (denseOutput || (! switchesHandler.isEmpty())) {
+      interpolator = new GraggBulirschStoerStepInterpolator(y, yDot0,
+                                                            y1, yDot1,
+                                                            yMidDots, forward);
+    } else {
+      interpolator = new DummyStepInterpolator(y, forward);
+    }
+    interpolator.storeTime(t0);
+
+    double  currentT         = t0;
+    double  hNew             = 0;
+    double  maxError         = Double.MAX_VALUE;
+    boolean previousRejected = false;
+    boolean firstTime        = true;
+    boolean newStep          = true;
+    boolean lastStep         = false;
+    boolean firstStepAlreadyComputed = false;
+    handler.reset();
+    costPerTimeUnit[0] = 0;
+    while (! lastStep) {
+
+      double h;
+      double error;
+      boolean reject = false;
+
+      if (newStep) {
+
+        interpolator.shift();
+
+        // first evaluation, at the beginning of the step
+        if (! firstStepAlreadyComputed) {
+          equations.computeDerivatives(currentT, y, yDot0);
+        }
+
+        if (firstTime) {
+
+          hNew = initializeStep(equations, forward,
+                                2 * targetIter + 1, scale,
+                                currentT, y, yDot0, yTmp, yTmpDot);
+
+          if (! forward) {
+            hNew = -hNew;
+          }
+
+        }
+
+        newStep = false;
+
+      }
+
+      h = hNew;
+
+      // step adjustment near bounds
+      if ((forward && (currentT + h > t))
+          || ((! forward) && (currentT + h < t))) {
+        h = t - currentT;
+      }
+      double nextT = currentT + h;
+      lastStep = forward ? (nextT >= t) : (nextT <= t);
+
+      // iterate over several substep sizes
+      int k = -1;
+      for (boolean loop = true; loop; ) {
+
+        ++k;
+
+        // modified midpoint integration with the current substep
+        if ( ! tryStep(equations, currentT, y, h, k, scale, fk[k],
+                       (k == 0) ? yMidDots[0] : diagonal[k-1],
+                       (k == 0) ? y1 : y1Diag[k-1],
+                       yTmp)) {
+
+          // the stability check failed, we reduce the global step
+          hNew   = Math.abs(filterStep(h * stabilityReduction, false));
+          reject = true;
+          loop   = false;
+
+        } else {
+
+          // the substep was computed successfully
+          if (k > 0) {
+
+            // extrapolate the state at the end of the step
+            // using last iteration data
+            extrapolate(0, k, y1Diag, y1);
+            rescale(y, y1, scale);
+
+            // estimate the error at the end of the step.
+            error = 0;
+            for (int j = 0; j < y0.length; ++j) {
+              double e = Math.abs(y1[j] - y1Diag[0][j]) / scale[j];
+              error += e * e;
+            }
+            error = Math.sqrt(error / y0.length);
+
+            if ((error > 1.0e15) || ((k > 1) && (error > maxError))) {
+              // error is too big, we reduce the global step
+              hNew   = Math.abs(filterStep(h * stabilityReduction, false));
+              reject = true;
+              loop   = false;
+            } else {
+
+              maxError = Math.max(4 * error, 1.0);
+
+              // compute optimal stepsize for this order
+              double exp = 1.0 / (2 * k + 1);
+              double fac = stepControl2 / Math.pow(error / stepControl1, exp);
+              double pow = Math.pow(stepControl3, exp);
+              fac = Math.max(pow / stepControl4, Math.min(1 / pow, fac));
+              optimalStep[k]     = Math.abs(filterStep(h * fac, true));
+              costPerTimeUnit[k] = costPerStep[k] / optimalStep[k];
+
+              // check convergence
+              switch (k - targetIter) {
+
+              case -1 :
+                if ((targetIter > 1) && ! previousRejected) {
+
+                  // check if we can stop iterations now
+                  if (error <= 1.0) {
+                    // convergence have been reached just before targetIter
+                    loop = false;
+                  } else {
+                    // estimate if there is a chance convergence will
+                    // be reached on next iteration, using the
+                    // asymptotic evolution of error
+                    double ratio = ((double) sequence [k] * sequence[k+1])
+                                 / (sequence[0] * sequence[0]);
+                    if (error > ratio * ratio) {
+                      // we don't expect to converge on next iteration
+                      // we reject the step immediately and reduce order
+                      reject = true;
+                      loop   = false;
+                      targetIter = k;
+                      if ((targetIter > 1)
+                          && (costPerTimeUnit[targetIter-1]
+                              < orderControl1 * costPerTimeUnit[targetIter])) {
+                        --targetIter;
+                      }
+                      hNew = optimalStep[targetIter];
+                    }
+                  }
+                }
+                break;
+
+              case 0:
+                if (error <= 1.0) {
+                  // convergence has been reached exactly at targetIter
+                  loop = false;
+                } else {
+                  // estimate if there is a chance convergence will
+                  // be reached on next iteration, using the
+                  // asymptotic evolution of error
+                  double ratio = ((double) sequence[k+1]) / sequence[0];
+                  if (error > ratio * ratio) {
+                    // we don't expect to converge on next iteration
+                    // we reject the step immediately
+                    reject = true;
+                    loop = false;
+                    if ((targetIter > 1)
+                        && (costPerTimeUnit[targetIter-1]
+                            < orderControl1 * costPerTimeUnit[targetIter])) {
+                      --targetIter;
+                    }
+                    hNew = optimalStep[targetIter];
+                  }
+                }
+                break;
+
+              case 1 :
+                if (error > 1.0) {
+                  reject = true;
+                  if ((targetIter > 1)
+                      && (costPerTimeUnit[targetIter-1]
+                          < orderControl1 * costPerTimeUnit[targetIter])) {
+                    --targetIter;
+                  }
+                  hNew = optimalStep[targetIter];
+                }
+                loop = false;
+                break;
+
+              default :
+                if ((firstTime || lastStep) && (error <= 1.0)) {
+                  loop = false;
+                }
+                break;
+
+              }
+
+            }
+          }
+        }
+      }
+
+      // dense output handling
+      double hInt = getMaxStep();
+      if (denseOutput && ! reject) {
+
+        // extrapolate state at middle point of the step
+        for (int j = 1; j <= k; ++j) {
+          extrapolate(0, j, diagonal, yMidDots[0]);
+        }
+
+        // derivative at end of step
+        equations.computeDerivatives(currentT + h, y1, yDot1);
+
+        int mu = 2 * k - mudif + 3;
+
+        for (int l = 0; l < mu; ++l) {
+
+          // derivative at middle point of the step
+          int l2 = l / 2;
+          double factor = Math.pow(0.5 * sequence[l2], l);
+          int middleIndex = fk[l2].length / 2;
+          for (int i = 0; i < y0.length; ++i) {
+            yMidDots[l+1][i] = factor * fk[l2][middleIndex + l][i];
+          }
+          for (int j = 1; j <= k - l2; ++j) {
+            factor = Math.pow(0.5 * sequence[j + l2], l);
+            middleIndex = fk[l2+j].length / 2;
+            for (int i = 0; i < y0.length; ++i) {
+              diagonal[j-1][i] = factor * fk[l2+j][middleIndex+l][i];
+            }
+            extrapolate(l2, j, diagonal, yMidDots[l+1]);
+          }
+          for (int i = 0; i < y0.length; ++i) {
+            yMidDots[l+1][i] *= h;
+          }
+
+          // compute centered differences to evaluate next derivatives
+          for (int j = (l + 1) / 2; j <= k; ++j) {
+            for (int m = fk[j].length - 1; m >= 2 * (l + 1); --m) {
+              for (int i = 0; i < y0.length; ++i) {
+                fk[j][m][i] -= fk[j][m-2][i];
+              }
+            }
+          }
+
+        }
+
+        if (mu >= 0) {
+
+          // estimate the dense output coefficients
+          GraggBulirschStoerStepInterpolator gbsInterpolator
+            = (GraggBulirschStoerStepInterpolator) interpolator;
+          gbsInterpolator.computeCoefficients(mu, h);
+
+          if (useInterpolationError) {
+            // use the interpolation error to limit stepsize
+            double interpError = gbsInterpolator.estimateError(scale);
+            hInt = Math.abs(h / Math.max(Math.pow(interpError, 1.0 / (mu+4)),
+                                         0.01));
+            if (interpError > 10.0) {
+              hNew = hInt;
+              reject = true;
+            }
+          }
+
+          // Switching functions handling
+          if (!reject) {
+            interpolator.storeTime(currentT + h);
+            if (switchesHandler.evaluateStep(interpolator)) {
+              reject = true;
+              hNew = Math.abs(switchesHandler.getEventTime() - currentT);
+            }
+          }
+
+        }
+
+        if (!reject) {
+          // we will reuse the slope for the beginning of next step
+          firstStepAlreadyComputed = true;
+          System.arraycopy(yDot1, 0, yDot0, 0, y0.length);
+        }
+
+      }
+
+      if (! reject) {
+
+        // store end of step state
+        currentT += h;
+        System.arraycopy(y1, 0, y, 0, y0.length);
+
+        switchesHandler.stepAccepted(currentT, y);
+        if (switchesHandler.stop()) {
+          lastStep = true;
+        }
+
+        // provide the step data to the step handler
+        interpolator.storeTime(currentT);
+        handler.handleStep(interpolator, lastStep);
+
+        if (switchesHandler.reset(currentT, y) && ! lastStep) {
+          // some switching function has triggered changes that
+          // invalidate the derivatives, we need to recompute them
+          firstStepAlreadyComputed = false;
+        }
+
+        int optimalIter;
+        if (k == 1) {
+          optimalIter = 2;
+          if (previousRejected) {
+            optimalIter = 1;
+          }
+        } else if (k <= targetIter) {
+          optimalIter = k;
+          if (costPerTimeUnit[k-1] < orderControl1 * costPerTimeUnit[k]) {
+            optimalIter = k-1;
+          } else if (costPerTimeUnit[k] < orderControl2 * costPerTimeUnit[k-1]) {
+            optimalIter = Math.min(k+1, sequence.length - 2);
+          }
+        } else {
+          optimalIter = k - 1;
+          if ((k > 2)
+              && (costPerTimeUnit[k-2] < orderControl1 * costPerTimeUnit[k-1])) {
+            optimalIter = k - 2;
+          }
+          if (costPerTimeUnit[k] < orderControl2 * costPerTimeUnit[optimalIter]) {
+            optimalIter = Math.min(k, sequence.length - 2);
+          }
+        }
+
+        if (previousRejected) {
+          // after a rejected step neither order nor stepsize
+          // should increase
+          targetIter = Math.min(optimalIter, k);
+          hNew = Math.min(Math.abs(h), optimalStep[targetIter]);
+        } else {
+          // stepsize control
+          if (optimalIter <= k) {
+            hNew = optimalStep[optimalIter];
+          } else {
+            if ((k < targetIter)
+                && (costPerTimeUnit[k] < orderControl2 * costPerTimeUnit[k-1])) {
+              hNew = filterStep(optimalStep[k]
+                                * costPerStep[optimalIter+1] / costPerStep[k],
+                                false);
+            } else {
+              hNew = filterStep(optimalStep[k]
+                                * costPerStep[optimalIter] / costPerStep[k],
+                                false);
+            }
+          }
+
+          targetIter = optimalIter;
+
+        }
+
+        newStep = true;
+
+      }
+
+      hNew = Math.min(hNew, hInt);
+      if (! forward) {
+        hNew = -hNew;
+      }
+
+      firstTime = false;
+
+      if (reject) {
+        lastStep = false;
+        previousRejected = true;
+      } else {
+        previousRejected = false;
+      }
+
+    }
+
+  }
+
+  /** maximal order. */
+  private int maxOrder;
+
+  /** step size sequence. */
+  private int[] sequence;
+
+  /** overall cost of applying step reduction up to iteration k+1,
+   *  in number of calls.
+   */
+  private int[] costPerStep;
+
+  /** cost per unit step. */
+  private double[] costPerTimeUnit;
+
+  /** optimal steps for each order. */
+  private double[] optimalStep;
+
+  /** extrapolation coefficients. */
+  private double[][] coeff;
+
+  /** stability check enabling parameter. */
+  private boolean performTest;
+
+  /** maximal number of checks for each iteration. */
+  private int maxChecks;
+
+  /** maximal number of iterations for which checks are performed. */
+  private int maxIter;
+
+  /** stepsize reduction factor in case of stability check failure. */
+  private double stabilityReduction;
+
+  /** first stepsize control factor. */
+  private double stepControl1;
+
+  /** second stepsize control factor. */
+  private double stepControl2;
+
+  /** third stepsize control factor. */
+  private double stepControl3;
+
+  /** fourth stepsize control factor. */
+  private double stepControl4;
+
+  /** first order control factor. */
+  private double orderControl1;
+
+  /** second order control factor. */
+  private double orderControl2;
+
+  /** dense outpute required. */
+  private boolean denseOutput;
+
+  /** use interpolation error in stepsize control. */
+  private boolean useInterpolationError;
+
+  /** interpolation order control parameter. */
+  private int mudif;
+
+}

Propchange: jakarta/commons/proper/math/trunk/src/java/org/apache/commons/math/ode/GraggBulirschStoerIntegrator.java
------------------------------------------------------------------------------
    svn:eol-style = native

Added: jakarta/commons/proper/math/trunk/src/java/org/apache/commons/math/ode/GraggBulirschStoerStepInterpolator.java
URL: http://svn.apache.org/viewvc/jakarta/commons/proper/math/trunk/src/java/org/apache/commons/math/ode/GraggBulirschStoerStepInterpolator.java?view=auto&rev=512066
==============================================================================
--- jakarta/commons/proper/math/trunk/src/java/org/apache/commons/math/ode/GraggBulirschStoerStepInterpolator.java (added)
+++ jakarta/commons/proper/math/trunk/src/java/org/apache/commons/math/ode/GraggBulirschStoerStepInterpolator.java Mon Feb 26 15:12:40 2007
@@ -0,0 +1,402 @@
+/*
+ * 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.math.ode;
+
+import java.io.ObjectInput;
+import java.io.ObjectOutput;
+import java.io.IOException;
+
+/**
+ * This class implements an interpolator for the Gragg-Bulirsch-Stoer
+ * integrator.
+
+ * <p>This interpolator compute dense output inside the last step
+ * produced by a Gragg-Bulirsch-Stoer integrator.</p>
+
+ * <p>
+ * This implementation is basically a reimplementation in Java of the
+ * <a
+ * href="http://www.unige.ch/math/folks/hairer/prog/nonstiff/odex.f">odex</a>
+ * fortran code by E. Hairer and G. Wanner. The redistribution policy
+ * for this code is available <a
+ * href="http://www.unige.ch/~hairer/prog/licence.txt">here</a>, for
+ * convenience, it is reproduced below.</p>
+ * </p>
+
+ * <table border="0" width="80%" cellpadding="10" align="center" bgcolor="#E0E0E0">
+ * <tr><td>Copyright (c) 2004, Ernst Hairer</td></tr>
+
+ * <tr><td>Redistribution and use in source and binary forms, with or
+ * without modification, are permitted provided that the following
+ * conditions are met:
+ * <ul>
+ *  <li>Redistributions of source code must retain the above copyright
+ *      notice, this list of conditions and the following disclaimer.</li>
+ *  <li>Redistributions in binary form must reproduce the above copyright
+ *      notice, this list of conditions and the following disclaimer in the
+ *      documentation and/or other materials provided with the distribution.</li>
+ * </ul></td></tr>
+
+ * <tr><td><strong>THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND
+ * CONTRIBUTORS "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING,
+ * BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS
+ * FOR A  PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE REGENTS OR
+ * CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
+ * EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
+ * PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
+ * PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
+ * LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
+ * NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
+ * SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.</strong></td></tr>
+ * </table>
+
+ * @see GraggBulirschStoerIntegrator
+
+ * @version $Id: GraggBulirschStoerStepInterpolator.java 1705 2006-09-17 19:57:39Z luc $
+ * @author E. Hairer and G. Wanner (fortran version)
+ * @author L. Maisonobe (Java port)
+
+ */
+
+class GraggBulirschStoerStepInterpolator
+  extends AbstractStepInterpolator {
+
+  /** Slope at the beginning of the step. */
+  private double[] y0Dot;
+
+  /** State at the end of the step. */
+  private double[] y1;
+
+  /** Slope at the end of the step. */
+  private double[] y1Dot;
+
+  /** Derivatives at the middle of the step.
+   * element 0 is state at midpoint, element 1 is first derivative ...
+   */
+  private double[][] yMidDots;
+
+  /** Interpolation polynoms. */
+  private double[][] polynoms;
+
+  /** Error coefficients for the interpolation. */
+  private double[] errfac;
+
+  /** Degree of the interpolation polynoms. */
+  private int currentDegree;
+
+  /** Reallocate the internal tables.
+   * Reallocate the internal tables in order to be able to handle
+   * interpolation polynoms up to the given degree
+   * @param maxDegree maximal degree to handle
+   */
+  private void resetTables(int maxDegree) {
+
+    if (maxDegree < 0) {
+      polynoms      = null;
+      errfac        = null;
+      currentDegree = -1;
+    } else {
+
+      double[][] newPols = new double[maxDegree + 1][];
+      if (polynoms != null) {
+        System.arraycopy(polynoms, 0, newPols, 0, polynoms.length);
+        for (int i = polynoms.length; i < newPols.length; ++i) {
+          newPols[i] = new double[currentState.length];
+        }
+      } else {
+        for (int i = 0; i < newPols.length; ++i) {
+          newPols[i] = new double[currentState.length];
+        }
+      }
+      polynoms = newPols;
+
+      // initialize the error factors array for interpolation
+      if (maxDegree <= 4) {
+        errfac = null;
+      } else {
+        errfac = new double[maxDegree - 4];
+        for (int i = 0; i < errfac.length; ++i) {
+          int ip5 = i + 5;
+          errfac[i] = 1.0 / (ip5 * ip5);
+          double e = 0.5 * Math.sqrt (((double) (i + 1)) / ip5);
+          for (int j = 0; j <= i; ++j) {
+            errfac[i] *= e / (j + 1);
+          }
+        }
+      }
+
+      currentDegree = 0;
+
+    }
+
+  }
+
+  /** Simple constructor.
+   * This constructor should not be used directly, it is only intended
+   * for the serialization process.
+   */
+  public GraggBulirschStoerStepInterpolator() {
+    y0Dot    = null;
+    y1       = null;
+    y1Dot    = null;
+    yMidDots = null;
+    resetTables(-1);
+  }
+
+  /** Simple constructor.
+   * @param y reference to the integrator array holding the current state
+   * @param y0Dot reference to the integrator array holding the slope
+   * at the beginning of the step
+   * @param y1 reference to the integrator array holding the state at
+   * the end of the step
+   * @param y1Dot reference to the integrator array holding the slope
+   * at theend of the step
+   * @param yMidDots reference to the integrator array holding the
+   * derivatives at the middle point of the step
+   * @param forward integration direction indicator
+   */
+  public GraggBulirschStoerStepInterpolator(double[] y, double[] y0Dot,
+                                            double[] y1, double[] y1Dot,
+                                            double[][] yMidDots,
+                                            boolean forward) {
+
+    super(y, forward);
+    this.y0Dot    = y0Dot;
+    this.y1       = y1;
+    this.y1Dot    = y1Dot;
+    this.yMidDots = yMidDots;
+
+    resetTables(yMidDots.length + 4);
+
+  }
+
+  /** Copy constructor.
+   * @param interpolator interpolator to copy from. The copy is a deep
+   * copy: its arrays are separated from the original arrays of the
+   * instance
+   */
+  public GraggBulirschStoerStepInterpolator
+    (GraggBulirschStoerStepInterpolator interpolator) {
+
+    super(interpolator);
+
+    int dimension = currentState.length;
+
+    // the interpolator has been finalized,
+    // the following arrays are not needed anymore
+    y0Dot    = null;
+    y1       = null;
+    y1Dot    = null;
+    yMidDots = null;
+
+    // copy the interpolation polynoms (up to the current degree only)
+    if (interpolator.polynoms == null) {
+      polynoms = null;
+      currentDegree = -1;
+    } else {
+      resetTables(interpolator.currentDegree);
+      for (int i = 0; i < polynoms.length; ++i) {
+        polynoms[i] = new double[dimension];
+        System.arraycopy(interpolator.polynoms[i], 0,
+                         polynoms[i], 0, dimension);
+      }
+      currentDegree = interpolator.currentDegree;
+    }
+
+  }
+
+  /**
+   * Clone the instance.
+   * the copy is a deep copy: its arrays are separated from the
+   * original arrays of the instance
+   * @return a copy of the instance
+   */
+  public Object clone() {
+    return new GraggBulirschStoerStepInterpolator(this);
+  }
+
+  /** Compute the interpolation coefficients for dense output.
+   * @param mu degree of the interpolation polynom
+   * @param h current step
+   */
+  public void computeCoefficients(int mu, double h) {
+
+    if ((polynoms == null) || (polynoms.length <= (mu + 4))) {
+      resetTables(mu + 4);
+    }
+
+    currentDegree = mu + 4;
+
+    for (int i = 0; i < currentState.length; ++i) {
+
+      double yp0   = h * y0Dot[i];
+      double yp1   = h * y1Dot[i];
+      double ydiff = y1[i] - currentState[i];
+      double aspl  = ydiff - yp1;
+      double bspl  = yp0 - ydiff;
+
+      polynoms[0][i] = currentState[i];
+      polynoms[1][i] = ydiff;
+      polynoms[2][i] = aspl;
+      polynoms[3][i] = bspl;
+
+      if (mu < 0) {
+        return;
+      }
+
+      // compute the remaining coefficients
+      double ph0 = 0.5 * (currentState[i] + y1[i]) + 0.125 * (aspl + bspl);
+      polynoms[4][i] = 16 * (yMidDots[0][i] - ph0);
+
+      if (mu > 0) {
+        double ph1 = ydiff + 0.25 * (aspl - bspl);
+        polynoms[5][i] = 16 * (yMidDots[1][i] - ph1);
+
+        if (mu > 1) {
+          double ph2 = yp1 - yp0;
+          polynoms[6][i] = 16 * (yMidDots[2][i] - ph2 + polynoms[4][i]);
+
+          if (mu > 2) {
+            double ph3 = 6 * (bspl - aspl);
+            polynoms[7][i] = 16 * (yMidDots[3][i] - ph3 + 3 * polynoms[5][i]);
+
+            for (int j = 4; j <= mu; ++j) {
+              double fac1 = 0.5 * j * (j - 1);
+              double fac2 = 2 * fac1 * (j - 2) * (j - 3);
+              polynoms[j+4][i] = 16 * (yMidDots[j][i]
+                                       + fac1 * polynoms[j+2][i]
+                                       - fac2 * polynoms[j][i]);
+            }
+
+          }
+        }
+      }
+    }
+
+  }
+
+  /** Estimate interpolation error.
+   * @param scale scaling array
+   * @return estimate of the interpolation error
+   */
+  public double estimateError(double[] scale) {
+    double error = 0;
+    if (currentDegree >= 5) {
+      for (int i = 0; i < currentState.length; ++i) {
+        double e = polynoms[currentDegree][i] / scale[i];
+        error += e * e;
+      }
+      error = Math.sqrt(error / currentState.length) * errfac[currentDegree-5];
+    }
+    return error;
+  }
+
+  /** Compute the state at the interpolated time.
+   * This is the main processing method that should be implemented by
+   * the derived classes to perform the interpolation.
+   * @param theta normalized interpolation abscissa within the step
+   * (theta is zero at the previous time step and one at the current time step)
+   * @param oneMinusThetaH time gap between the interpolated time and
+   * the current time
+   * @throws DerivativeException this exception is propagated to the caller if the
+   * underlying user function triggers one
+   */
+  protected void computeInterpolatedState(double theta,
+                                          double oneMinusThetaH)
+    throws DerivativeException {
+
+    int dimension = currentState.length;
+
+    double oneMinusTheta = 1.0 - theta;
+    double theta05       = theta - 0.5;
+    double t4            = theta * oneMinusTheta;
+    t4 = t4 * t4;
+
+    for (int i = 0; i < dimension; ++i) {
+      interpolatedState[i] = polynoms[0][i]
+        + theta * (polynoms[1][i]
+                   + oneMinusTheta * (polynoms[2][i] * theta
+                                      + polynoms[3][i] * oneMinusTheta));
+
+      if (currentDegree > 3) {
+        double c = polynoms[currentDegree][i];
+        for (int j = currentDegree - 1; j > 3; --j) {
+          c = polynoms[j][i] + c * theta05 / (j - 3);
+        }
+        interpolatedState[i] += t4 * c;
+      }
+    }
+
+  }
+    
+  /** Save the state of the instance.
+   * @param out stream where to save the state
+   * @exception IOException in case of write error
+   */
+  public void writeExternal(ObjectOutput out)
+    throws IOException {
+
+    int dimension = currentState.length;
+
+    // save the state of the base class
+    writeBaseExternal(out);
+
+    // save the local attributes (but not the temporary vectors)
+    out.writeInt(currentDegree);
+    for (int k = 0; k <= currentDegree; ++k) {
+      for (int l = 0; l < dimension; ++l) {
+        out.writeDouble(polynoms[k][l]);
+      }
+    }
+
+  }
+
+  /** Read the state of the instance.
+   * @param in stream where to read the state from
+   * @exception IOException in case of read error
+   */
+  public void readExternal(ObjectInput in)
+    throws IOException {
+
+    // read the base class 
+    double t = readBaseExternal(in);
+    int dimension = currentState.length;
+
+    // read the local attributes
+    int degree = in.readInt();
+    resetTables(degree);
+    currentDegree = degree;
+
+    for (int k = 0; k <= currentDegree; ++k) {
+      for (int l = 0; l < dimension; ++l) {
+        polynoms[k][l] = in.readDouble();
+      }
+    }
+
+    try {
+      // we can now set the interpolated time and state
+      setInterpolatedTime(t);
+    } catch (DerivativeException e) {
+      throw new IOException(e.getMessage());
+    }
+
+  }
+
+  private static final long serialVersionUID = 7320613236731409847L;
+
+}

Propchange: jakarta/commons/proper/math/trunk/src/java/org/apache/commons/math/ode/GraggBulirschStoerStepInterpolator.java
------------------------------------------------------------------------------
    svn:eol-style = native

Added: jakarta/commons/proper/math/trunk/src/java/org/apache/commons/math/ode/HighamHall54Integrator.java
URL: http://svn.apache.org/viewvc/jakarta/commons/proper/math/trunk/src/java/org/apache/commons/math/ode/HighamHall54Integrator.java?view=auto&rev=512066
==============================================================================
--- jakarta/commons/proper/math/trunk/src/java/org/apache/commons/math/ode/HighamHall54Integrator.java (added)
+++ jakarta/commons/proper/math/trunk/src/java/org/apache/commons/math/ode/HighamHall54Integrator.java Mon Feb 26 15:12:40 2007
@@ -0,0 +1,139 @@
+// 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.math.ode;
+
+/**
+ * This class implements the 5(4) Higham and Hall integrator for
+ * Ordinary Differential Equations.
+
+ * <p>This integrator is an embedded Runge-Kutta-Fehlberg integrator
+ * of order 5(4) used in local extrapolation mode (i.e. the solution
+ * is computed using the high order formula) with stepsize control
+ * (and automatic step initialization) and continuous output. This
+ * method uses 7 functions evaluations per step.</p>
+
+ * @version $Id: HighamHall54Integrator.java 1705 2006-09-17 19:57:39Z luc $
+ * @author L. Maisonobe
+
+ */
+
+public class HighamHall54Integrator
+  extends RungeKuttaFehlbergIntegrator {
+
+  private static final String methodName = "Higham-Hall 5(4)";
+
+  private static final double[] c = {
+    2.0/9.0, 1.0/3.0, 1.0/2.0, 3.0/5.0, 1.0, 1.0
+  };
+
+  private static final double[][] a = {
+    {2.0/9.0},
+    {1.0/12.0, 1.0/4.0},
+    {1.0/8.0, 0.0, 3.0/8.0},
+    {91.0/500.0, -27.0/100.0, 78.0/125.0, 8.0/125.0},
+    {-11.0/20.0, 27.0/20.0, 12.0/5.0, -36.0/5.0, 5.0},
+    {1.0/12.0, 0.0, 27.0/32.0, -4.0/3.0, 125.0/96.0, 5.0/48.0}
+  };
+
+  private static final double[] b = {
+    1.0/12.0, 0.0, 27.0/32.0, -4.0/3.0, 125.0/96.0, 5.0/48.0, 0.0
+  };
+
+  private static final double[] e = {
+    -1.0/20.0, 0.0, 81.0/160.0, -6.0/5.0, 25.0/32.0, 1.0/16.0, -1.0/10.0
+  };
+
+  /** Simple constructor.
+   * Build a fifth order Higham and Hall integrator with the given step bounds
+   * @param minStep minimal step (must be positive even for backward
+   * integration), the last step can be smaller than this
+   * @param maxStep maximal step (must be positive even for backward
+   * integration)
+   * @param scalAbsoluteTolerance allowed absolute error
+   * @param scalRelativeTolerance allowed relative error
+   */
+  public HighamHall54Integrator(double minStep, double maxStep,
+                                double scalAbsoluteTolerance,
+                                double scalRelativeTolerance) {
+    super(false, c, a, b, new HighamHall54StepInterpolator(),
+          minStep, maxStep, scalAbsoluteTolerance, scalRelativeTolerance);
+  }
+
+  /** Simple constructor.
+   * Build a fifth order Higham and Hall integrator with the given step bounds
+   * @param minStep minimal step (must be positive even for backward
+   * integration), the last step can be smaller than this
+   * @param maxStep maximal step (must be positive even for backward
+   * integration)
+   * @param vecAbsoluteTolerance allowed absolute error
+   * @param vecRelativeTolerance allowed relative error
+   */
+  public HighamHall54Integrator(double minStep, double maxStep,
+                                double[] vecAbsoluteTolerance,
+                                double[] vecRelativeTolerance) {
+    super(false, c, a, b, new HighamHall54StepInterpolator(),
+          minStep, maxStep, vecAbsoluteTolerance, vecRelativeTolerance);
+  }
+
+  /** Get the name of the method.
+   * @return name of the method
+   */
+  public String getName() {
+    return methodName;
+  }
+
+  /** Get the order of the method.
+   * @return order of the method
+   */
+  public int getOrder() {
+    return 5;
+  }
+
+  /** Compute the error ratio.
+   * @param yDotK derivatives computed during the first stages
+   * @param y0 estimate of the step at the start of the step
+   * @param y1 estimate of the step at the end of the step
+   * @param h  current step
+   * @return error ratio, greater than 1 if step should be rejected
+   */
+  protected double estimateError(double[][] yDotK,
+                                 double[] y0, double[] y1,
+                                 double h) {
+
+    double error = 0;
+
+    for (int j = 0; j < y0.length; ++j) {
+      double errSum = e[0] * yDotK[0][j];
+      for (int l = 1; l < e.length; ++l) {
+        errSum += e[l] * yDotK[l][j];
+      }
+
+      double yScale = Math.max(Math.abs(y0[j]), Math.abs(y1[j]));
+      double tol = (vecAbsoluteTolerance == null)
+        ? (scalAbsoluteTolerance + scalRelativeTolerance * yScale)
+        : (vecAbsoluteTolerance[j] + vecRelativeTolerance[j] * yScale);
+      double ratio  = h * errSum / tol;
+      error += ratio * ratio;
+
+    }
+
+    return Math.sqrt(error / y0.length);
+
+  }
+
+}

Propchange: jakarta/commons/proper/math/trunk/src/java/org/apache/commons/math/ode/HighamHall54Integrator.java
------------------------------------------------------------------------------
    svn:eol-style = native

Added: jakarta/commons/proper/math/trunk/src/java/org/apache/commons/math/ode/HighamHall54StepInterpolator.java
URL: http://svn.apache.org/viewvc/jakarta/commons/proper/math/trunk/src/java/org/apache/commons/math/ode/HighamHall54StepInterpolator.java?view=auto&rev=512066
==============================================================================
--- jakarta/commons/proper/math/trunk/src/java/org/apache/commons/math/ode/HighamHall54StepInterpolator.java (added)
+++ jakarta/commons/proper/math/trunk/src/java/org/apache/commons/math/ode/HighamHall54StepInterpolator.java Mon Feb 26 15:12:40 2007
@@ -0,0 +1,96 @@
+// 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.math.ode;
+
+/**
+ * This class represents an interpolator over the last step during an
+ * ODE integration for the 5(4) Higham and Hall integrator.
+ *
+ * @see HighamHall54Integrator
+ *
+ * @version $Id: HighamHall54StepInterpolator.java 1705 2006-09-17 19:57:39Z luc $
+ * @author L. Maisonobe
+ *
+ */
+
+class HighamHall54StepInterpolator
+  extends RungeKuttaStepInterpolator {
+
+  /** Simple constructor.
+   * This constructor builds an instance that is not usable yet, the
+   * {@link AbstractStepInterpolator#reinitialize} method should be called
+   * before using the instance in order to initialize the internal arrays. This
+   * constructor is used only in order to delay the initialization in
+   * some cases. The {@link RungeKuttaFehlbergIntegrator} uses the
+   * prototyping design pattern to create the step interpolators by
+   * cloning an uninitialized model and latter initializing the copy.
+   */
+  public HighamHall54StepInterpolator() {
+    super();
+  }
+
+  /** Copy constructor.
+   * @param interpolator interpolator to copy from. The copy is a deep
+   * copy: its arrays are separated from the original arrays of the
+   * instance
+   */
+  public HighamHall54StepInterpolator(HighamHall54StepInterpolator interpolator) {
+    super(interpolator);
+  }
+
+  /**
+   * Clone the instance.
+   * the copy is a deep copy: its arrays are separated from the
+   * original arrays of the instance
+   * @return a copy of the instance
+   */
+  public Object clone() {
+    return new HighamHall54StepInterpolator(this);
+  }
+
+  /** Compute the state at the interpolated time.
+   * @param theta normalized interpolation abscissa within the step
+   * (theta is zero at the previous time step and one at the current time step)
+   * @param oneMinusThetaH time gap between the interpolated time and
+   * the current time
+   * @throws DerivativeException this exception is propagated to the caller if the
+   * underlying user function triggers one
+   */
+  protected void computeInterpolatedState(double theta,
+                                          double oneMinusThetaH)
+    throws DerivativeException {
+
+    double theta2 = theta * theta;
+
+    double b0 = h * (-1.0/12.0 + theta * (1.0 + theta * (-15.0/4.0 + theta * (16.0/3.0 + theta * -5.0/2.0))));
+    double b2 = h * (-27.0/32.0 + theta2 * (459.0/32.0 + theta * (-243.0/8.0 + theta * 135.0/8.0)));
+    double b3 = h * (4.0/3.0 + theta2 * (-22.0 + theta * (152.0/3.0  + theta * -30.0)));
+    double b4 = h * (-125.0/96.0 + theta2 * (375.0/32.0 + theta * (-625.0/24.0 + theta * 125.0/8.0)));
+    double b5 = h * (-5.0/48.0 + theta2 * (-5.0/16.0 + theta * 5.0/12.0));
+
+    for (int i = 0; i < interpolatedState.length; ++i) {
+      interpolatedState[i] = currentState[i]
+                           + b0 * yDotK[0][i] + b2 * yDotK[2][i] + b3 * yDotK[3][i]
+                           + b4 * yDotK[4][i] + b5 * yDotK[5][i];
+    }
+
+  }
+
+  private static final long serialVersionUID = -3583240427587318654L;
+
+}

Propchange: jakarta/commons/proper/math/trunk/src/java/org/apache/commons/math/ode/HighamHall54StepInterpolator.java
------------------------------------------------------------------------------
    svn:eol-style = native

Added: jakarta/commons/proper/math/trunk/src/java/org/apache/commons/math/ode/IntegratorException.java
URL: http://svn.apache.org/viewvc/jakarta/commons/proper/math/trunk/src/java/org/apache/commons/math/ode/IntegratorException.java?view=auto&rev=512066
==============================================================================
--- jakarta/commons/proper/math/trunk/src/java/org/apache/commons/math/ode/IntegratorException.java (added)
+++ jakarta/commons/proper/math/trunk/src/java/org/apache/commons/math/ode/IntegratorException.java Mon Feb 26 15:12:40 2007
@@ -0,0 +1,42 @@
+// 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.math.ode;
+
+import org.apache.commons.math.MathException;
+
+/**
+ * This exception is made available to users to report
+ * the error conditions that are triggered during integration
+ * @author Luc Maisonobe
+ * @version $Id: IntegratorException.java 1705 2006-09-17 19:57:39Z luc $
+ */
+public class IntegratorException
+  extends MathException {
+    
+  /** Simple constructor.
+   * Build an exception by translating and formating a message
+   * @param specifier format specifier (to be translated)
+   * @param parts to insert in the format (no translation)
+   */
+  public IntegratorException(String specifier, String[] parts) {
+    super(specifier, parts);
+  }
+
+  private static final long serialVersionUID = -1390328069787882608L;
+
+}

Propchange: jakarta/commons/proper/math/trunk/src/java/org/apache/commons/math/ode/IntegratorException.java
------------------------------------------------------------------------------
    svn:eol-style = native

Added: jakarta/commons/proper/math/trunk/src/java/org/apache/commons/math/ode/MidpointIntegrator.java
URL: http://svn.apache.org/viewvc/jakarta/commons/proper/math/trunk/src/java/org/apache/commons/math/ode/MidpointIntegrator.java?view=auto&rev=512066
==============================================================================
--- jakarta/commons/proper/math/trunk/src/java/org/apache/commons/math/ode/MidpointIntegrator.java (added)
+++ jakarta/commons/proper/math/trunk/src/java/org/apache/commons/math/ode/MidpointIntegrator.java Mon Feb 26 15:12:40 2007
@@ -0,0 +1,75 @@
+// 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.math.ode;
+
+/**
+ * This class implements a second order Runge-Kutta integrator for
+ * Ordinary Differential Equations.
+
+ * <p>This method is an explicit Runge-Kutta method, its Butcher-array
+ * is the following one :
+ * <pre>
+ *    0  |  0    0
+ *   1/2 | 1/2   0
+ *       |----------
+ *       |  0    1
+ * </pre>
+ * </p>
+
+ * @see EulerIntegrator
+ * @see ClassicalRungeKuttaIntegrator
+ * @see GillIntegrator
+
+ * @version $Id: MidpointIntegrator.java 1705 2006-09-17 19:57:39Z luc $
+ * @author L. Maisonobe
+
+ */
+
+public class MidpointIntegrator
+  extends RungeKuttaIntegrator {
+
+  private static final String methodName = "midpoint";
+
+  private static final double[] c = {
+    1.0 / 2.0
+  };
+
+  private static final double[][] a = {
+    { 1.0 / 2.0 }
+  };
+
+  private static final double[] b = {
+    0.0, 1.0
+  };
+
+  /** Simple constructor.
+   * Build a midpoint integrator with the given step.
+   * @param step integration step
+   */
+  public MidpointIntegrator(double step) {
+    super(false, c, a, b, new MidpointStepInterpolator(), step);
+  }
+
+  /** Get the name of the method.
+   * @return name of the method
+   */
+  public String getName() {
+    return methodName;
+  }
+
+}

Propchange: jakarta/commons/proper/math/trunk/src/java/org/apache/commons/math/ode/MidpointIntegrator.java
------------------------------------------------------------------------------
    svn:eol-style = native

Added: jakarta/commons/proper/math/trunk/src/java/org/apache/commons/math/ode/MidpointStepInterpolator.java
URL: http://svn.apache.org/viewvc/jakarta/commons/proper/math/trunk/src/java/org/apache/commons/math/ode/MidpointStepInterpolator.java?view=auto&rev=512066
==============================================================================
--- jakarta/commons/proper/math/trunk/src/java/org/apache/commons/math/ode/MidpointStepInterpolator.java (added)
+++ jakarta/commons/proper/math/trunk/src/java/org/apache/commons/math/ode/MidpointStepInterpolator.java Mon Feb 26 15:12:40 2007
@@ -0,0 +1,103 @@
+// 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.math.ode;
+
+/**
+ * This class implements a step interpolator for second order
+ * Runge-Kutta integrator.
+
+ * <p>This interpolator allow to compute dense output inside the last
+ * step computed. The interpolation equation is consistent with the
+ * integration scheme :
+
+ * <pre>
+ *   y(t_n + theta h) = y (t_n + h) + (1-theta) h [theta y'_1 - (1+theta) y'_2]
+ * </pre>
+
+ * where theta belongs to [0 ; 1] and where y'_1 and y'_2 are the two
+ * evaluations of the derivatives already computed during the
+ * step.</p>
+
+ * @see MidpointIntegrator
+
+ * @version $Id: MidpointStepInterpolator.java 1705 2006-09-17 19:57:39Z luc $
+ * @author L. Maisonobe
+
+ */
+
+class MidpointStepInterpolator
+  extends RungeKuttaStepInterpolator {
+    
+  /** Simple constructor.
+   * This constructor builds an instance that is not usable yet, the
+   * {@link AbstractStepInterpolator#reinitialize} method should be called
+   * before using the instance in order to initialize the internal arrays. This
+   * constructor is used only in order to delay the initialization in
+   * some cases. The {@link RungeKuttaIntegrator} class uses the
+   * prototyping design pattern to create the step interpolators by
+   * cloning an uninitialized model and latter initializing the copy.
+   */
+  public MidpointStepInterpolator() {
+  }
+
+  /** Copy constructor.
+   * @param interpolator interpolator to copy from. The copy is a deep
+   * copy: its arrays are separated from the original arrays of the
+   * instance
+   */
+  public MidpointStepInterpolator(MidpointStepInterpolator interpolator) {
+    super(interpolator);
+  }
+
+  /**
+   * Clone the instance.
+   * the copy is a deep copy: its arrays are separated from the
+   * original arrays of the instance
+   * @return a copy of the instance
+   */
+  public Object clone() {
+    return new MidpointStepInterpolator(this);
+  }
+
+  /** Compute the state at the interpolated time.
+   * This is the main processing method that should be implemented by
+   * the derived classes to perform the interpolation.
+   * @param theta normalized interpolation abscissa within the step
+   * (theta is zero at the previous time step and one at the current time step)
+   * @param oneMinusThetaH time gap between the interpolated time and
+   * the current time
+   * @throws DerivativeException this exception is propagated to the caller if the
+   * underlying user function triggers one
+   */
+  protected void computeInterpolatedState(double theta,
+                                          double oneMinusThetaH)
+    throws DerivativeException {
+
+    double coeff1 = oneMinusThetaH * theta;
+    double coeff2 = oneMinusThetaH * (1.0 + theta);
+
+    for (int i = 0; i < interpolatedState.length; ++i) {
+      interpolatedState[i] = currentState[i]
+                           + coeff1 * yDotK[0][i] - coeff2 * yDotK[1][i];
+    }
+
+  }
+
+  private static final long serialVersionUID = -865524111506042509L;
+
+}

Propchange: jakarta/commons/proper/math/trunk/src/java/org/apache/commons/math/ode/MidpointStepInterpolator.java
------------------------------------------------------------------------------
    svn:eol-style = native

Added: jakarta/commons/proper/math/trunk/src/java/org/apache/commons/math/ode/RungeKuttaFehlbergIntegrator.java
URL: http://svn.apache.org/viewvc/jakarta/commons/proper/math/trunk/src/java/org/apache/commons/math/ode/RungeKuttaFehlbergIntegrator.java?view=auto&rev=512066
==============================================================================
--- jakarta/commons/proper/math/trunk/src/java/org/apache/commons/math/ode/RungeKuttaFehlbergIntegrator.java (added)
+++ jakarta/commons/proper/math/trunk/src/java/org/apache/commons/math/ode/RungeKuttaFehlbergIntegrator.java Mon Feb 26 15:12:40 2007
@@ -0,0 +1,401 @@
+// 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.math.ode;
+
+/**
+ * This class implements the common part of all Runge-Kutta-Fehlberg
+ * integrators for Ordinary Differential Equations.
+
+ * <p>These methods are embedded explicit Runge-Kutta methods with two
+ * sets of coefficients allowing to estimate the error, their Butcher
+ * arrays are as follows :
+ * <pre>
+ *    0  |
+ *   c2  | a21
+ *   c3  | a31  a32
+ *   ... |        ...
+ *   cs  | as1  as2  ...  ass-1
+ *       |--------------------------
+ *       |  b1   b2  ...   bs-1  bs
+ *       |  b'1  b'2 ...   b's-1 b's
+ * </pre>
+ * </p>
+
+ * <p>In fact, we rather use the array defined by ej = bj - b'j to
+ * compute directly the error rather than computing two estimates and
+ * then comparing them.</p>
+
+ * <p>Some methods are qualified as <i>fsal</i> (first same as last)
+ * methods. This means the last evaluation of the derivatives in one
+ * step is the same as the first in the next step. Then, this
+ * evaluation can be reused from one step to the next one and the cost
+ * of such a method is really s-1 evaluations despite the method still
+ * has s stages. This behaviour is true only for successful steps, if
+ * the step is rejected after the error estimation phase, no
+ * evaluation is saved. For an <i>fsal</i> method, we have cs = 1 and
+ * asi = bi for all i.</p>
+
+ * @version $Id: RungeKuttaFehlbergIntegrator.java 1705 2006-09-17 19:57:39Z luc $
+ * @author L. Maisonobe
+
+ */
+
+public abstract class RungeKuttaFehlbergIntegrator
+  extends AdaptiveStepsizeIntegrator {
+
+  /** Build a Runge-Kutta integrator with the given Butcher array.
+   * @param fsal indicate that the method is an <i>fsal</i>
+   * @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 external weights for the high order method from Butcher array
+   * @param prototype prototype of the step interpolator to use
+   * @param minStep minimal step (must be positive even for backward
+   * integration), the last step can be smaller than this
+   * @param maxStep maximal step (must be positive even for backward
+   * integration)
+   * @param scalAbsoluteTolerance allowed absolute error
+   * @param scalRelativeTolerance allowed relative error
+   */
+  protected RungeKuttaFehlbergIntegrator(boolean fsal,
+                                         double[] c, double[][] a, double[] b,
+                                         RungeKuttaStepInterpolator prototype,
+                                         double minStep, double maxStep,
+                                         double scalAbsoluteTolerance,
+                                         double scalRelativeTolerance) {
+
+    super(minStep, maxStep, scalAbsoluteTolerance, scalRelativeTolerance);
+
+    this.fsal      = fsal;
+    this.c         = c;
+    this.a         = a;
+    this.b         = b;
+    this.prototype = prototype;
+
+    exp = -1.0 / getOrder();
+
+    this.safety = 0.9;
+
+    // set the default values of the algorithm control parameters
+    setMinReduction(0.2);
+    setMaxGrowth(10.0);
+
+  }
+
+  /** Build a Runge-Kutta integrator with the given Butcher array.
+   * @param fsal indicate that the method is an <i>fsal</i>
+   * @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 external weights for the high order method from Butcher array
+   * @param prototype prototype of the step interpolator to use
+   * @param minStep minimal step (must be positive even for backward
+   * integration), the last step can be smaller than this
+   * @param maxStep maximal step (must be positive even for backward
+   * integration)
+   * @param vecAbsoluteTolerance allowed absolute error
+   * @param vecRelativeTolerance allowed relative error
+   */
+  protected RungeKuttaFehlbergIntegrator(boolean fsal,
+                                         double[] c, double[][] a, double[] b,
+                                         RungeKuttaStepInterpolator prototype,
+                                         double   minStep, double maxStep,
+                                         double[] vecAbsoluteTolerance,
+                                         double[] vecRelativeTolerance) {
+
+    super(minStep, maxStep, vecAbsoluteTolerance, vecRelativeTolerance);
+
+    this.fsal      = fsal;
+    this.c         = c;
+    this.a         = a;
+    this.b         = b;
+    this.prototype = prototype;
+
+    exp = -1.0 / getOrder();
+
+    this.safety = 0.9;
+
+    // set the default values of the algorithm control parameters
+    setMinReduction(0.2);
+    setMaxGrowth(10.0);
+
+  }
+
+  /** Get the name of the method.
+   * @return name of the method
+   */
+  public abstract String getName();
+
+  /** Get the order of the method.
+   * @return order of the method
+   */
+  public abstract int getOrder();
+
+  /** Get the safety factor for stepsize control.
+   * @return safety factor
+   */
+  public double getSafety() {
+    return safety;
+  }
+
+  /** Set the safety factor for stepsize control.
+   * @param safety safety factor
+   */
+  public void setSafety(double safety) {
+    this.safety = safety;
+  }
+
+  public void integrate(FirstOrderDifferentialEquations equations,
+                        double t0, double[] y0,
+                        double t, double[] y)
+  throws DerivativeException, IntegratorException {
+
+    // sanity check
+    if (equations.getDimension() != y0.length) {
+      throw new IntegratorException("dimensions mismatch: ODE problem has dimension {0},"
+                                    + " state vector has dimension {1}",
+                                    new String[] {
+                                      Integer.toString(equations.getDimension()),
+                                      Integer.toString(y0.length)
+                                    });
+    }
+    if (Math.abs(t - t0) <= 1.0e-12 * Math.max(Math.abs(t0), Math.abs(t))) {
+      throw new IntegratorException("too small integration interval: length = {0}",
+                                    new String[] {
+                                      Double.toString(Math.abs(t - t0))
+                                    });
+    }
+    
+    boolean forward = (t > t0);
+
+    // create some internal working arrays
+    int stages = c.length + 1;
+    if (y != y0) {
+      System.arraycopy(y0, 0, y, 0, y0.length);
+    }
+    double[][] yDotK = new double[stages][];
+    for (int i = 0; i < stages; ++i) {
+      yDotK [i] = new double[y0.length];
+    }
+    double[] yTmp = new double[y0.length];
+
+    // set up an interpolator sharing the integrator arrays
+    AbstractStepInterpolator interpolator;
+    if (handler.requiresDenseOutput() || (! switchesHandler.isEmpty())) {
+      RungeKuttaStepInterpolator rki = (RungeKuttaStepInterpolator) prototype.clone();
+      rki.reinitialize(equations, yTmp, yDotK, forward);
+      interpolator = rki;
+    } else {
+      interpolator = new DummyStepInterpolator(yTmp, forward);
+    }
+    interpolator.storeTime(t0);
+
+    double  currentT  = t0;
+    double  hNew      = 0;
+    boolean firstTime = true;
+    boolean lastStep;
+    handler.reset();
+    do {
+
+      interpolator.shift();
+
+      double h     = 0;
+      double error = 0;
+      for (boolean loop = true; loop;) {
+
+        if (firstTime || !fsal) {
+          // first stage
+          equations.computeDerivatives(currentT, y, yDotK[0]);
+        }
+
+        if (firstTime) {
+          double[] scale;
+          if (vecAbsoluteTolerance != null) {
+            scale = vecAbsoluteTolerance;
+          } else {
+            scale = new double[y0.length];
+            for (int i = 0; i < scale.length; ++i) {
+              scale[i] = scalAbsoluteTolerance;
+            }
+          }
+          hNew = initializeStep(equations, forward, getOrder(), scale,
+                                currentT, y, yDotK[0], yTmp, yDotK[1]);
+          firstTime = false;
+        }
+
+        h = hNew;
+
+        // step adjustment near bounds
+        if ((forward && (currentT + h > t))
+            || ((! forward) && (currentT + h < t))) {
+          h = t - currentT;
+        }
+
+        // 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(currentT + 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];
+          }
+          yTmp[j] = y[j] + h * sum;
+        }
+
+        // estimate the error at the end of the step
+        error = estimateError(yDotK, y, yTmp, h);
+        if (error <= 1.0) {
+
+          // Switching functions handling
+          interpolator.storeTime(currentT + h);
+          if (switchesHandler.evaluateStep(interpolator)) {
+            // reject the step to match exactly the next switch time
+            hNew = switchesHandler.getEventTime() - currentT;
+          } else {
+            // accept the step
+            loop = false;
+          }
+
+        } else {
+          // reject the step and attempt to reduce error by stepsize control
+          double factor = Math.min(maxGrowth,
+                                   Math.max(minReduction,
+                                            safety * Math.pow(error, exp)));
+          hNew = filterStep(h * factor, false);
+        }
+
+      }
+
+      // the step has been accepted
+      currentT += h;
+      System.arraycopy(yTmp, 0, y, 0, y0.length);
+      switchesHandler.stepAccepted(currentT, y);
+      if (switchesHandler.stop()) {
+        lastStep = true;
+      } else {
+        lastStep = forward ? (currentT >= t) : (currentT <= t);
+      }
+
+      // provide the step data to the step handler
+      interpolator.storeTime(currentT);
+      handler.handleStep(interpolator, lastStep);
+
+      if (fsal) {
+        // save the last evaluation for the next step
+        System.arraycopy(yDotK[stages - 1], 0, yDotK[0], 0, y0.length);
+      }
+
+      if (switchesHandler.reset(currentT, y) && ! lastStep) {
+        // some switching function has triggered changes that
+        // invalidate the derivatives, we need to recompute them
+        equations.computeDerivatives(currentT, y, yDotK[0]);
+      }
+
+      if (! lastStep) {
+        // stepsize control for next step
+        double  factor     = Math.min(maxGrowth,
+                                      Math.max(minReduction,
+                                               safety * Math.pow(error, exp)));
+        double  scaledH    = h * factor;
+        double  nextT      = currentT + scaledH;
+        boolean nextIsLast = forward ? (nextT >= t) : (nextT <= t);
+        hNew = filterStep(scaledH, nextIsLast);
+      }
+
+    } while (! lastStep);
+
+  }
+
+  /** Get the minimal reduction factor for stepsize control.
+   * @return minimal reduction factor
+   */
+  public double getMinReduction() {
+    return minReduction;
+  }
+
+  /** Set the minimal reduction factor for stepsize control.
+   * @param minReduction minimal reduction factor
+   */
+  public void setMinReduction(double minReduction) {
+    this.minReduction = minReduction;
+  }
+
+  /** Get the maximal growth factor for stepsize control.
+   * @return maximal growth factor
+   */
+  public double getMaxGrowth() {
+    return maxGrowth;
+  }
+
+  /** Set the maximal growth factor for stepsize control.
+   * @param maxGrowth maximal growth factor
+   */
+  public void setMaxGrowth(double maxGrowth) {
+    this.maxGrowth = maxGrowth;
+  }
+
+  /** Compute the error ratio.
+   * @param yDotK derivatives computed during the first stages
+   * @param y0 estimate of the step at the start of the step
+   * @param y1 estimate of the step at the end of the step
+   * @param h  current step
+   * @return error ratio, greater than 1 if step should be rejected
+   */
+  protected abstract double estimateError(double[][] yDotK,
+                                          double[] y0, double[] y1,
+                                          double h);
+
+  /** Indicator for <i>fsal</i> methods. */
+  private boolean fsal;
+
+  /** Time steps from Butcher array (without the first zero). */
+  private double[] c;
+
+  /** Internal weights from Butcher array (without the first empty row). */
+  private double[][] a;
+
+  /** External weights for the high order method from Butcher array. */
+  private double[] b;
+
+  /** Prototype of the step interpolator. */
+  private RungeKuttaStepInterpolator prototype;
+                                         
+  /** Stepsize control exponent. */
+  private double exp;
+
+  /** Safety factor for stepsize control. */
+  private double safety;
+
+  /** Minimal reduction factor for stepsize control. */
+  private double minReduction;
+
+  /** Maximal growth factor for stepsize control. */
+  private double maxGrowth;
+
+}

Propchange: jakarta/commons/proper/math/trunk/src/java/org/apache/commons/math/ode/RungeKuttaFehlbergIntegrator.java
------------------------------------------------------------------------------
    svn:eol-style = native



---------------------------------------------------------------------
To unsubscribe, e-mail: commons-dev-unsubscribe@jakarta.apache.org
For additional commands, e-mail: commons-dev-help@jakarta.apache.org


Mime
View raw message