commons-dev mailing list archives

Site index · List index
Message view « Date » · « Thread »
Top « Date » · « Thread »
From l..@apache.org
Subject svn commit: r512066 [2/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/DormandPrince853Integrator.java
URL: http://svn.apache.org/viewvc/jakarta/commons/proper/math/trunk/src/java/org/apache/commons/math/ode/DormandPrince853Integrator.java?view=auto&rev=512066
==============================================================================
--- jakarta/commons/proper/math/trunk/src/java/org/apache/commons/math/ode/DormandPrince853Integrator.java (added)
+++ jakarta/commons/proper/math/trunk/src/java/org/apache/commons/math/ode/DormandPrince853Integrator.java Mon Feb 26 15:12:40 2007
@@ -0,0 +1,258 @@
+// 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 8(5,3) Dormand-Prince integrator for Ordinary
+ * Differential Equations.
+
+ * <p>This integrator is an embedded Runge-Kutta-Fehlberg integrator
+ * of order 8(5,3) 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 12 functions evaluations per step for integration and 4
+ * evaluations for interpolation. However, since the first
+ * interpolation evaluation is the same as the first integration
+ * evaluation of the next step, we have included it in the integrator
+ * rather than in the interpolator and specified the method was an
+ * <i>fsal</i>. Hence, despite we have 13 stages here, the cost is
+ * really 12 evaluations per step even if no interpolation is done,
+ * and the overcost of interpolation is only 3 evaluations.</p>
+
+ * <p>This method is based on an 8(6) method by Dormand and Prince
+ * (i.e. order 8 for the integration and order 6 for error estimation)
+ * modified by Hairer and Wanner to use a 5th order error estimator
+ * with 3rd order correction. This modification was introduced because
+ * the original method failed in some cases (wrong steps can be
+ * accepted when step size is too large, for example in the
+ * Brusselator problem) and also had <i>severe difficulties when
+ * applied to problems with discontinuities</i>. This modification is
+ * explained in the second edition of the first volume (Nonstiff
+ * Problems) of the reference book by Hairer, Norsett and Wanner:
+ * <i>Solving Ordinary Differential Equations</i> (Springer-Verlag,
+ * ISBN 3-540-56670-8).</p>
+
+ * @version $Id: DormandPrince853Integrator.java 1705 2006-09-17 19:57:39Z luc $
+ * @author L. Maisonobe
+
+ */
+
+public class DormandPrince853Integrator
+  extends RungeKuttaFehlbergIntegrator {
+
+  private static final String methodName = "Dormand-Prince 8 (5, 3)";
+
+  private static final double sqrt6 = Math.sqrt(6.0);
+
+  private static final double[] c = {
+    (12.0 - 2.0 * sqrt6) / 135.0, (6.0 - sqrt6) / 45.0, (6.0 - sqrt6) / 30.0,
+    (6.0 + sqrt6) / 30.0, 1.0/3.0, 1.0/4.0, 4.0/13.0, 127.0/195.0, 3.0/5.0,
+    6.0/7.0, 1.0, 1.0
+  };
+
+  private static final double[][] a = {
+
+    // k2
+    {(12.0 - 2.0 * sqrt6) / 135.0},
+
+    // k3
+    {(6.0 - sqrt6) / 180.0, (6.0 - sqrt6) / 60.0},
+
+    // k4
+    {(6.0 - sqrt6) / 120.0, 0.0, (6.0 - sqrt6) / 40.0},
+
+    // k5
+    {(462.0 + 107.0 * sqrt6) / 3000.0, 0.0,
+     (-402.0 - 197.0 * sqrt6) / 1000.0, (168.0 + 73.0 * sqrt6) / 375.0},
+
+    // k6
+    {1.0 / 27.0, 0.0, 0.0, (16.0 + sqrt6) / 108.0, (16.0 - sqrt6) / 108.0},
+
+    // k7
+    {19.0 / 512.0, 0.0, 0.0, (118.0 + 23.0 * sqrt6) / 1024.0,
+     (118.0 - 23.0 * sqrt6) / 1024.0, -9.0 / 512.0},
+
+    // k8
+    {13772.0 / 371293.0, 0.0, 0.0, (51544.0 + 4784.0 * sqrt6) / 371293.0,
+     (51544.0 - 4784.0 * sqrt6) / 371293.0, -5688.0 / 371293.0, 3072.0 / 371293.0},
+
+    // k9
+    {58656157643.0 / 93983540625.0, 0.0, 0.0,
+     (-1324889724104.0 - 318801444819.0 * sqrt6) / 626556937500.0,
+     (-1324889724104.0 + 318801444819.0 * sqrt6) / 626556937500.0,
+     96044563816.0 / 3480871875.0, 5682451879168.0 / 281950621875.0,
+     -165125654.0 / 3796875.0},
+
+    // k10
+    {8909899.0 / 18653125.0, 0.0, 0.0,
+     (-4521408.0 - 1137963.0 * sqrt6) / 2937500.0,
+     (-4521408.0 + 1137963.0 * sqrt6) / 2937500.0,
+     96663078.0 / 4553125.0, 2107245056.0 / 137915625.0,
+     -4913652016.0 / 147609375.0, -78894270.0 / 3880452869.0},
+
+    // k11
+    {-20401265806.0 / 21769653311.0, 0.0, 0.0,
+     (354216.0 + 94326.0 * sqrt6) / 112847.0,
+     (354216.0 - 94326.0 * sqrt6) / 112847.0,
+     -43306765128.0 / 5313852383.0, -20866708358144.0 / 1126708119789.0,
+     14886003438020.0 / 654632330667.0, 35290686222309375.0 / 14152473387134411.0,
+     -1477884375.0 / 485066827.0},
+
+    // k12
+    {39815761.0 / 17514443.0, 0.0, 0.0,
+     (-3457480.0 - 960905.0 * sqrt6) / 551636.0,
+     (-3457480.0 + 960905.0 * sqrt6) / 551636.0,
+     -844554132.0 / 47026969.0, 8444996352.0 / 302158619.0,
+     -2509602342.0 / 877790785.0, -28388795297996250.0 / 3199510091356783.0,
+     226716250.0 / 18341897.0, 1371316744.0 / 2131383595.0},
+
+    // k13 should be for interpolation only, but since it is the same
+    // stage as the first evaluation of the next step, we perform it
+    // here at no cost by specifying this is an fsal method
+    {104257.0/1920240.0, 0.0, 0.0, 0.0, 0.0, 3399327.0/763840.0,
+     66578432.0/35198415.0, -1674902723.0/288716400.0,
+     54980371265625.0/176692375811392.0, -734375.0/4826304.0,
+     171414593.0/851261400.0, 137909.0/3084480.0}
+
+  };
+
+  private static final double[] b = {
+      104257.0/1920240.0,
+      0.0,
+      0.0,
+      0.0,
+      0.0,
+      3399327.0/763840.0,
+      66578432.0/35198415.0,
+      -1674902723.0/288716400.0,
+      54980371265625.0/176692375811392.0,
+      -734375.0/4826304.0,
+      171414593.0/851261400.0,
+      137909.0/3084480.0,
+      0.0
+  };
+
+  private static final double e1_01 =         116092271.0 / 8848465920.0;
+  private static final double e1_06 =          -1871647.0 / 1527680.0;
+  private static final double e1_07 =         -69799717.0 / 140793660.0;
+  private static final double e1_08 =     1230164450203.0 / 739113984000.0;
+  private static final double e1_09 = -1980813971228885.0 / 5654156025964544.0;
+  private static final double e1_10 =         464500805.0 / 1389975552.0;
+  private static final double e1_11 =     1606764981773.0 / 19613062656000.0;
+  private static final double e1_12 =           -137909.0 / 6168960.0;
+
+  private static final double e2_01 =           -364463.0 / 1920240.0;
+  private static final double e2_06 =           3399327.0 / 763840.0;
+  private static final double e2_07 =          66578432.0 / 35198415.0;
+  private static final double e2_08 =       -1674902723.0 / 288716400.0;
+  private static final double e2_09 =   -74684743568175.0 / 176692375811392.0;
+  private static final double e2_10 =           -734375.0 / 4826304.0;
+  private static final double e2_11 =         171414593.0 / 851261400.0;
+  private static final double e2_12 =             69869.0 / 3084480.0;
+
+  /** Simple constructor.
+   * Build an eighth order Dormand-Prince 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 DormandPrince853Integrator(double minStep, double maxStep,
+                                    double scalAbsoluteTolerance,
+                                    double scalRelativeTolerance) {
+    super(true, c, a, b,
+          new DormandPrince853StepInterpolator(),
+          minStep, maxStep, scalAbsoluteTolerance, scalRelativeTolerance);
+  }
+
+  /** Simple constructor.
+   * Build an eighth order Dormand-Prince 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 DormandPrince853Integrator(double minStep, double maxStep,
+                                    double[] vecAbsoluteTolerance,
+                                    double[] vecRelativeTolerance) {
+    super(true, c, a, b,
+          new DormandPrince853StepInterpolator(),
+          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 8;
+  }
+
+  /** 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 error1 = 0;
+    double error2 = 0;
+
+    for (int j = 0; j < y0.length; ++j) {
+      double errSum1 = e1_01 * yDotK[0][j]  + e1_06 * yDotK[5][j]
+                     + e1_07 * yDotK[6][j]  + e1_08 * yDotK[7][j]
+                     + e1_09 * yDotK[8][j]  + e1_10 * yDotK[9][j]
+                     + e1_11 * yDotK[10][j] + e1_12 * yDotK[11][j];
+      double errSum2 = e2_01 * yDotK[0][j]  + e2_06 * yDotK[5][j]
+                     + e2_07 * yDotK[6][j]  + e2_08 * yDotK[7][j]
+                     + e2_09 * yDotK[8][j]  + e2_10 * yDotK[9][j]
+                     + e2_11 * yDotK[10][j] + e2_12 * yDotK[11][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 ratio1  = errSum1 / tol;
+      error1        += ratio1 * ratio1;
+      double ratio2  = errSum2 / tol;
+      error2        += ratio2 * ratio2;
+    }
+
+    double den = error1 + 0.01 * error2;
+    if (den <= 0.0) {
+      den = 1.0;
+    }
+
+    return Math.abs(h) * error1 / Math.sqrt(y0.length * den);
+
+  }
+
+}

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

Added: jakarta/commons/proper/math/trunk/src/java/org/apache/commons/math/ode/DormandPrince853StepInterpolator.java
URL: http://svn.apache.org/viewvc/jakarta/commons/proper/math/trunk/src/java/org/apache/commons/math/ode/DormandPrince853StepInterpolator.java?view=auto&rev=512066
==============================================================================
--- jakarta/commons/proper/math/trunk/src/java/org/apache/commons/math/ode/DormandPrince853StepInterpolator.java (added)
+++ jakarta/commons/proper/math/trunk/src/java/org/apache/commons/math/ode/DormandPrince853StepInterpolator.java Mon Feb 26 15:12:40 2007
@@ -0,0 +1,413 @@
+// 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.ObjectOutput;
+import java.io.ObjectInput;
+import java.io.IOException;
+
+/**
+ * This class represents an interpolator over the last step during an
+ * ODE integration for the 8(5,3) Dormand-Prince integrator.
+ *
+ * @see DormandPrince853Integrator
+ *
+ * @version $Id: DormandPrince853StepInterpolator.java 1705 2006-09-17 19:57:39Z luc $
+ * @author L. Maisonobe
+ *
+ */
+
+class DormandPrince853StepInterpolator
+  extends RungeKuttaStepInterpolator {
+
+  /** Simple constructor.
+   * This constructor builds an instance that is not usable yet, the
+   * {@link #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 DormandPrince853StepInterpolator() {
+    super();
+    yDotKLast = null;
+    yTmp      = null;
+    v         = null;
+    vectorsInitialized = false;
+  }
+
+  /** 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 DormandPrince853StepInterpolator(DormandPrince853StepInterpolator interpolator) {
+
+    super(interpolator);
+
+    if (interpolator.currentState == null) {
+
+      yDotKLast = null;
+      v         = null;
+      vectorsInitialized = false;
+
+    } else {
+
+      int dimension = interpolator.currentState.length;
+
+      yDotKLast    = new double[3][];
+      for (int k = 0; k < yDotKLast.length; ++k) {
+        yDotKLast[k] = new double[dimension];
+        System.arraycopy(interpolator.yDotKLast[k], 0, yDotKLast[k], 0,
+                         dimension);
+      }
+
+      v = new double[7][];
+      for (int k = 0; k < v.length; ++k) {
+        v[k] = new double[dimension];
+        System.arraycopy(interpolator.v[k], 0, v[k], 0, dimension);
+      }
+
+      vectorsInitialized = interpolator.vectorsInitialized;
+
+    }
+
+    // the step has been finalized, we don't need this anymore
+    yTmp = null;
+
+  }
+
+  /**
+   * 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 DormandPrince853StepInterpolator(this);
+  }
+
+  /** Reinitialize the instance
+   * Some Runge-Kutta-Fehlberg integrators need fewer functions
+   * evaluations than their counterpart step interpolators. So the
+   * interpolator should perform the last evaluations they need by
+   * themselves. The {@link RungeKuttaFehlbergIntegrator
+   * RungeKuttaFehlbergIntegrator} abstract class calls this method in
+   * order to let the step interpolator perform the evaluations it
+   * needs. These evaluations will be performed during the call to
+   * <code>doFinalize</code> if any, i.e. only if the step handler
+   * either calls the {@link AbstractStepInterpolator#finalizeStep
+   * finalizeStep} method or the {@link
+   * AbstractStepInterpolator#getInterpolatedState getInterpolatedState}
+   * method (for an interpolator which needs a finalization) or if it clones
+   * the step interpolator.
+   * @param equations set of differential equations being integrated
+   * @param y reference to the integrator array holding the state at
+   * the end of the step
+   * @param yDotK reference to the integrator array holding all the
+   * intermediate slopes
+   * @param forward integration direction indicator
+   */
+  public void reinitialize(FirstOrderDifferentialEquations equations,
+                           double[] y, double[][] yDotK, boolean forward) {
+
+    super.reinitialize(equations, y, yDotK, forward);
+
+    int dimension = currentState.length;
+
+    yDotKLast = new double[3][];
+    for (int k = 0; k < yDotKLast.length; ++k) {
+      yDotKLast[k] = new double[dimension];
+    }
+
+    yTmp = new double[dimension];
+
+    v = new double[7][];
+    for (int k = 0; k < v.length; ++k) {
+      v[k]  = new double[dimension];
+    }
+
+    vectorsInitialized = false;
+
+  }
+
+  /** Store the current step time.
+   * @param t current time
+   */
+  public void storeTime(double t) {
+    super.storeTime(t);
+    vectorsInitialized = false;
+  }
+
+  /** 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 {
+
+    if (! vectorsInitialized) {
+
+      if (v == null) {
+        v = new double[7][];
+        for (int k = 0; k < 7; ++k) {
+          v[k] = new double[interpolatedState.length];
+        }
+      }
+
+      // perform the last evaluations if they have not been done yet
+      finalizeStep();
+
+      // compute the interpolation vectors for this time step
+      for (int i = 0; i < interpolatedState.length; ++i) {
+        v[0][i] = h * (b_01 * yDotK[0][i]  + b_06 * yDotK[5][i] + b_07 * yDotK[6][i]
+                     + b_08 * yDotK[7][i]  + b_09 * yDotK[8][i] + b_10 * yDotK[9][i]
+                     + b_11 * yDotK[10][i] + b_12 * yDotK[11][i]);
+        v[1][i] = h * yDotK[0][i] - v[0][i];
+        v[2][i] = v[0][i] - v[1][i] - h * yDotK[12][i];
+        for (int k = 0; k < d.length; ++k) {
+          v[k+3][i] = h * (d[k][0] * yDotK[0][i]  + d[k][1] * yDotK[5][i]  + d[k][2] * yDotK[6][i]
+                         + d[k][3] * yDotK[7][i]  + d[k][4] * yDotK[8][i]  + d[k][5] * yDotK[9][i]
+                         + d[k][6] * yDotK[10][i] + d[k][7] * yDotK[11][i] + d[k][8] * yDotK[12][i]
+                         + d[k][9]  * yDotKLast[0][i]
+                         + d[k][10] * yDotKLast[1][i]
+                         + d[k][11] * yDotKLast[2][i]);
+        }
+      }
+
+      vectorsInitialized = true;
+
+    }
+
+    double eta = oneMinusThetaH / h;
+
+    for (int i = 0; i < interpolatedState.length; ++i) {
+      interpolatedState[i] = currentState[i]
+                           -   eta * (v[0][i]
+                           - theta * (v[1][i]
+                           + theta * (v[2][i]
+                           +   eta * (v[3][i]
+                           + theta * (v[4][i]
+                           +   eta * (v[5][i]
+                           + theta * (v[6][i])))))));
+    }
+
+  }
+ 
+  /**
+   * Really finalize the step.
+   * Perform the last 3 functions evaluations (k14, k15, k16)
+   * @throws DerivativeException this exception is propagated to the caller if the
+   * underlying user function triggers one
+   */
+  protected void doFinalize()
+    throws DerivativeException {
+
+    double s;
+
+    // k14
+    for (int j = 0; j < currentState.length; ++j) {
+      s = k14_01 * yDotK[0][j]  + k14_06 * yDotK[5][j]  + k14_07 * yDotK[6][j]
+        + k14_08 * yDotK[7][j]  + k14_09 * yDotK[8][j]  + k14_10 * yDotK[9][j]
+        + k14_11 * yDotK[10][j] + k14_12 * yDotK[11][j] + k14_13 * yDotK[12][j];
+      yTmp[j] = currentState[j] + h * s;
+    }
+    equations.computeDerivatives(previousTime + c14 * h, yTmp, yDotKLast[0]);
+
+    // k15
+    for (int j = 0; j < currentState.length; ++j) {
+     s = k15_01 * yDotK[0][j]  + k15_06 * yDotK[5][j]  + k15_07 * yDotK[6][j]
+       + k15_08 * yDotK[7][j]  + k15_09 * yDotK[8][j]  + k15_10 * yDotK[9][j]
+       + k15_11 * yDotK[10][j] + k15_12 * yDotK[11][j] + k15_13 * yDotK[12][j]
+       + k15_14 * yDotKLast[0][j];
+     yTmp[j] = currentState[j] + h * s;
+    }
+    equations.computeDerivatives(previousTime + c15 * h, yTmp, yDotKLast[1]);
+
+    // k16
+    for (int j = 0; j < currentState.length; ++j) {
+      s = k16_01 * yDotK[0][j]  + k16_06 * yDotK[5][j]  + k16_07 * yDotK[6][j]
+        + k16_08 * yDotK[7][j]  + k16_09 * yDotK[8][j]  + k16_10 * yDotK[9][j]
+        + k16_11 * yDotK[10][j] + k16_12 * yDotK[11][j] + k16_13 * yDotK[12][j]
+        + k16_14 * yDotKLast[0][j] +  k16_15 * yDotKLast[1][j];
+      yTmp[j] = currentState[j] + h * s;
+    }
+    equations.computeDerivatives(previousTime + c16 * h, yTmp, yDotKLast[2]);
+
+  }
+
+  /** 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 {
+
+    try {
+      // save the local attributes
+      finalizeStep();
+    } catch (DerivativeException e) {
+      throw new IOException(e.getMessage());
+    }
+    out.writeInt(currentState.length);
+    for (int i = 0; i < currentState.length; ++i) {
+      out.writeDouble(yDotKLast[0][i]);
+      out.writeDouble(yDotKLast[1][i]);
+      out.writeDouble(yDotKLast[2][i]);
+    }
+
+    // save the state of the base class
+    super.writeExternal(out);
+
+  }
+
+  /** 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 local attributes
+    yDotKLast = new double[3][];
+    int dimension = in.readInt();
+    yDotKLast[0] = new double[dimension];
+    yDotKLast[1] = new double[dimension];
+    yDotKLast[2] = new double[dimension];
+
+    for (int i = 0; i < dimension; ++i) {
+      yDotKLast[0][i] = in.readDouble();
+      yDotKLast[1][i] = in.readDouble();
+      yDotKLast[2][i] = in.readDouble();
+    }
+
+    // read the base state
+    super.readExternal(in);
+
+  }
+
+  /** Last evaluations. */
+  private double[][] yDotKLast;
+
+  /** Temporary state vector. */
+  private double[] yTmp;
+
+  /** Vectors for interpolation. */
+  private double[][] v;
+
+  /** Initialization indicator for the interpolation vectors. */
+  private boolean vectorsInitialized;
+
+  // external weights of the integrator,
+  // note that b_02 through b_05 are null
+  private static double b_01 =         104257.0 / 1920240.0;
+  private static double b_06 =        3399327.0 / 763840.0;
+  private static double b_07 =       66578432.0 / 35198415.0;
+  private static double b_08 =    -1674902723.0 / 288716400.0;
+  private static double b_09 = 54980371265625.0 / 176692375811392.0;
+  private static double b_10 =        -734375.0 / 4826304.0;
+  private static double b_11 =      171414593.0 / 851261400.0;
+  private static double b_12 =         137909.0 / 3084480.0;
+
+  // k14 for interpolation only
+  private static double c14    = 1.0 / 10.0;
+
+  private static double k14_01 =       13481885573.0 / 240030000000.0      - b_01;
+  private static double k14_06 =                 0.0                       - b_06;
+  private static double k14_07 =      139418837528.0 / 549975234375.0      - b_07;
+  private static double k14_08 =   -11108320068443.0 / 45111937500000.0    - b_08;
+  private static double k14_09 = -1769651421925959.0 / 14249385146080000.0 - b_09;
+  private static double k14_10 =          57799439.0 / 377055000.0         - b_10;
+  private static double k14_11 =      793322643029.0 / 96734250000000.0    - b_11;
+  private static double k14_12 =        1458939311.0 / 192780000000.0      - b_12;
+  private static double k14_13 =             -4149.0 / 500000.0;
+
+  // k15 for interpolation only
+  private static double c15    = 1.0 / 5.0;
+
+  private static double k15_01 =     1595561272731.0 / 50120273500000.0    - b_01;
+  private static double k15_06 =      975183916491.0 / 34457688031250.0    - b_06;
+  private static double k15_07 =    38492013932672.0 / 718912673015625.0   - b_07;
+  private static double k15_08 = -1114881286517557.0 / 20298710767500000.0 - b_08;
+  private static double k15_09 =                 0.0                       - b_09;
+  private static double k15_10 =                 0.0                       - b_10;
+  private static double k15_11 =    -2538710946863.0 / 23431227861250000.0 - b_11;
+  private static double k15_12 =        8824659001.0 / 23066716781250.0    - b_12;
+  private static double k15_13 =      -11518334563.0 / 33831184612500.0;
+  private static double k15_14 =        1912306948.0 / 13532473845.0;
+
+  // k16 for interpolation only
+  private static double c16    = 7.0 / 9.0;
+
+  private static double k16_01 =      -13613986967.0 / 31741908048.0       - b_01;
+  private static double k16_06 =       -4755612631.0 / 1012344804.0        - b_06;
+  private static double k16_07 =    42939257944576.0 / 5588559685701.0     - b_07;
+  private static double k16_08 =    77881972900277.0 / 19140370552944.0    - b_08;
+  private static double k16_09 =    22719829234375.0 / 63689648654052.0    - b_09;
+  private static double k16_10 =                 0.0                       - b_10;
+  private static double k16_11 =                 0.0                       - b_11;
+  private static double k16_12 =                 0.0                       - b_12;
+  private static double k16_13 =       -1199007803.0 / 857031517296.0;
+  private static double k16_14 =      157882067000.0 / 53564469831.0;
+  private static double k16_15 =     -290468882375.0 / 31741908048.0;
+
+  // interpolation weights
+  // (beware that only the non-null values are in the table)
+  private static double[][] d = {
+
+    {        -17751989329.0 / 2106076560.0,               4272954039.0 / 7539864640.0,
+            -118476319744.0 / 38604839385.0,            755123450731.0 / 316657731600.0,
+      3692384461234828125.0 / 1744130441634250432.0,     -4612609375.0 / 5293382976.0,
+            2091772278379.0 / 933644586600.0,             2136624137.0 / 3382989120.0,
+                  -126493.0 / 1421424.0,                    98350000.0 / 5419179.0,
+                -18878125.0 / 2053168.0,                 -1944542619.0 / 438351368.0},
+
+    {         32941697297.0 / 3159114840.0,             456696183123.0 / 1884966160.0,
+           19132610714624.0 / 115814518155.0,       -177904688592943.0 / 474986597400.0,
+     -4821139941836765625.0 / 218016305204281304.0,      30702015625.0 / 3970037232.0,
+          -85916079474274.0 / 2800933759800.0,           -5919468007.0 / 634310460.0,
+                  2479159.0 / 157936.0,                    -18750000.0 / 602131.0,
+                -19203125.0 / 2053168.0,                 15700361463.0 / 438351368.0},
+
+    {         12627015655.0 / 631822968.0,              -72955222965.0 / 188496616.0,
+          -13145744952320.0 / 69488710893.0,          30084216194513.0 / 56998391688.0,
+      -296858761006640625.0 / 25648977082856624.0,         569140625.0 / 82709109.0,
+             -18684190637.0 / 18672891732.0,                69644045.0 / 89549712.0,
+                -11847025.0 / 4264272.0,                  -978650000.0 / 16257537.0,
+                519371875.0 / 6159504.0,                  5256837225.0 / 438351368.0},
+
+    {          -450944925.0 / 17550638.0,               -14532122925.0 / 94248308.0,
+            -595876966400.0 / 2573655959.0,             188748653015.0 / 527762886.0,
+      2545485458115234375.0 / 27252038150535163.0,       -1376953125.0 / 36759604.0,
+              53995596795.0 / 518691437.0,                 210311225.0 / 7047894.0,
+                 -1718875.0 / 39484.0,                      58000000.0 / 602131.0,
+                 -1546875.0 / 39484.0,                   -1262172375.0 / 8429834.0}
+
+  };
+
+  private static final long serialVersionUID = 4165537490327432186L;
+
+}

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

Added: jakarta/commons/proper/math/trunk/src/java/org/apache/commons/math/ode/DummyStepHandler.java
URL: http://svn.apache.org/viewvc/jakarta/commons/proper/math/trunk/src/java/org/apache/commons/math/ode/DummyStepHandler.java?view=auto&rev=512066
==============================================================================
--- jakarta/commons/proper/math/trunk/src/java/org/apache/commons/math/ode/DummyStepHandler.java (added)
+++ jakarta/commons/proper/math/trunk/src/java/org/apache/commons/math/ode/DummyStepHandler.java Mon Feb 26 15:12:40 2007
@@ -0,0 +1,92 @@
+// 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 is a step handler that do nothing.
+
+ * <p>This class is provided as a convenience for users who are only
+ * interested in the final state of an integration and not in the
+ * intermediate steps. Its handleStep method does nothing.</p>
+
+ * <p>Since this class has no internal state, it is implemented using
+ * the Singleton design pattern. This means that only one instance is
+ * ever created, which can be retrieved using the getInstance
+ * method. This explains why there is no public constructor.</p>
+
+ * @see StepHandler
+
+ * @version $Id: DummyStepHandler.java 1705 2006-09-17 19:57:39Z luc $
+ * @author L. Maisonobe
+
+ */
+
+public class DummyStepHandler
+  implements StepHandler {
+
+  /** Private constructor.
+   * The constructor is private to prevent users from creating
+   * instances (Singleton design-pattern).
+   */
+  private DummyStepHandler() {
+  }
+
+  /** Get the only instance.
+   * @return the only instance
+   */
+  public static DummyStepHandler getInstance() {
+    if (instance == null) {
+      instance = new DummyStepHandler();
+    }
+    return instance;
+  }
+
+  /** Determines whether this handler needs dense output.
+   * Since this handler does nothing, it does not require dense output.
+   * @return always false
+   */
+  public boolean requiresDenseOutput() {
+    return false;
+  }
+
+  /** Reset the step handler.
+   * Initialize the internal data as required before the first step is
+   * handled.
+   */
+  public void reset() {
+  }
+
+  /**
+   * Handle the last accepted step.
+   * This method does nothing in this class.
+   * @param interpolator interpolator for the last accepted step. For
+   * efficiency purposes, the various integrators reuse the same
+   * object on each call, so if the instance wants to keep it across
+   * all calls (for example to provide at the end of the integration a
+   * continuous model valid throughout the integration range), it
+   * should build a local copy using the clone method and store this
+   * copy.
+   * @param isLast true if the step is the last one
+   */
+  public void handleStep(StepInterpolator interpolator, boolean isLast) {
+  }
+
+  /** The only instance. */
+  private static DummyStepHandler instance = null;
+
+}

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

Added: jakarta/commons/proper/math/trunk/src/java/org/apache/commons/math/ode/DummyStepInterpolator.java
URL: http://svn.apache.org/viewvc/jakarta/commons/proper/math/trunk/src/java/org/apache/commons/math/ode/DummyStepInterpolator.java?view=auto&rev=512066
==============================================================================
--- jakarta/commons/proper/math/trunk/src/java/org/apache/commons/math/ode/DummyStepInterpolator.java (added)
+++ jakarta/commons/proper/math/trunk/src/java/org/apache/commons/math/ode/DummyStepInterpolator.java Mon Feb 26 15:12:40 2007
@@ -0,0 +1,123 @@
+// 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 is a step interpolator that does nothing.
+ *
+ * <p>This class is used when the {@link StepHandler "step handler"}
+ * set up by the user does not need step interpolation. It does not
+ * recompute the state when {@link AbstractStepInterpolator#setInterpolatedTime
+ * setInterpolatedTime} is called. This implies the interpolated state
+ * is always the state at the end of the current step.</p>
+ *
+ * @see StepHandler
+ *
+ * @version $Id: DummyStepInterpolator.java 1705 2006-09-17 19:57:39Z luc $
+ * @author L. Maisonobe
+ *
+ */
+
+public class DummyStepInterpolator
+  extends AbstractStepInterpolator {
+
+  /** 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. As an example, the {@link
+   * RungeKuttaFehlbergIntegrator} uses the prototyping design pattern
+   * to create the step interpolators by cloning an uninitialized
+   * model and latter initializing the copy.
+   */
+  public DummyStepInterpolator() {
+    super();
+  }
+
+  /** Simple constructor.
+   * @param y reference to the integrator array holding the state at
+   * the end of the step
+   * @param forward integration direction indicator
+   */
+  protected DummyStepInterpolator(double[] y, boolean forward) {
+    super(y, forward);
+  }
+
+  /** Copy constructor.
+
+   * <p>The copied interpolator should have been finalized before the
+   * copy, otherwise the copy will not be able to perform correctly
+   * any interpolation and will throw a {@link NullPointerException}
+   * later. Since we don't want this constructor to throw the
+   * exceptions finalization may involve and since we don't want this
+   * method to modify the state of the copied interpolator,
+   * finalization is <strong>not</strong> done automatically, it
+   * remains under user control.</p>
+
+   * <p>The copy is a deep copy: its arrays are separated from the
+   * original arrays of the instance.</p>
+
+   * @param interpolator interpolator to copy from.
+
+   */
+  protected DummyStepInterpolator(DummyStepInterpolator interpolator) {
+    super(interpolator);
+  }
+
+  /** Compute the state at the interpolated time.
+   * In this class, this method does nothing: the interpolated state
+   * is always the state at the end of the current step.
+   * @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 {
+  }
+    
+  public void writeExternal(ObjectOutput out)
+    throws IOException {
+    // save the state of the base class
+    writeBaseExternal(out);
+  }
+
+  public void readExternal(ObjectInput in)
+    throws IOException {
+
+    // read the base class 
+    double t = readBaseExternal(in);
+
+    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 = 1708010296707839488L;
+
+}

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

Added: jakarta/commons/proper/math/trunk/src/java/org/apache/commons/math/ode/EulerIntegrator.java
URL: http://svn.apache.org/viewvc/jakarta/commons/proper/math/trunk/src/java/org/apache/commons/math/ode/EulerIntegrator.java?view=auto&rev=512066
==============================================================================
--- jakarta/commons/proper/math/trunk/src/java/org/apache/commons/math/ode/EulerIntegrator.java (added)
+++ jakarta/commons/proper/math/trunk/src/java/org/apache/commons/math/ode/EulerIntegrator.java Mon Feb 26 15:12:40 2007
@@ -0,0 +1,80 @@
+// 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 simple Euler integrator for Ordinary
+ * Differential Equations.
+
+ * <p>The Euler algorithm is the simplest one that can be used to
+ * integrate ordinary differential equations. It is a simple inversion
+ * of the forward difference expression :
+ * <code>f'=(f(t+h)-f(t))/h</code> which leads to
+ * <code>f(t+h)=f(t)+hf'</code>. The interpolation scheme used for
+ * dense output is the linear scheme already used for integration.</p>
+
+ * <p>This algorithm looks cheap because it needs only one function
+ * evaluation per step. However, as it uses linear estimates, it needs
+ * very small steps to achieve high accuracy, and small steps lead to
+ * numerical errors and instabilities.</p>
+
+ * <p>This algorithm is almost never used and has been included in
+ * this package only as a comparison reference for more useful
+ * integrators.</p>
+
+ * @see MidpointIntegrator
+ * @see ClassicalRungeKuttaIntegrator
+ * @see GillIntegrator
+ * @see ThreeEighthesIntegrator
+
+ * @version $Id: EulerIntegrator.java 1705 2006-09-17 19:57:39Z luc $
+ * @author L. Maisonobe
+
+ */
+
+public class EulerIntegrator
+  extends RungeKuttaIntegrator {
+
+  private static final String methodName = "Euler";
+
+  private static final double[] c = {
+  };
+
+  private static final double[][] a = {
+  };
+
+  private static final double[] b = {
+    1.0
+  };
+
+  /** Simple constructor.
+   * Build an Euler integrator with the given step.
+   * @param step integration step
+   */
+  public EulerIntegrator(double step) {
+    super(false, c, a, b, new EulerStepInterpolator(), 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/EulerIntegrator.java
------------------------------------------------------------------------------
    svn:eol-style = native

Added: jakarta/commons/proper/math/trunk/src/java/org/apache/commons/math/ode/EulerStepInterpolator.java
URL: http://svn.apache.org/viewvc/jakarta/commons/proper/math/trunk/src/java/org/apache/commons/math/ode/EulerStepInterpolator.java?view=auto&rev=512066
==============================================================================
--- jakarta/commons/proper/math/trunk/src/java/org/apache/commons/math/ode/EulerStepInterpolator.java (added)
+++ jakarta/commons/proper/math/trunk/src/java/org/apache/commons/math/ode/EulerStepInterpolator.java Mon Feb 26 15:12:40 2007
@@ -0,0 +1,97 @@
+// 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 linear interpolator for step.
+
+ * <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 y'
+ * </pre>
+
+ * where theta belongs to [0 ; 1] and where y' is the evaluation of
+ * the derivatives already computed during the step.</p>
+
+ * @see EulerIntegrator
+
+ * @version $Id: EulerStepInterpolator.java 1705 2006-09-17 19:57:39Z luc $
+ * @author L. Maisonobe
+
+ */
+
+class EulerStepInterpolator
+  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 EulerStepInterpolator() {
+  }
+
+  /** 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 EulerStepInterpolator(EulerStepInterpolator 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 EulerStepInterpolator(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 {
+
+    for (int i = 0; i < interpolatedState.length; ++i) {
+      interpolatedState[i] = currentState[i] - oneMinusThetaH * yDotK[0][i];
+    }
+
+  }
+
+  private static final long serialVersionUID = -7179861704951334960L;
+
+}

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

Added: jakarta/commons/proper/math/trunk/src/java/org/apache/commons/math/ode/FirstOrderConverter.java
URL: http://svn.apache.org/viewvc/jakarta/commons/proper/math/trunk/src/java/org/apache/commons/math/ode/FirstOrderConverter.java?view=auto&rev=512066
==============================================================================
--- jakarta/commons/proper/math/trunk/src/java/org/apache/commons/math/ode/FirstOrderConverter.java (added)
+++ jakarta/commons/proper/math/trunk/src/java/org/apache/commons/math/ode/FirstOrderConverter.java Mon Feb 26 15:12:40 2007
@@ -0,0 +1,108 @@
+// 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 converts second order differential equations to first
+ * order ones.
+
+ * <p>This class is a wrapper around a {@link
+ * SecondOrderDifferentialEquations} which allow to use a {@link
+ * FirstOrderIntegrator} to integrate it.</p>
+
+ * <p>The transformation is done by changing the n dimension state
+ * vector to a 2n dimension vector, where the first n components are
+ * the initial state variables and the n last components are their
+ * first time derivative. The first time derivative of this state
+ * vector then really contains both the first and second time
+ * derivative of the initial state vector, which can be handled by the
+ * underlying second order equations set.</p>
+
+ * <p>One should be aware that the data is duplicated during the
+ * transformation process and that for each call to {@link
+ * #computeDerivatives computeDerivatives}, this wrapper does copy 4n
+ * scalars : 2n before the call to {@link
+ * SecondOrderDifferentialEquations#computeSecondDerivatives
+ * computeSecondDerivatives} in order to dispatch the y state vector
+ * into z and zDot, and 2n after the call to gather zDot and zDDot
+ * into yDot. Since the underlying problem by itself perhaps also
+ * needs to copy data and dispatch the arrays into domain objects,
+ * this has an impact on both memory and CPU usage. The only way to
+ * avoid this duplication is to perform the transformation at the
+ * problem level, i.e. to implement the problem as a first order one
+ * and then avoid using this class.</p>
+
+ * @see FirstOrderIntegrator
+ * @see FirstOrderDifferentialEquations
+ * @see SecondOrderDifferentialEquations
+
+ * @version $Id: FirstOrderConverter.java 1705 2006-09-17 19:57:39Z luc $
+ * @author L. Maisonobe
+
+ */
+
+public class FirstOrderConverter
+  implements FirstOrderDifferentialEquations {
+
+  /** Simple constructor.
+   * Build a converter around a second order equations set.
+   * @param equations second order equations set to convert
+   */
+  public FirstOrderConverter (SecondOrderDifferentialEquations equations) {
+      this.equations = equations;
+      dimension      = equations.getDimension();
+      z              = new double[dimension];
+      zDot           = new double[dimension];
+      zDDot          = new double[dimension];
+  }
+
+  public int getDimension() {
+    return 2 * dimension;
+  }
+
+  public void computeDerivatives(double t, double[] y, double[] yDot)
+  throws DerivativeException {
+
+    // split the state vector in two
+    System.arraycopy(y, 0,         z,    0, dimension);
+    System.arraycopy(y, dimension, zDot, 0, dimension);
+
+    // apply the underlying equations set
+    equations.computeSecondDerivatives(t, z, zDot, zDDot);
+
+    // build the result state derivative
+    System.arraycopy(zDot,  0, yDot, 0,         dimension);
+    System.arraycopy(zDDot, 0, yDot, dimension, dimension);
+    
+  }
+
+  /** Underlying second order equations set. */
+  private SecondOrderDifferentialEquations equations;
+
+  /** second order problem dimension. */
+  private int dimension;
+
+  /** state vector. */
+  private double[] z;
+
+  /** first time derivative of the state vector. */
+  private double[] zDot;
+
+  /** second time derivative of the state vector. */
+  private double[] zDDot;
+
+}

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

Added: jakarta/commons/proper/math/trunk/src/java/org/apache/commons/math/ode/FirstOrderDifferentialEquations.java
URL: http://svn.apache.org/viewvc/jakarta/commons/proper/math/trunk/src/java/org/apache/commons/math/ode/FirstOrderDifferentialEquations.java?view=auto&rev=512066
==============================================================================
--- jakarta/commons/proper/math/trunk/src/java/org/apache/commons/math/ode/FirstOrderDifferentialEquations.java (added)
+++ jakarta/commons/proper/math/trunk/src/java/org/apache/commons/math/ode/FirstOrderDifferentialEquations.java Mon Feb 26 15:12:40 2007
@@ -0,0 +1,66 @@
+// 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 interface represents a first order differential equations set.
+ *
+ * <p>This interface should be implemented by all real first order
+ * differential equation problems before they can be handled by the
+ * integrators {@link FirstOrderIntegrator#integrate} method.</p>
+ *
+ * <p>A first order differential equations problem, as seen by an
+ * integrator is the time derivative <code>dY/dt</code> of a state
+ * vector <code>Y</code>, both being one dimensional arrays. From the
+ * integrator point of view, this derivative depends only on the
+ * current time <code>t</code> and on the state vector
+ * <code>Y</code>.</p>
+ *
+ * <p>For real problems, the derivative depends also on parameters
+ * that do not belong to the state vector (dynamical model constants
+ * for example). These constants are completely outside of the scope
+ * of this interface, the classes that implement it are allowed to
+ * handle them as they want.</p>
+ *
+ * @see FirstOrderIntegrator
+ * @see FirstOrderConverter
+ * @see SecondOrderDifferentialEquations
+ * @see org.spaceroots.mantissa.utilities.ArraySliceMappable
+ *
+ * @version $Id: FirstOrderDifferentialEquations.java 1705 2006-09-17 19:57:39Z luc $
+ * @author L. Maisonobe
+ *
+ */
+
+public interface FirstOrderDifferentialEquations {
+    
+    /** Get the dimension of the problem.
+     * @return dimension of the problem
+     */
+    public int getDimension();
+    
+    /** Get the current time derivative of the state vector.
+     * @param t current value of the independant <I>time</I> variable
+     * @param y array containing the current value of the state vector
+     * @param yDot placeholder array where to put the time derivative of the state vector
+     * @throws DerivativeException this exception is propagated to the caller if the
+     * underlying user function triggers one
+     */
+    public void computeDerivatives(double t, double[] y, double[] yDot)
+    throws DerivativeException;
+    
+}

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

Added: jakarta/commons/proper/math/trunk/src/java/org/apache/commons/math/ode/FirstOrderIntegrator.java
URL: http://svn.apache.org/viewvc/jakarta/commons/proper/math/trunk/src/java/org/apache/commons/math/ode/FirstOrderIntegrator.java?view=auto&rev=512066
==============================================================================
--- jakarta/commons/proper/math/trunk/src/java/org/apache/commons/math/ode/FirstOrderIntegrator.java (added)
+++ jakarta/commons/proper/math/trunk/src/java/org/apache/commons/math/ode/FirstOrderIntegrator.java Mon Feb 26 15:12:40 2007
@@ -0,0 +1,84 @@
+// 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 interface represents a first order integrator for
+ * differential equations.
+
+ * <p>The classes which are devoted to solve first order differential
+ * equations should implement this interface. The problems which can
+ * be handled should implement the {@link
+ * FirstOrderDifferentialEquations} interface.</p>
+
+ * @see FirstOrderDifferentialEquations
+ * @see StepHandler
+ * @see SwitchingFunction
+
+ * @version $Id: FirstOrderIntegrator.java 1705 2006-09-17 19:57:39Z luc $
+ * @author L. Maisonobe
+
+ */
+
+public interface FirstOrderIntegrator {
+
+  /** Get the name of the method.
+   * @return name of the method
+   */
+  public String getName();
+
+  /** 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);
+
+  /** Get the step handler for this integrator.
+   * @return the step handler for this integrator
+   */
+  public StepHandler getStepHandler();
+
+  /** 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);
+
+  /** Integrate the differential equations up to the given time
+   * @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 thant <code>t0</code> for backward integration)
+   * @param y placeholder where to put the state vector at each successful
+   *  step (and hence at the end of integration), can be the same object as y0
+   * @throws IntegratorException if the integrator cannot perform integration
+   * @throws DerivativeException this exception is propagated to the caller if
+   * the underlying user function triggers one
+   */
+  public void integrate (FirstOrderDifferentialEquations equations,
+                         double t0, double[] y0,
+                         double t, double[] y)
+    throws DerivativeException, IntegratorException;
+
+}

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

Added: jakarta/commons/proper/math/trunk/src/java/org/apache/commons/math/ode/FixedStepHandler.java
URL: http://svn.apache.org/viewvc/jakarta/commons/proper/math/trunk/src/java/org/apache/commons/math/ode/FixedStepHandler.java?view=auto&rev=512066
==============================================================================
--- jakarta/commons/proper/math/trunk/src/java/org/apache/commons/math/ode/FixedStepHandler.java (added)
+++ jakarta/commons/proper/math/trunk/src/java/org/apache/commons/math/ode/FixedStepHandler.java Mon Feb 26 15:12:40 2007
@@ -0,0 +1,57 @@
+// 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 interface represents a handler that should be called after
+ * each successful fixed step.
+
+ * <p>This interface should be implemented by anyone who is interested
+ * in getting the solution of an ordinary differential equation at
+ * fixed time steps. Objects implementing this interface should be
+ * wrapped within an instance of {@link StepNormalizer} that itself
+ * is used as the general {@link StepHandler} by the integrator. The
+ * {@link StepNormalizer} object is called according to the integrator
+ * internal algorithms and it calls objects implementing this
+ * interface as necessary at fixed time steps.</p>
+
+ * @see StepHandler
+ * @see StepNormalizer
+
+ * @version $Id: FixedStepHandler.java 1705 2006-09-17 19:57:39Z luc $
+ * @author L. Maisonobe
+
+ */
+
+public interface FixedStepHandler {
+
+  /**
+   * Handle the last accepted step
+   * @param t time of the current step
+
+   * @param y state vector at t. For efficiency purposes, the {@link
+   * StepNormalizer} class reuse the same array on each call, so if
+   * the instance wants to keep it across all calls (for example to
+   * provide at the end of the integration a complete array of all
+   * steps), it should build a local copy store this copy.
+
+   * @param isLast true if the step is the last one
+   */
+  public void handleStep(double t, double[] y, boolean isLast);
+    
+}

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

Added: jakarta/commons/proper/math/trunk/src/java/org/apache/commons/math/ode/GillIntegrator.java
URL: http://svn.apache.org/viewvc/jakarta/commons/proper/math/trunk/src/java/org/apache/commons/math/ode/GillIntegrator.java?view=auto&rev=512066
==============================================================================
--- jakarta/commons/proper/math/trunk/src/java/org/apache/commons/math/ode/GillIntegrator.java (added)
+++ jakarta/commons/proper/math/trunk/src/java/org/apache/commons/math/ode/GillIntegrator.java Mon Feb 26 15:12:40 2007
@@ -0,0 +1,82 @@
+// 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 Gill fourth 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       0      0
+ *   1/2 |   1/2       0       0      0
+ *   1/2 | (q-1)/2  (2-q)/2    0      0
+ *    1  |    0       -q/2  (2+q)/2   0
+ *       |-------------------------------
+ *       |   1/6    (2-q)/6 (2+q)/6  1/6
+ * </pre>
+ * where q = sqrt(2)</p>
+
+ * @see EulerIntegrator
+ * @see ClassicalRungeKuttaIntegrator
+ * @see MidpointIntegrator
+ * @see ThreeEighthesIntegrator
+
+ * @version $Id: GillIntegrator.java 1705 2006-09-17 19:57:39Z luc $
+ * @author L. Maisonobe
+
+ */
+
+public class GillIntegrator
+  extends RungeKuttaIntegrator {
+
+  private static final String methodName = "Gill";
+
+  private static final double sqrt2 = Math.sqrt(2.0);
+
+  private static final double[] c = {
+    1.0 / 2.0, 1.0 / 2.0, 1.0
+  };
+
+  private static final double[][] a = {
+    { 1.0 / 2.0 },
+    { (sqrt2 - 1.0) / 2.0, (2.0 - sqrt2) / 2.0 },
+    { 0.0, -sqrt2 / 2.0, (2.0 + sqrt2) / 2.0 }
+  };
+
+  private static final double[] b = {
+    1.0 / 6.0, (2.0 - sqrt2) / 6.0, (2.0 + sqrt2) / 6.0, 1.0 / 6.0
+  };
+
+  /** Simple constructor.
+   * Build a fourth-order Gill integrator with the given step.
+   * @param step integration step
+   */
+  public GillIntegrator(double step) {
+    super(false, c, a, b, new GillStepInterpolator(), 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/GillIntegrator.java
------------------------------------------------------------------------------
    svn:eol-style = native

Added: jakarta/commons/proper/math/trunk/src/java/org/apache/commons/math/ode/GillStepInterpolator.java
URL: http://svn.apache.org/viewvc/jakarta/commons/proper/math/trunk/src/java/org/apache/commons/math/ode/GillStepInterpolator.java?view=auto&rev=512066
==============================================================================
--- jakarta/commons/proper/math/trunk/src/java/org/apache/commons/math/ode/GillStepInterpolator.java (added)
+++ jakarta/commons/proper/math/trunk/src/java/org/apache/commons/math/ode/GillStepInterpolator.java Mon Feb 26 15:12:40 2007
@@ -0,0 +1,119 @@
+// 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 the Gill fourth
+ * order Runge-Kutta integrator.
+
+ * <p>This interpolator allows 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/6) [ (1 - theta) (1 - 4 theta) y'_1
+ *                                        + (1 - theta) (1 + 2 theta) ((2-q) y'_2 + (2+q) y'_3)
+ *                                        + (1 + theta + 4 theta^2) y'_4
+ *                                        ]
+ * </pre>
+ * where theta belongs to [0 ; 1], q = sqrt(2) and where y'_1 to y'_4
+ * are the four evaluations of the derivatives already computed during
+ * the step.</p>
+
+ * @see GillIntegrator
+
+ * @version $Id: GillStepInterpolator.java 1705 2006-09-17 19:57:39Z luc $
+ * @author L. Maisonobe
+
+ */
+
+class GillStepInterpolator
+  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 GillStepInterpolator() {
+  }
+
+  /** 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 GillStepInterpolator(GillStepInterpolator 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 GillStepInterpolator(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 fourTheta = 4 * theta;
+    double s         = oneMinusThetaH / 6.0;
+    double soMt      = s * (1 - theta);
+    double c23       = soMt * (1 + 2 * theta);
+    double coeff1    = soMt * (1 - fourTheta);
+    double coeff2    = c23  * tMq;
+    double coeff3    = c23  * tPq;
+    double coeff4    = s * (1 + theta * (1 + fourTheta));
+
+    for (int i = 0; i < interpolatedState.length; ++i) {
+      interpolatedState[i] = currentState[i]
+                            - coeff1 * yDotK[0][i] - coeff2 * yDotK[1][i]
+                            - coeff3 * yDotK[2][i] - coeff4 * yDotK[3][i];
+     }
+
+  }
+
+  /** First Gill coefficient. */
+  private static final double tMq = 2 - Math.sqrt(2.0);
+
+  /** Second Gill coefficient. */
+  private static final double tPq = 2 + Math.sqrt(2.0);
+
+  private static final long serialVersionUID = -107804074496313322L;
+
+}

Propchange: jakarta/commons/proper/math/trunk/src/java/org/apache/commons/math/ode/GillStepInterpolator.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