commons-commits mailing list archives

Site index · List index
Message view « Date » · « Thread »
Top « Date » · « Thread »
From l..@apache.org
Subject svn commit: r580753 - in /commons/proper/math/trunk/src/java/org/apache/commons/math/ode: AdaptiveStepsizeIntegrator.java FirstOrderIntegrator.java GraggBulirschStoerIntegrator.java RungeKuttaFehlbergIntegrator.java RungeKuttaIntegrator.java
Date Sun, 30 Sep 2007 16:58:51 GMT
Author: luc
Date: Sun Sep 30 09:58:50 2007
New Revision: 580753

URL: http://svn.apache.org/viewvc?rev=580753&view=rev
Log:
added the two methods getCurrentStepStart and getCurrentStepsize to interface FirstOrderIntegrator
allowing ODE problems to retrieve the current step start and size during integration for each
step trial (i.e. even before the step is accepted)

Modified:
    commons/proper/math/trunk/src/java/org/apache/commons/math/ode/AdaptiveStepsizeIntegrator.java
    commons/proper/math/trunk/src/java/org/apache/commons/math/ode/FirstOrderIntegrator.java
    commons/proper/math/trunk/src/java/org/apache/commons/math/ode/GraggBulirschStoerIntegrator.java
    commons/proper/math/trunk/src/java/org/apache/commons/math/ode/RungeKuttaFehlbergIntegrator.java
    commons/proper/math/trunk/src/java/org/apache/commons/math/ode/RungeKuttaIntegrator.java

Modified: commons/proper/math/trunk/src/java/org/apache/commons/math/ode/AdaptiveStepsizeIntegrator.java
URL: http://svn.apache.org/viewvc/commons/proper/math/trunk/src/java/org/apache/commons/math/ode/AdaptiveStepsizeIntegrator.java?rev=580753&r1=580752&r2=580753&view=diff
==============================================================================
--- commons/proper/math/trunk/src/java/org/apache/commons/math/ode/AdaptiveStepsizeIntegrator.java
(original)
+++ commons/proper/math/trunk/src/java/org/apache/commons/math/ode/AdaptiveStepsizeIntegrator.java
Sun Sep 30 09:58:50 2007
@@ -75,6 +75,8 @@
 
     switchesHandler = new SwitchingFunctionsHandler();
 
+    resetInternalState();
+
   }
 
   /** Build an integrator with the given stepsize bounds.
@@ -104,6 +106,8 @@
 
     switchesHandler = new SwitchingFunctionsHandler();
 
+    resetInternalState();
+
   }
 
   /** Set the initial step size.
@@ -324,6 +328,20 @@
                                   double t, double[] y)
     throws DerivativeException, IntegratorException;
 
+  public double getCurrentStepStart() {
+    return stepStart;
+  }
+
+  public double getCurrentStepsize() {
+    return stepSize;
+  }
+
+  /** Reset internal state to dummy values. */
+  protected void resetInternalState() {
+    stepStart = Double.NaN;
+    stepSize  = Math.sqrt(minStep * maxStep);
+  }
+
   /** Get the minimal step.
    * @return minimal step
    */
@@ -364,5 +382,11 @@
 
   /** Switching functions handler. */
   protected SwitchingFunctionsHandler switchesHandler;
+
+  /** Current step start time. */
+  protected double stepStart;
+
+  /** Current stepsize. */
+  protected double stepSize;
 
 }

Modified: commons/proper/math/trunk/src/java/org/apache/commons/math/ode/FirstOrderIntegrator.java
URL: http://svn.apache.org/viewvc/commons/proper/math/trunk/src/java/org/apache/commons/math/ode/FirstOrderIntegrator.java?rev=580753&r1=580752&r2=580753&view=diff
==============================================================================
--- commons/proper/math/trunk/src/java/org/apache/commons/math/ode/FirstOrderIntegrator.java
(original)
+++ commons/proper/math/trunk/src/java/org/apache/commons/math/ode/FirstOrderIntegrator.java
Sun Sep 30 09:58:50 2007
@@ -63,12 +63,16 @@
                                    double maxCheckInterval,
                                    double convergence);
 
-  /** Integrate the differential equations up to the given time
+  /** Integrate the differential equations up to the given time.
+   * <p>This method solves an Initial Value Problem (IVP).</p>
+   * <p>Since this method stores some internal state variables made
+   * available in its public interface during integration ({@link
+   * #getCurrentStepsize()}), it is <em>not</em> thread-safe.</p>
    * @param equations differential equations to integrate
    * @param t0 initial time
    * @param y0 initial value of the state vector at t0
    * @param t target time for the integration
-   * (can be set to a value smaller thant <code>t0</code> for backward integration)
+   * (can be set to a value smaller than <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
@@ -79,5 +83,27 @@
                          double t0, double[] y0,
                          double t, double[] y)
     throws DerivativeException, IntegratorException;
+
+  /** Get the current value of the step start time t<sub>i</sub>.
+   * <p>This method can be called during integration (typically by
+   * the object implementing the {@link FirstOrderDifferentialEquations
+   * differential equations} problem) if the value of the current step that
+   * is attempted is needed.</p>
+   * <p>The result is undefined if the method is called outside of
+   * calls to {@link #integrate}</p>
+   * @return current value of the step start time t<sub>i</sub>
+   */
+  public double getCurrentStepStart();
+
+  /** Get the current value of the integration stepsize.
+   * <p>This method can be called during integration (typically by
+   * the object implementing the {@link FirstOrderDifferentialEquations
+   * differential equations} problem) if the value of the current stepsize
+   * that is tried is needed.</p>
+   * <p>The result is undefined if the method is called outside of
+   * calls to {@link #integrate}</p>
+   * @return current value of the stepsize
+   */
+  public double getCurrentStepsize();
 
 }

Modified: commons/proper/math/trunk/src/java/org/apache/commons/math/ode/GraggBulirschStoerIntegrator.java
URL: http://svn.apache.org/viewvc/commons/proper/math/trunk/src/java/org/apache/commons/math/ode/GraggBulirschStoerIntegrator.java?rev=580753&r1=580752&r2=580753&view=diff
==============================================================================
--- commons/proper/math/trunk/src/java/org/apache/commons/math/ode/GraggBulirschStoerIntegrator.java
(original)
+++ commons/proper/math/trunk/src/java/org/apache/commons/math/ode/GraggBulirschStoerIntegrator.java
Sun Sep 30 09:58:50 2007
@@ -579,7 +579,7 @@
     }
     interpolator.storeTime(t0);
 
-    double  currentT         = t0;
+    stepStart = t0;
     double  hNew             = 0;
     double  maxError         = Double.MAX_VALUE;
     boolean previousRejected = false;
@@ -591,7 +591,6 @@
     costPerTimeUnit[0] = 0;
     while (! lastStep) {
 
-      double h;
       double error;
       boolean reject = false;
 
@@ -601,14 +600,14 @@
 
         // first evaluation, at the beginning of the step
         if (! firstStepAlreadyComputed) {
-          equations.computeDerivatives(currentT, y, yDot0);
+          equations.computeDerivatives(stepStart, y, yDot0);
         }
 
         if (firstTime) {
 
           hNew = initializeStep(equations, forward,
                                 2 * targetIter + 1, scale,
-                                currentT, y, yDot0, yTmp, yTmpDot);
+                                stepStart, y, yDot0, yTmp, yTmpDot);
 
           if (! forward) {
             hNew = -hNew;
@@ -620,14 +619,14 @@
 
       }
 
-      h = hNew;
+      stepSize = hNew;
 
       // step adjustment near bounds
-      if ((forward && (currentT + h > t))
-          || ((! forward) && (currentT + h < t))) {
-        h = t - currentT;
+      if ((forward && (stepStart + stepSize > t))
+          || ((! forward) && (stepStart + stepSize < t))) {
+        stepSize = t - stepStart;
       }
-      double nextT = currentT + h;
+      double nextT = stepStart + stepSize;
       lastStep = forward ? (nextT >= t) : (nextT <= t);
 
       // iterate over several substep sizes
@@ -637,13 +636,13 @@
         ++k;
 
         // modified midpoint integration with the current substep
-        if ( ! tryStep(equations, currentT, y, h, k, scale, fk[k],
+        if ( ! tryStep(equations, stepStart, y, stepSize, 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));
+          hNew   = Math.abs(filterStep(stepSize * stabilityReduction, false));
           reject = true;
           loop   = false;
 
@@ -667,7 +666,7 @@
 
             if ((error > 1.0e15) || ((k > 1) && (error > maxError))) {
               // error is too big, we reduce the global step
-              hNew   = Math.abs(filterStep(h * stabilityReduction, false));
+              hNew   = Math.abs(filterStep(stepSize * stabilityReduction, false));
               reject = true;
               loop   = false;
             } else {
@@ -679,7 +678,7 @@
               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));
+              optimalStep[k]     = Math.abs(filterStep(stepSize * fac, true));
               costPerTimeUnit[k] = costPerStep[k] / optimalStep[k];
 
               // check convergence
@@ -775,7 +774,7 @@
         }
 
         // derivative at end of step
-        equations.computeDerivatives(currentT + h, y1, yDot1);
+        equations.computeDerivatives(stepStart + stepSize, y1, yDot1);
 
         int mu = 2 * k - mudif + 3;
 
@@ -797,7 +796,7 @@
             extrapolate(l2, j, diagonal, yMidDots[l+1]);
           }
           for (int i = 0; i < y0.length; ++i) {
-            yMidDots[l+1][i] *= h;
+            yMidDots[l+1][i] *= stepSize;
           }
 
           // compute centered differences to evaluate next derivatives
@@ -816,13 +815,13 @@
           // estimate the dense output coefficients
           GraggBulirschStoerStepInterpolator gbsInterpolator
             = (GraggBulirschStoerStepInterpolator) interpolator;
-          gbsInterpolator.computeCoefficients(mu, h);
+          gbsInterpolator.computeCoefficients(mu, stepSize);
 
           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));
+            hInt = Math.abs(stepSize / Math.max(Math.pow(interpError, 1.0 / (mu+4)),
+                                                0.01));
             if (interpError > 10.0) {
               hNew = hInt;
               reject = true;
@@ -831,10 +830,10 @@
 
           // Switching functions handling
           if (!reject) {
-            interpolator.storeTime(currentT + h);
+            interpolator.storeTime(stepStart + stepSize);
             if (switchesHandler.evaluateStep(interpolator)) {
               reject = true;
-              hNew = Math.abs(switchesHandler.getEventTime() - currentT);
+              hNew = Math.abs(switchesHandler.getEventTime() - stepStart);
             }
           }
 
@@ -851,19 +850,19 @@
       if (! reject) {
 
         // store end of step state
-        currentT += h;
+        stepStart += stepSize;
         System.arraycopy(y1, 0, y, 0, y0.length);
 
-        switchesHandler.stepAccepted(currentT, y);
+        switchesHandler.stepAccepted(stepStart, y);
         if (switchesHandler.stop()) {
           lastStep = true;
         }
 
         // provide the step data to the step handler
-        interpolator.storeTime(currentT);
+        interpolator.storeTime(stepStart);
         handler.handleStep(interpolator, lastStep);
 
-        if (switchesHandler.reset(currentT, y) && ! lastStep) {
+        if (switchesHandler.reset(stepStart, y) && ! lastStep) {
           // some switching function has triggered changes that
           // invalidate the derivatives, we need to recompute them
           firstStepAlreadyComputed = false;
@@ -897,7 +896,7 @@
           // after a rejected step neither order nor stepsize
           // should increase
           targetIter = Math.min(optimalIter, k);
-          hNew = Math.min(Math.abs(h), optimalStep[targetIter]);
+          hNew = Math.min(Math.abs(stepSize), optimalStep[targetIter]);
         } else {
           // stepsize control
           if (optimalIter <= k) {

Modified: commons/proper/math/trunk/src/java/org/apache/commons/math/ode/RungeKuttaFehlbergIntegrator.java
URL: http://svn.apache.org/viewvc/commons/proper/math/trunk/src/java/org/apache/commons/math/ode/RungeKuttaFehlbergIntegrator.java?rev=580753&r1=580752&r2=580753&view=diff
==============================================================================
--- commons/proper/math/trunk/src/java/org/apache/commons/math/ode/RungeKuttaFehlbergIntegrator.java
(original)
+++ commons/proper/math/trunk/src/java/org/apache/commons/math/ode/RungeKuttaFehlbergIntegrator.java
Sun Sep 30 09:58:50 2007
@@ -187,7 +187,7 @@
     }
     interpolator.storeTime(t0);
 
-    double  currentT  = t0;
+    stepStart  = t0;
     double  hNew      = 0;
     boolean firstTime = true;
     boolean lastStep;
@@ -196,13 +196,12 @@
 
       interpolator.shift();
 
-      double h     = 0;
       double error = 0;
       for (boolean loop = true; loop;) {
 
         if (firstTime || !fsal) {
           // first stage
-          equations.computeDerivatives(currentT, y, yDotK[0]);
+          equations.computeDerivatives(stepStart, y, yDotK[0]);
         }
 
         if (firstTime) {
@@ -216,16 +215,16 @@
             }
           }
           hNew = initializeStep(equations, forward, getOrder(), scale,
-                                currentT, y, yDotK[0], yTmp, yDotK[1]);
+                                stepStart, y, yDotK[0], yTmp, yDotK[1]);
           firstTime = false;
         }
 
-        h = hNew;
+        stepSize = hNew;
 
         // step adjustment near bounds
-        if ((forward && (currentT + h > t))
-            || ((! forward) && (currentT + h < t))) {
-          h = t - currentT;
+        if ((forward && (stepStart + stepSize > t))
+            || ((! forward) && (stepStart + stepSize < t))) {
+          stepSize = t - stepStart;
         }
 
         // next stages
@@ -236,10 +235,10 @@
             for (int l = 1; l < k; ++l) {
               sum += a[k-1][l] * yDotK[l][j];
             }
-            yTmp[j] = y[j] + h * sum;
+            yTmp[j] = y[j] + stepSize * sum;
           }
 
-          equations.computeDerivatives(currentT + c[k-1] * h, yTmp, yDotK[k]);
+          equations.computeDerivatives(stepStart + c[k-1] * stepSize, yTmp, yDotK[k]);
 
         }
 
@@ -249,18 +248,18 @@
           for (int l = 1; l < stages; ++l) {
             sum    += b[l] * yDotK[l][j];
           }
-          yTmp[j] = y[j] + h * sum;
+          yTmp[j] = y[j] + stepSize * sum;
         }
 
         // estimate the error at the end of the step
-        error = estimateError(yDotK, y, yTmp, h);
+        error = estimateError(yDotK, y, yTmp, stepSize);
         if (error <= 1.0) {
 
           // Switching functions handling
-          interpolator.storeTime(currentT + h);
+          interpolator.storeTime(stepStart + stepSize);
           if (switchesHandler.evaluateStep(interpolator)) {
             // reject the step to match exactly the next switch time
-            hNew = switchesHandler.getEventTime() - currentT;
+            hNew = switchesHandler.getEventTime() - stepStart;
           } else {
             // accept the step
             loop = false;
@@ -271,23 +270,23 @@
           double factor = Math.min(maxGrowth,
                                    Math.max(minReduction,
                                             safety * Math.pow(error, exp)));
-          hNew = filterStep(h * factor, false);
+          hNew = filterStep(stepSize * factor, false);
         }
 
       }
 
       // the step has been accepted
-      currentT += h;
+      stepStart += stepSize;
       System.arraycopy(yTmp, 0, y, 0, y0.length);
-      switchesHandler.stepAccepted(currentT, y);
+      switchesHandler.stepAccepted(stepStart, y);
       if (switchesHandler.stop()) {
         lastStep = true;
       } else {
-        lastStep = forward ? (currentT >= t) : (currentT <= t);
+        lastStep = forward ? (stepStart >= t) : (stepStart <= t);
       }
 
       // provide the step data to the step handler
-      interpolator.storeTime(currentT);
+      interpolator.storeTime(stepStart);
       handler.handleStep(interpolator, lastStep);
 
       if (fsal) {
@@ -295,10 +294,10 @@
         System.arraycopy(yDotK[stages - 1], 0, yDotK[0], 0, y0.length);
       }
 
-      if (switchesHandler.reset(currentT, y) && ! lastStep) {
+      if (switchesHandler.reset(stepStart, y) && ! lastStep) {
         // some switching function has triggered changes that
         // invalidate the derivatives, we need to recompute them
-        equations.computeDerivatives(currentT, y, yDotK[0]);
+        equations.computeDerivatives(stepStart, y, yDotK[0]);
       }
 
       if (! lastStep) {
@@ -306,13 +305,15 @@
         double  factor     = Math.min(maxGrowth,
                                       Math.max(minReduction,
                                                safety * Math.pow(error, exp)));
-        double  scaledH    = h * factor;
-        double  nextT      = currentT + scaledH;
+        double  scaledH    = stepSize * factor;
+        double  nextT      = stepStart + scaledH;
         boolean nextIsLast = forward ? (nextT >= t) : (nextT <= t);
         hNew = filterStep(scaledH, nextIsLast);
       }
 
     } while (! lastStep);
+
+    resetInternalState();
 
   }
 

Modified: commons/proper/math/trunk/src/java/org/apache/commons/math/ode/RungeKuttaIntegrator.java
URL: http://svn.apache.org/viewvc/commons/proper/math/trunk/src/java/org/apache/commons/math/ode/RungeKuttaIntegrator.java?rev=580753&r1=580752&r2=580753&view=diff
==============================================================================
--- commons/proper/math/trunk/src/java/org/apache/commons/math/ode/RungeKuttaIntegrator.java
(original)
+++ commons/proper/math/trunk/src/java/org/apache/commons/math/ode/RungeKuttaIntegrator.java
Sun Sep 30 09:58:50 2007
@@ -78,6 +78,7 @@
     this.step       = step;
     handler         = DummyStepHandler.getInstance();
     switchesHandler = new SwitchingFunctionsHandler();
+    resetInternalState();
   }
 
   /** Get the name of the method.
@@ -180,11 +181,11 @@
     interpolator.storeTime(t0);
 
     // recompute the step
-    double  currentT  = t0;
     long    nbStep    = Math.max(1l, Math.abs(Math.round((t - t0) / step)));
-    double  h         = (t - t0) / nbStep;
     boolean firstTime = true;
     boolean lastStep  = false;
+    stepStart = t0;
+    stepSize  = (t - t0) / nbStep;
     handler.reset();
     for (long i = 0; ! lastStep; ++i) {
 
@@ -195,7 +196,7 @@
 
         if (firstTime || !fsal) {
           // first stage
-          equations.computeDerivatives(currentT, y, yDotK[0]);
+          equations.computeDerivatives(stepStart, y, yDotK[0]);
           firstTime = false;
         }
 
@@ -207,10 +208,10 @@
             for (int l = 1; l < k; ++l) {
               sum += a[k-1][l] * yDotK[l][j];
             }
-            yTmp[j] = y[j] + h * sum;
+            yTmp[j] = y[j] + stepSize * sum;
           }
 
-          equations.computeDerivatives(currentT + c[k-1] * h, yTmp, yDotK[k]);
+          equations.computeDerivatives(stepStart + c[k-1] * stepSize, yTmp, yDotK[k]);
 
         }
 
@@ -220,14 +221,14 @@
           for (int l = 1; l < stages; ++l) {
             sum    += b[l] * yDotK[l][j];
           }
-          yTmp[j] = y[j] + h * sum;
+          yTmp[j] = y[j] + stepSize * sum;
         }
 
         // Switching functions handling
-        interpolator.storeTime(currentT + h);
+        interpolator.storeTime(stepStart + stepSize);
         if (switchesHandler.evaluateStep(interpolator)) {
           needUpdate = true;
-          h = switchesHandler.getEventTime() - currentT;
+          stepSize = switchesHandler.getEventTime() - stepStart;
         } else {
           loop = false;
         }
@@ -235,9 +236,9 @@
       }
 
       // the step has been accepted
-      currentT += h;
+      stepStart += stepSize;
       System.arraycopy(yTmp, 0, y, 0, y0.length);
-      switchesHandler.stepAccepted(currentT, y);
+      switchesHandler.stepAccepted(stepStart, y);
       if (switchesHandler.stop()) {
         lastStep = true;
       } else {
@@ -245,7 +246,7 @@
       }
 
       // provide the step data to the step handler
-      interpolator.storeTime(currentT);
+      interpolator.storeTime(stepStart);
       handler.handleStep(interpolator, lastStep);
 
       if (fsal) {
@@ -253,22 +254,38 @@
         System.arraycopy(yDotK[stages - 1], 0, yDotK[0], 0, y0.length);
       }
 
-      if (switchesHandler.reset(currentT, y) && ! lastStep) {
+      if (switchesHandler.reset(stepStart, y) && ! lastStep) {
         // some switching function has triggered changes that
         // invalidate the derivatives, we need to recompute them
-        equations.computeDerivatives(currentT, y, yDotK[0]);
+        equations.computeDerivatives(stepStart, y, yDotK[0]);
       }
 
       if (needUpdate) {
         // a switching function has changed the step
         // we need to recompute stepsize
-        nbStep = Math.max(1l, Math.abs(Math.round((t - currentT) / step)));
-        h = (t - currentT) / nbStep;
+        nbStep = Math.max(1l, Math.abs(Math.round((t - stepStart) / step)));
+        stepSize = (t - stepStart) / nbStep;
         i = -1;
       }
 
     }
 
+    resetInternalState();
+
+  }
+
+  public double getCurrentStepStart() {
+    return stepStart;
+  }
+
+  public double getCurrentStepsize() {
+    return stepSize;
+  }
+
+  /** Reset internal state to dummy values. */
+  private void resetInternalState() {
+    stepStart = Double.NaN;
+    stepSize  = Double.NaN;
   }
 
   /** Indicator for <i>fsal</i> methods. */
@@ -294,5 +311,11 @@
 
   /** Switching functions handler. */
   protected SwitchingFunctionsHandler switchesHandler;
+
+  /** Current step start time. */
+  private double stepStart;
+
+  /** Current stepsize. */
+  private double stepSize;
 
 }



Mime
View raw message