commons-commits mailing list archives

Site index · List index
Message view « Date » · « Thread »
Top « Date » · « Thread »
From l..@apache.org
Subject svn commit: r1588753 - in /commons/proper/math/trunk/src: changes/ main/java/org/apache/commons/math3/ode/ main/java/org/apache/commons/math3/ode/nonstiff/ site/xdoc/userguide/ test/java/org/apache/commons/math3/ode/nonstiff/
Date Sun, 20 Apr 2014 13:25:12 GMT
Author: luc
Date: Sun Apr 20 13:25:11 2014
New Revision: 1588753

URL: http://svn.apache.org/r1588753
Log:
Added an order 6 fixed-step ODE integrator.

The integrator was designed by H. A. Luther in 1968. We have added a
corresponding step interpolator by solving the order conditions provided
by the rkcheck tool.

Added:
    commons/proper/math/trunk/src/main/java/org/apache/commons/math3/ode/nonstiff/LutherIntegrator.java
  (with props)
    commons/proper/math/trunk/src/main/java/org/apache/commons/math3/ode/nonstiff/LutherStepInterpolator.java
  (with props)
    commons/proper/math/trunk/src/test/java/org/apache/commons/math3/ode/nonstiff/LutherIntegratorTest.java
  (with props)
    commons/proper/math/trunk/src/test/java/org/apache/commons/math3/ode/nonstiff/LutherStepInterpolatorTest.java
  (with props)
Modified:
    commons/proper/math/trunk/src/changes/changes.xml
    commons/proper/math/trunk/src/main/java/org/apache/commons/math3/ode/package-info.java
    commons/proper/math/trunk/src/site/xdoc/userguide/ode.xml

Modified: commons/proper/math/trunk/src/changes/changes.xml
URL: http://svn.apache.org/viewvc/commons/proper/math/trunk/src/changes/changes.xml?rev=1588753&r1=1588752&r2=1588753&view=diff
==============================================================================
--- commons/proper/math/trunk/src/changes/changes.xml (original)
+++ commons/proper/math/trunk/src/changes/changes.xml Sun Apr 20 13:25:11 2014
@@ -56,6 +56,9 @@ If the output is not quite correct, chec
         Added new methods for testing floating-point equality between the real
         (resp. imaginary) parts of two complex numbers.
       </action>
+      <action dev="luc" type="add" >
+        Added an order 6 fixed-step ODE integrator designed by H. A. Luther in 1968.
+      </action>
       <action dev="luc" type="update" >
         Bracketing utility for univariate root solvers returns a tighter interval than before.
         It also allows choosing the search interval expansion rate, supporting both linear

Added: commons/proper/math/trunk/src/main/java/org/apache/commons/math3/ode/nonstiff/LutherIntegrator.java
URL: http://svn.apache.org/viewvc/commons/proper/math/trunk/src/main/java/org/apache/commons/math3/ode/nonstiff/LutherIntegrator.java?rev=1588753&view=auto
==============================================================================
--- commons/proper/math/trunk/src/main/java/org/apache/commons/math3/ode/nonstiff/LutherIntegrator.java
(added)
+++ commons/proper/math/trunk/src/main/java/org/apache/commons/math3/ode/nonstiff/LutherIntegrator.java
Sun Apr 20 13:25:11 2014
@@ -0,0 +1,90 @@
+/*
+ * Licensed to the Apache Software Foundation (ASF) under one or more
+ * contributor license agreements.  See the NOTICE file distributed with
+ * this work for additional information regarding copyright ownership.
+ * The ASF licenses this file to You under the Apache License, Version 2.0
+ * (the "License"); you may not use this file except in compliance with
+ * the License.  You may obtain a copy of the License at
+ *
+ *      http://www.apache.org/licenses/LICENSE-2.0
+ *
+ * Unless required by applicable law or agreed to in writing, software
+ * distributed under the License is distributed on an "AS IS" BASIS,
+ * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+ * See the License for the specific language governing permissions and
+ * limitations under the License.
+ */
+
+package org.apache.commons.math3.ode.nonstiff;
+
+import org.apache.commons.math3.util.FastMath;
+
+
+/**
+ * This class implements the Luther sixth order Runge-Kutta
+ * integrator for Ordinary Differential Equations.
+
+ * <p>
+ * This method is described in H. A. Luther 1968 paper <a
+ * href="http://www.ams.org/journals/mcom/1968-22-102/S0025-5718-68-99876-1/S0025-5718-68-99876-1.pdf">
+ * An explicit Sixth-Order Runge-Kutta Formula</a>.
+ * </p>
+
+ * <p>This method is an explicit Runge-Kutta method, its Butcher-array
+ * is the following one :
+ * <pre>
+ *        0   |               0                     0                     0             
       0                     0                     0
+ *        1   |               1                     0                     0             
       0                     0                     0
+ *       1/2  |              3/8                   1/8                    0             
       0                     0                     0
+ *       2/3  |              8/27                  2/27                  8/27           
       0                     0                     0
+ *   (7-q)/14 | (  -21 +   9q)/392    (  -56 +   8q)/392    (  336 -  48q)/392    (  -63
+   3q)/392                  0                     0
+ *   (7+q)/14 | (-1155 - 255q)/1960   ( -280 -  40q)/1960   (    0 - 320q)/1960   (   63
+ 363q)/1960   ( 2352 + 392q)/1960                 0
+ *        1   | (  330 + 105q)/180    (  120 +   0q)/180    ( -200 + 280q)/180    (  126
- 189q)/180    ( -686 - 126q)/180     ( 490 -  70q)/180
+ *            |--------------------------------------------------------------------------------------------------------------------------------------------------
+ *            |              1/20                   0                   16/45           
      0                   49/180                 49/180         1/20
+ * </pre>
+ * where q = &radic;21</p>
+ *
+ * @see EulerIntegrator
+ * @see ClassicalRungeKuttaIntegrator
+ * @see GillIntegrator
+ * @see MidpointIntegrator
+ * @see ThreeEighthesIntegrator
+ * @version $Id$
+ * @since 3.3
+ */
+
+public class LutherIntegrator extends RungeKuttaIntegrator {
+
+    /** Square root. */
+    private static final double Q = FastMath.sqrt(21);
+
+    /** Time steps Butcher array. */
+    private static final double[] STATIC_C = {
+        1.0, 1.0 / 2.0, 2.0 / 3.0, (7.0 - Q) / 14.0, (7.0 + Q) / 14.0, 1.0
+    };
+
+    /** Internal weights Butcher array. */
+    private static final double[][] STATIC_A = {
+        {                      1.0        },
+        {                   3.0 /   8.0,                  1.0 /   8.0  },
+        {                   8.0 /   27.0,                 2.0 /   27.0,                 
8.0 /   27.0  },
+        { (  -21.0 +   9.0 * Q) /  392.0, ( -56.0 +  8.0 * Q) /  392.0, ( 336.0 -  48.0 *
Q) /  392.0, (-63.0 +   3.0 * Q) /  392.0 },
+        { (-1155.0 - 255.0 * Q) / 1960.0, (-280.0 - 40.0 * Q) / 1960.0, (   0.0 - 320.0 *
Q) / 1960.0, ( 63.0 + 363.0 * Q) / 1960.0,   (2352.0 + 392.0 * Q) / 1960.0 },
+        { (  330.0 + 105.0 * Q) /  180.0, ( 120.0 +  0.0 * Q) /  180.0, (-200.0 + 280.0 *
Q) /  180.0, (126.0 - 189.0 * Q) /  180.0,   (-686.0 - 126.0 * Q) /  180.0,   (490.0 -  70.0
* Q) / 180.0 }
+    };
+
+    /** Propagation weights Butcher array. */
+    private static final double[] STATIC_B = {
+        1.0 / 20.0, 0, 16.0 / 45.0, 0, 49.0 / 180.0, 49.0 / 180.0, 1.0 / 20.0
+    };
+
+    /** Simple constructor.
+     * Build a fourth-order Luther integrator with the given step.
+     * @param step integration step
+     */
+    public LutherIntegrator(final double step) {
+        super("Luther", STATIC_C, STATIC_A, STATIC_B, new LutherStepInterpolator(), step);
+    }
+
+}

Propchange: commons/proper/math/trunk/src/main/java/org/apache/commons/math3/ode/nonstiff/LutherIntegrator.java
------------------------------------------------------------------------------
    svn:eol-style = native

Propchange: commons/proper/math/trunk/src/main/java/org/apache/commons/math3/ode/nonstiff/LutherIntegrator.java
------------------------------------------------------------------------------
    svn:keywords = "Author Date Id Revision"

Added: commons/proper/math/trunk/src/main/java/org/apache/commons/math3/ode/nonstiff/LutherStepInterpolator.java
URL: http://svn.apache.org/viewvc/commons/proper/math/trunk/src/main/java/org/apache/commons/math3/ode/nonstiff/LutherStepInterpolator.java?rev=1588753&view=auto
==============================================================================
--- commons/proper/math/trunk/src/main/java/org/apache/commons/math3/ode/nonstiff/LutherStepInterpolator.java
(added)
+++ commons/proper/math/trunk/src/main/java/org/apache/commons/math3/ode/nonstiff/LutherStepInterpolator.java
Sun Apr 20 13:25:11 2014
@@ -0,0 +1,180 @@
+/*
+ * Licensed to the Apache Software Foundation (ASF) under one or more
+ * contributor license agreements.  See the NOTICE file distributed with
+ * this work for additional information regarding copyright ownership.
+ * The ASF licenses this file to You under the Apache License, Version 2.0
+ * (the "License"); you may not use this file except in compliance with
+ * the License.  You may obtain a copy of the License at
+ *
+ *      http://www.apache.org/licenses/LICENSE-2.0
+ *
+ * Unless required by applicable law or agreed to in writing, software
+ * distributed under the License is distributed on an "AS IS" BASIS,
+ * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+ * See the License for the specific language governing permissions and
+ * limitations under the License.
+ */
+
+package org.apache.commons.math3.ode.nonstiff;
+
+import org.apache.commons.math3.ode.sampling.StepInterpolator;
+import org.apache.commons.math3.util.FastMath;
+
+/**
+ * This class implements a step interpolator for second order
+ * Runge-Kutta integrator.
+ *
+ * <p>This interpolator computes dense output inside the last
+ * step computed. The interpolation equation is consistent with the
+ * integration scheme.</p>
+ *
+ * @see LutherIntegrator
+ * @version $Id$
+ * @since 3.3
+ */
+
+class LutherStepInterpolator extends RungeKuttaStepInterpolator {
+
+    /** Serializable version identifier */
+    private static final long serialVersionUID = 20140416L;
+
+    /** Square root. */
+    private static final double Q = FastMath.sqrt(21);
+
+    /** Simple constructor.
+     * This constructor builds an instance that is not usable yet, the
+     * {@link
+     * org.apache.commons.math3.ode.sampling.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 later initializing the copy.
+     */
+    public LutherStepInterpolator() {
+    }
+
+    /** 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 LutherStepInterpolator(final LutherStepInterpolator interpolator) {
+        super(interpolator);
+    }
+
+    /** {@inheritDoc} */
+    @Override
+    protected StepInterpolator doCopy() {
+        return new LutherStepInterpolator(this);
+    }
+
+
+    /** {@inheritDoc} */
+    @Override
+    protected void computeInterpolatedStateAndDerivatives(final double theta,
+                                                          final double oneMinusThetaH) {
+
+        // the coefficients below have been computed by solving the
+        // order conditions from a theorem from Butcher (1963), using
+        // the method explained in Folkmar Bornemann paper "Runge-Kutta
+        // Methods, Trees, and Maple", Center of Mathematical Sciences, Munich
+        // University of Technology, February 9, 2001
+        //<http://wwwzenger.informatik.tu-muenchen.de/selcuk/sjam012101.html>
+
+        // the method is implemented in the rkcheck tool
+        // <https://www.spaceroots.org/software/rkcheck/index.html>.
+        // Running it for order 5 gives the following order conditions
+        // for an interpolator:
+        // order 1 conditions
+        // \sum_{i=1}^{i=s}\left(b_{i} \right) =1
+        // order 2 conditions
+        // \sum_{i=1}^{i=s}\left(b_{i} c_{i}\right) = \frac{\theta}{2}
+        // order 3 conditions
+        // \sum_{i=2}^{i=s}\left(b_{i} \sum_{j=1}^{j=i-1}{\left(a_{i,j} c_{j} \right)}\right)
= \frac{\theta^{2}}{6}
+        // \sum_{i=1}^{i=s}\left(b_{i} c_{i}^{2}\right) = \frac{\theta^{2}}{3}
+        // order 4 conditions
+        // \sum_{i=3}^{i=s}\left(b_{i} \sum_{j=2}^{j=i-1}{\left(a_{i,j} \sum_{k=1}^{k=j-1}{\left(a_{j,k}
c_{k} \right)} \right)}\right) = \frac{\theta^{3}}{24}
+        // \sum_{i=2}^{i=s}\left(b_{i} \sum_{j=1}^{j=i-1}{\left(a_{i,j} c_{j}^{2} \right)}\right)
= \frac{\theta^{3}}{12}
+        // \sum_{i=2}^{i=s}\left(b_{i} c_{i}\sum_{j=1}^{j=i-1}{\left(a_{i,j} c_{j} \right)}\right)
= \frac{\theta^{3}}{8}
+        // \sum_{i=1}^{i=s}\left(b_{i} c_{i}^{3}\right) = \frac{\theta^{3}}{4}
+        // order 5 conditions
+        // \sum_{i=4}^{i=s}\left(b_{i} \sum_{j=3}^{j=i-1}{\left(a_{i,j} \sum_{k=2}^{k=j-1}{\left(a_{j,k}
\sum_{l=1}^{l=k-1}{\left(a_{k,l} c_{l} \right)} \right)} \right)}\right) = \frac{\theta^{4}}{120}
+        // \sum_{i=3}^{i=s}\left(b_{i} \sum_{j=2}^{j=i-1}{\left(a_{i,j} \sum_{k=1}^{k=j-1}{\left(a_{j,k}
c_{k}^{2} \right)} \right)}\right) = \frac{\theta^{4}}{60}
+        // \sum_{i=3}^{i=s}\left(b_{i} \sum_{j=2}^{j=i-1}{\left(a_{i,j} c_{j}\sum_{k=1}^{k=j-1}{\left(a_{j,k}
c_{k} \right)} \right)}\right) = \frac{\theta^{4}}{40}
+        // \sum_{i=2}^{i=s}\left(b_{i} \sum_{j=1}^{j=i-1}{\left(a_{i,j} c_{j}^{3} \right)}\right)
= \frac{\theta^{4}}{20}
+        // \sum_{i=3}^{i=s}\left(b_{i} c_{i}\sum_{j=2}^{j=i-1}{\left(a_{i,j} \sum_{k=1}^{k=j-1}{\left(a_{j,k}
c_{k} \right)} \right)}\right) = \frac{\theta^{4}}{30}
+        // \sum_{i=2}^{i=s}\left(b_{i} c_{i}\sum_{j=1}^{j=i-1}{\left(a_{i,j} c_{j}^{2} \right)}\right)
= \frac{\theta^{4}}{15}
+        // \sum_{i=2}^{i=s}\left(b_{i} \left(\sum_{j=1}^{j=i-1}{\left(a_{i,j} c_{j} \right)}
\right)^{2}\right) = \frac{\theta^{4}}{20}
+        // \sum_{i=2}^{i=s}\left(b_{i} c_{i}^{2}\sum_{j=1}^{j=i-1}{\left(a_{i,j} c_{j} \right)}\right)
= \frac{\theta^{4}}{10}
+        // \sum_{i=1}^{i=s}\left(b_{i} c_{i}^{4}\right) = \frac{\theta^{4}}{5}
+
+        // The a_{j,k} and c_{k} are given by the integrator Butcher arrays. What remains
to solve
+        // are the b_i for the interpolator. They are found by solving the above equations.
+        // For a given interpolator, some equations are redundant, so in our case when we
select
+        // all equations from order 1 to 4, we still don't have enough independent equations
+        // to solve from b_1 to b_7. We need to also select one equation from order 5. Here,
+        // we selected the last equation. It appears this choice implied at least the last
3 equations
+        // are fulfilled, but some of the former ones are not, so the resulting interpolator
is order 5.
+        // At the end, we get the b_i as polynomials in theta.
+
+        final double coeffDot1 =  1 + theta * ( -54            /   5.0 + theta * (   36 
                 + theta * ( -47                   + theta *   21)));
+        final double coeffDot2 =  0;
+        final double coeffDot3 =      theta * (-208            /  15.0 + theta * (  320 
          / 3.0  + theta * (-608            /  3.0 + theta *  112)));
+        final double coeffDot4 =      theta * ( 324            /  25.0 + theta * ( -486 
          / 5.0  + theta * ( 972            /  5.0 + theta * -567           /  5.0)));
+        final double coeffDot5 =      theta * ((833 + 343 * Q) / 150.0 + theta * ((-637 -
357 * Q) / 30.0 + theta * ((392 + 287 * Q) / 15.0 + theta * (-49 - 49 * Q) /  5.0)));
+        final double coeffDot6 =      theta * ((833 - 343 * Q) / 150.0 + theta * ((-637 +
357 * Q) / 30.0 + theta * ((392 - 287 * Q) / 15.0 + theta * (-49 + 49 * Q) /  5.0)));
+        final double coeffDot7 =      theta * (   3            /   5.0 + theta * (   -3 
                 + theta *     3));
+
+        if ((previousState != null) && (theta <= 0.5)) {
+
+            final double coeff1    =  1 + theta * ( -27            /   5.0 + theta * (  
12                   + theta * ( -47            /  4.0 + theta *   21           /  5.0)));
+            final double coeff2    =  0;
+            final double coeff3    =      theta * (-104            /  15.0 + theta * (  320
           / 9.0  + theta * (-152            /  3.0 + theta *  112           /  5.0)));
+            final double coeff4    =      theta * ( 162            /  25.0 + theta * ( -162
           / 5.0  + theta * ( 243            /  5.0 + theta * -567           / 25.0)));
+            final double coeff5    =      theta * ((833 + 343 * Q) / 300.0 + theta * ((-637
- 357 * Q) / 90.0 + theta * ((392 + 287 * Q) / 60.0 + theta * (-49 - 49 * Q) / 25.0)));
+            final double coeff6    =      theta * ((833 - 343 * Q) / 300.0 + theta * ((-637
+ 357 * Q) / 90.0 + theta * ((392 - 287 * Q) / 60.0 + theta * (-49 + 49 * Q) / 25.0)));
+            final double coeff7    =      theta * (   3            /  10.0 + theta * (  
-1                   + theta * (   3            /  4.0)));
+            for (int i = 0; i < interpolatedState.length; ++i) {
+                final double yDot1 = yDotK[0][i];
+                final double yDot2 = yDotK[1][i];
+                final double yDot3 = yDotK[2][i];
+                final double yDot4 = yDotK[3][i];
+                final double yDot5 = yDotK[4][i];
+                final double yDot6 = yDotK[5][i];
+                final double yDot7 = yDotK[6][i];
+                interpolatedState[i] = previousState[i] +
+                        theta * h * (coeff1 * yDot1 + coeff2 * yDot2 + coeff3 * yDot3 +
+                                     coeff4 * yDot4 + coeff5 * yDot5 + coeff6 * yDot6 + coeff7
* yDot7);
+                interpolatedDerivatives[i] = coeffDot1 * yDot1 + coeffDot2 * yDot2 + coeffDot3
* yDot3 +
+                        coeffDot4 * yDot4 + coeffDot5 * yDot5 + coeffDot6 * yDot6 + coeffDot7
* yDot7;
+            }
+        } else {
+
+            final double coeff1    =  -1 /  20.0 + theta * (  19            /  20.0 + theta
* (  -89             /  20.0  + theta * (   151            /  20.0 + theta *  -21        
  /   5.0)));
+            final double coeff2    =  0;
+            final double coeff3    = -16 /  45.0 + theta * ( -16            /  45.0 + theta
* ( -328             /  45.0  + theta * (   424            /  15.0 + theta * -112        
  /   5.0)));
+            final double coeff4    =               theta * (                          theta
* (  162             /  25.0  + theta * (  -648            /  25.0 + theta *  567        
  /  25.0)));
+            final double coeff5    = -49 / 180.0 + theta * ( -49            / 180.0 + theta
* ((2254 + 1029 * Q) / 900.0  + theta * ((-1372 - 847 * Q) / 300.0 + theta * ( 49 + 49 * Q)
/  25.0)));
+            final double coeff6    = -49 / 180.0 + theta * ( -49            / 180.0 + theta
* ((2254 - 1029 * Q) / 900.0  + theta * ((-1372 + 847 * Q) / 300.0 + theta * ( 49 - 49 * Q)
/  25.0)));
+            final double coeff7    =  -1 /  20.0 + theta * (  -1            /  20.0 + theta
* (    1             /   4.0  + theta * (    -3            /   4.0)));
+            for (int i = 0; i < interpolatedState.length; ++i) {
+                final double yDot1 = yDotK[0][i];
+                final double yDot2 = yDotK[1][i];
+                final double yDot3 = yDotK[2][i];
+                final double yDot4 = yDotK[3][i];
+                final double yDot5 = yDotK[4][i];
+                final double yDot6 = yDotK[5][i];
+                final double yDot7 = yDotK[6][i];
+                interpolatedState[i] = currentState[i] +
+                        oneMinusThetaH * (coeff1 * yDot1 + coeff2 * yDot2 + coeff3 * yDot3
+
+                                          coeff4 * yDot4 + coeff5 * yDot5 + coeff6 * yDot6
+ coeff7 * yDot7);
+                interpolatedDerivatives[i] = coeffDot1 * yDot1 + coeffDot2 * yDot2 + coeffDot3
* yDot3 +
+                        coeffDot4 * yDot4 + coeffDot5 * yDot5 + coeffDot6 * yDot6 + coeffDot7
* yDot7;
+            }
+        }
+
+    }
+
+}

Propchange: commons/proper/math/trunk/src/main/java/org/apache/commons/math3/ode/nonstiff/LutherStepInterpolator.java
------------------------------------------------------------------------------
    svn:eol-style = native

Propchange: commons/proper/math/trunk/src/main/java/org/apache/commons/math3/ode/nonstiff/LutherStepInterpolator.java
------------------------------------------------------------------------------
    svn:keywords = "Author Date Id Revision"

Modified: commons/proper/math/trunk/src/main/java/org/apache/commons/math3/ode/package-info.java
URL: http://svn.apache.org/viewvc/commons/proper/math/trunk/src/main/java/org/apache/commons/math3/ode/package-info.java?rev=1588753&r1=1588752&r2=1588753&view=diff
==============================================================================
--- commons/proper/math/trunk/src/main/java/org/apache/commons/math3/ode/package-info.java
(original)
+++ commons/proper/math/trunk/src/main/java/org/apache/commons/math3/ode/package-info.java
Sun Apr 20 13:25:11 2014
@@ -136,6 +136,7 @@
  * <tr><td>{@link org.apache.commons.math3.ode.nonstiff.ClassicalRungeKuttaIntegrator
Classical Runge-Kutta}</td><td>4</td></tr>
  * <tr><td>{@link org.apache.commons.math3.ode.nonstiff.GillIntegrator Gill}</td><td>4</td></tr>
  * <tr><td>{@link org.apache.commons.math3.ode.nonstiff.ThreeEighthesIntegrator
3/8}</td><td>4</td></tr>
+ * <tr><td>{@link org.apache.commons.math3.ode.nonstiff.LutherIntegrator Luther}</td><td>6</td></tr>
  * </table>
  * </p>
  *

Modified: commons/proper/math/trunk/src/site/xdoc/userguide/ode.xml
URL: http://svn.apache.org/viewvc/commons/proper/math/trunk/src/site/xdoc/userguide/ode.xml?rev=1588753&r1=1588752&r2=1588753&view=diff
==============================================================================
--- commons/proper/math/trunk/src/site/xdoc/userguide/ode.xml (original)
+++ commons/proper/math/trunk/src/site/xdoc/userguide/ode.xml Sun Apr 20 13:25:11 2014
@@ -265,6 +265,7 @@ public int eventOccurred(double t, doubl
           <tr><td><a href="../apidocs/org/apache/commons/math3/ode/nonstiff/ClassicalRungeKuttaIntegrator.html">Classical
Runge-Kutta</a></td><td>4</td></tr>
           <tr><td><a href="../apidocs/org/apache/commons/math3/ode/nonstiff/GillIntegrator.html">Gill</a></td><td>4</td></tr>
           <tr><td><a href="../apidocs/org/apache/commons/math3/ode/nonstiff/ThreeEighthesIntegrator.html">3/8</a></td><td>4</td></tr>
+          <tr><td><a href="../apidocs/org/apache/commons/math3/ode/nonstiff/LutherIntegrator.html">Luther</a></td><td>6</td></tr>
           </table>
         </p>
         <p>

Added: commons/proper/math/trunk/src/test/java/org/apache/commons/math3/ode/nonstiff/LutherIntegratorTest.java
URL: http://svn.apache.org/viewvc/commons/proper/math/trunk/src/test/java/org/apache/commons/math3/ode/nonstiff/LutherIntegratorTest.java?rev=1588753&view=auto
==============================================================================
--- commons/proper/math/trunk/src/test/java/org/apache/commons/math3/ode/nonstiff/LutherIntegratorTest.java
(added)
+++ commons/proper/math/trunk/src/test/java/org/apache/commons/math3/ode/nonstiff/LutherIntegratorTest.java
Sun Apr 20 13:25:11 2014
@@ -0,0 +1,309 @@
+/*
+ * Licensed to the Apache Software Foundation (ASF) under one or more
+ * contributor license agreements.  See the NOTICE file distributed with
+ * this work for additional information regarding copyright ownership.
+ * The ASF licenses this file to You under the Apache License, Version 2.0
+ * (the "License"); you may not use this file except in compliance with
+ * the License.  You may obtain a copy of the License at
+ *
+ *      http://www.apache.org/licenses/LICENSE-2.0
+ *
+ * Unless required by applicable law or agreed to in writing, software
+ * distributed under the License is distributed on an "AS IS" BASIS,
+ * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+ * See the License for the specific language governing permissions and
+ * limitations under the License.
+ */
+
+package org.apache.commons.math3.ode.nonstiff;
+
+
+import org.apache.commons.math3.exception.DimensionMismatchException;
+import org.apache.commons.math3.exception.MaxCountExceededException;
+import org.apache.commons.math3.exception.NoBracketingException;
+import org.apache.commons.math3.exception.NumberIsTooSmallException;
+import org.apache.commons.math3.ode.FirstOrderDifferentialEquations;
+import org.apache.commons.math3.ode.FirstOrderIntegrator;
+import org.apache.commons.math3.ode.TestProblem1;
+import org.apache.commons.math3.ode.TestProblem3;
+import org.apache.commons.math3.ode.TestProblem5;
+import org.apache.commons.math3.ode.TestProblemAbstract;
+import org.apache.commons.math3.ode.TestProblemFactory;
+import org.apache.commons.math3.ode.TestProblemHandler;
+import org.apache.commons.math3.ode.events.EventHandler;
+import org.apache.commons.math3.ode.sampling.StepHandler;
+import org.apache.commons.math3.ode.sampling.StepInterpolator;
+import org.apache.commons.math3.util.FastMath;
+import org.junit.Assert;
+import org.junit.Test;
+
+public class LutherIntegratorTest {
+
+    @Test
+    public void testMissedEndEvent()
+            throws DimensionMismatchException, NumberIsTooSmallException,
+            MaxCountExceededException, NoBracketingException {
+        final double   t0     = 1878250320.0000029;
+        final double   tEvent = 1878250379.9999986;
+        final double[] k      = { 1.0e-4, 1.0e-5, 1.0e-6 };
+        FirstOrderDifferentialEquations ode = new FirstOrderDifferentialEquations() {
+
+            public int getDimension() {
+                return k.length;
+            }
+
+            public void computeDerivatives(double t, double[] y, double[] yDot) {
+                for (int i = 0; i < y.length; ++i) {
+                    yDot[i] = k[i] * y[i];
+                }
+            }
+        };
+
+        LutherIntegrator integrator = new LutherIntegrator(60.0);
+
+        double[] y0   = new double[k.length];
+        for (int i = 0; i < y0.length; ++i) {
+            y0[i] = i + 1;
+        }
+        double[] y    = new double[k.length];
+
+        double finalT = integrator.integrate(ode, t0, y0, tEvent, y);
+        Assert.assertEquals(tEvent, finalT, 1.0e-15);
+        for (int i = 0; i < y.length; ++i) {
+            Assert.assertEquals(y0[i] * FastMath.exp(k[i] * (finalT - t0)), y[i], 1.0e-15);
+        }
+
+        integrator.addEventHandler(new EventHandler() {
+
+            public void init(double t0, double[] y0, double t) {
+            }
+
+            public void resetState(double t, double[] y) {
+            }
+
+            public double g(double t, double[] y) {
+                return t - tEvent;
+            }
+
+            public Action eventOccurred(double t, double[] y, boolean increasing) {
+                Assert.assertEquals(tEvent, t, 1.0e-15);
+                return Action.CONTINUE;
+            }
+        }, Double.POSITIVE_INFINITY, 1.0e-20, 100);
+        finalT = integrator.integrate(ode, t0, y0, tEvent + 120, y);
+        Assert.assertEquals(tEvent + 120, finalT, 1.0e-15);
+        for (int i = 0; i < y.length; ++i) {
+            Assert.assertEquals(y0[i] * FastMath.exp(k[i] * (finalT - t0)), y[i], 1.0e-15);
+        }
+
+    }
+
+    @Test
+    public void testSanityChecks()
+            throws DimensionMismatchException, NumberIsTooSmallException,
+            MaxCountExceededException, NoBracketingException {
+        try  {
+            TestProblem1 pb = new TestProblem1();
+            new LutherIntegrator(0.01).integrate(pb,
+                                                 0.0, new double[pb.getDimension()+10],
+                                                 1.0, new double[pb.getDimension()]);
+            Assert.fail("an exception should have been thrown");
+        } catch(DimensionMismatchException ie) {
+        }
+        try  {
+            TestProblem1 pb = new TestProblem1();
+            new LutherIntegrator(0.01).integrate(pb,
+                                                 0.0, new double[pb.getDimension()],
+                                                 1.0, new double[pb.getDimension()+10]);
+            Assert.fail("an exception should have been thrown");
+        } catch(DimensionMismatchException ie) {
+        }
+        try  {
+            TestProblem1 pb = new TestProblem1();
+            new LutherIntegrator(0.01).integrate(pb,
+                                                 0.0, new double[pb.getDimension()],
+                                                 0.0, new double[pb.getDimension()]);
+            Assert.fail("an exception should have been thrown");
+        } catch(NumberIsTooSmallException ie) {
+        }
+    }
+
+    @Test
+    public void testDecreasingSteps()
+            throws DimensionMismatchException, NumberIsTooSmallException,
+            MaxCountExceededException, NoBracketingException {
+
+        TestProblemAbstract[] problems = TestProblemFactory.getProblems();
+        for (int k = 0; k < problems.length; ++k) {
+
+            double previousValueError = Double.NaN;
+            double previousTimeError = Double.NaN;
+            for (int i = 4; i < 10; ++i) {
+
+                TestProblemAbstract pb = problems[k].copy();
+                double step = (pb.getFinalTime() - pb.getInitialTime()) * FastMath.pow(2.0,
-i);
+
+                FirstOrderIntegrator integ = new LutherIntegrator(step);
+                TestProblemHandler handler = new TestProblemHandler(pb, integ);
+                integ.addStepHandler(handler);
+                EventHandler[] functions = pb.getEventsHandlers();
+                for (int l = 0; l < functions.length; ++l) {
+                    integ.addEventHandler(functions[l],
+                                          Double.POSITIVE_INFINITY, 1.0e-6 * step, 1000);
+                }
+                Assert.assertEquals(functions.length, integ.getEventHandlers().size());
+                double stopTime = integ.integrate(pb, pb.getInitialTime(), pb.getInitialState(),
+                                                  pb.getFinalTime(), new double[pb.getDimension()]);
+                if (functions.length == 0) {
+                    Assert.assertEquals(pb.getFinalTime(), stopTime, 1.0e-10);
+                }
+
+                double error = handler.getMaximalValueError();
+                if (i > 4) {
+                    Assert.assertTrue(error < 1.01 * FastMath.abs(previousValueError));
+                }
+                previousValueError = error;
+
+                double timeError = handler.getMaximalTimeError();
+                if (i > 4) {
+                    Assert.assertTrue(timeError <= FastMath.abs(previousTimeError));
+                }
+                previousTimeError = timeError;
+
+                integ.clearEventHandlers();
+                Assert.assertEquals(0, integ.getEventHandlers().size());
+            }
+
+        }
+
+    }
+
+    @Test
+    public void testSmallStep()
+            throws DimensionMismatchException, NumberIsTooSmallException,
+            MaxCountExceededException, NoBracketingException {
+
+        TestProblem1 pb = new TestProblem1();
+        double step = (pb.getFinalTime() - pb.getInitialTime()) * 0.001;
+
+        FirstOrderIntegrator integ = new LutherIntegrator(step);
+        TestProblemHandler handler = new TestProblemHandler(pb, integ);
+        integ.addStepHandler(handler);
+        integ.integrate(pb, pb.getInitialTime(), pb.getInitialState(),
+                        pb.getFinalTime(), new double[pb.getDimension()]);
+
+        Assert.assertTrue(handler.getLastError() < 9.0e-17);
+        Assert.assertTrue(handler.getMaximalValueError() < 4.0e-15);
+        Assert.assertEquals(0, handler.getMaximalTimeError(), 1.0e-12);
+        Assert.assertEquals("Luther", integ.getName());
+    }
+
+    @Test
+    public void testBigStep()
+            throws DimensionMismatchException, NumberIsTooSmallException,
+            MaxCountExceededException, NoBracketingException {
+
+        TestProblem1 pb = new TestProblem1();
+        double step = (pb.getFinalTime() - pb.getInitialTime()) * 0.2;
+
+        FirstOrderIntegrator integ = new LutherIntegrator(step);
+        TestProblemHandler handler = new TestProblemHandler(pb, integ);
+        integ.addStepHandler(handler);
+        integ.integrate(pb, pb.getInitialTime(), pb.getInitialState(),
+                        pb.getFinalTime(), new double[pb.getDimension()]);
+
+        Assert.assertTrue(handler.getLastError() > 0.00002);
+        Assert.assertTrue(handler.getMaximalValueError() > 0.001);
+        Assert.assertEquals(0, handler.getMaximalTimeError(), 1.0e-12);
+
+    }
+
+    @Test
+    public void testBackward()
+            throws DimensionMismatchException, NumberIsTooSmallException,
+            MaxCountExceededException, NoBracketingException {
+
+        TestProblem5 pb = new TestProblem5();
+        double step = FastMath.abs(pb.getFinalTime() - pb.getInitialTime()) * 0.001;
+
+        FirstOrderIntegrator integ = new LutherIntegrator(step);
+        TestProblemHandler handler = new TestProblemHandler(pb, integ);
+        integ.addStepHandler(handler);
+        integ.integrate(pb, pb.getInitialTime(), pb.getInitialState(),
+                        pb.getFinalTime(), new double[pb.getDimension()]);
+
+        Assert.assertTrue(handler.getLastError() < 3.0e-13);
+        Assert.assertTrue(handler.getMaximalValueError() < 5.0e-13);
+        Assert.assertEquals(0, handler.getMaximalTimeError(), 1.0e-12);
+        Assert.assertEquals("Luther", integ.getName());
+    }
+
+    @Test
+    public void testKepler()
+            throws DimensionMismatchException, NumberIsTooSmallException,
+            MaxCountExceededException, NoBracketingException {
+
+        final TestProblem3 pb  = new TestProblem3(0.9);
+        double step = (pb.getFinalTime() - pb.getInitialTime()) * 0.0003;
+
+        FirstOrderIntegrator integ = new LutherIntegrator(step);
+        integ.addStepHandler(new KeplerHandler(pb));
+        integ.integrate(pb,
+                        pb.getInitialTime(), pb.getInitialState(),
+                        pb.getFinalTime(), new double[pb.getDimension()]);
+    }
+
+    private static class KeplerHandler implements StepHandler {
+        public KeplerHandler(TestProblem3 pb) {
+            this.pb = pb;
+            maxError = 0;
+        }
+        public void init(double t0, double[] y0, double t) {
+            maxError = 0;
+        }
+        public void handleStep(StepInterpolator interpolator, boolean isLast) {
+
+            double[] interpolatedY = interpolator.getInterpolatedState ();
+            double[] theoreticalY  = pb.computeTheoreticalState(interpolator.getCurrentTime());
+            double dx = interpolatedY[0] - theoreticalY[0];
+            double dy = interpolatedY[1] - theoreticalY[1];
+            double error = dx * dx + dy * dy;
+            if (error > maxError) {
+                maxError = error;
+            }
+            if (isLast) {
+                Assert.assertTrue(maxError < 2.2e-7);
+            }
+        }
+        private double maxError = 0;
+        private TestProblem3 pb;
+    }
+
+    @Test
+    public void testStepSize()
+            throws DimensionMismatchException, NumberIsTooSmallException,
+            MaxCountExceededException, NoBracketingException {
+        final double step = 1.23456;
+        FirstOrderIntegrator integ = new LutherIntegrator(step);
+        integ.addStepHandler(new StepHandler() {
+            public void handleStep(StepInterpolator interpolator, boolean isLast) {
+                if (! isLast) {
+                    Assert.assertEquals(step,
+                                        interpolator.getCurrentTime() - interpolator.getPreviousTime(),
+                                        1.0e-12);
+                }
+            }
+            public void init(double t0, double[] y0, double t) {
+            }
+        });
+        integ.integrate(new FirstOrderDifferentialEquations() {
+            public void computeDerivatives(double t, double[] y, double[] dot) {
+                dot[0] = 1.0;
+            }
+            public int getDimension() {
+                return 1;
+            }
+        }, 0.0, new double[] { 0.0 }, 5.0, new double[1]);
+    }
+
+}

Propchange: commons/proper/math/trunk/src/test/java/org/apache/commons/math3/ode/nonstiff/LutherIntegratorTest.java
------------------------------------------------------------------------------
    svn:eol-style = native

Propchange: commons/proper/math/trunk/src/test/java/org/apache/commons/math3/ode/nonstiff/LutherIntegratorTest.java
------------------------------------------------------------------------------
    svn:keywords = "Author Date Id Revision"

Added: commons/proper/math/trunk/src/test/java/org/apache/commons/math3/ode/nonstiff/LutherStepInterpolatorTest.java
URL: http://svn.apache.org/viewvc/commons/proper/math/trunk/src/test/java/org/apache/commons/math3/ode/nonstiff/LutherStepInterpolatorTest.java?rev=1588753&view=auto
==============================================================================
--- commons/proper/math/trunk/src/test/java/org/apache/commons/math3/ode/nonstiff/LutherStepInterpolatorTest.java
(added)
+++ commons/proper/math/trunk/src/test/java/org/apache/commons/math3/ode/nonstiff/LutherStepInterpolatorTest.java
Sun Apr 20 13:25:11 2014
@@ -0,0 +1,98 @@
+/*
+ * Licensed to the Apache Software Foundation (ASF) under one or more
+ * contributor license agreements.  See the NOTICE file distributed with
+ * this work for additional information regarding copyright ownership.
+ * The ASF licenses this file to You under the Apache License, Version 2.0
+ * (the "License"); you may not use this file except in compliance with
+ * the License.  You may obtain a copy of the License at
+ *
+ *      http://www.apache.org/licenses/LICENSE-2.0
+ *
+ * Unless required by applicable law or agreed to in writing, software
+ * distributed under the License is distributed on an "AS IS" BASIS,
+ * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+ * See the License for the specific language governing permissions and
+ * limitations under the License.
+ */
+
+package org.apache.commons.math3.ode.nonstiff;
+
+
+import java.io.ByteArrayInputStream;
+import java.io.ByteArrayOutputStream;
+import java.io.IOException;
+import java.io.ObjectInputStream;
+import java.io.ObjectOutputStream;
+import java.util.Random;
+
+import org.apache.commons.math3.exception.DimensionMismatchException;
+import org.apache.commons.math3.exception.MaxCountExceededException;
+import org.apache.commons.math3.exception.NoBracketingException;
+import org.apache.commons.math3.exception.NumberIsTooSmallException;
+import org.apache.commons.math3.ode.ContinuousOutputModel;
+import org.apache.commons.math3.ode.TestProblem3;
+import org.apache.commons.math3.ode.sampling.StepHandler;
+import org.apache.commons.math3.ode.sampling.StepInterpolatorTestUtils;
+import org.junit.Assert;
+import org.junit.Test;
+
+public class LutherStepInterpolatorTest {
+
+    @Test
+    public void derivativesConsistency()
+            throws DimensionMismatchException, NumberIsTooSmallException,
+            MaxCountExceededException, NoBracketingException {
+        TestProblem3 pb = new TestProblem3();
+        double step = (pb.getFinalTime() - pb.getInitialTime()) * 0.001;
+        LutherIntegrator integ = new LutherIntegrator(step);
+        StepInterpolatorTestUtils.checkDerivativesConsistency(integ, pb, 1.0e-10);
+    }
+
+    @Test
+    public void serialization()
+            throws IOException, ClassNotFoundException,
+            DimensionMismatchException, NumberIsTooSmallException,
+            MaxCountExceededException, NoBracketingException  {
+
+        TestProblem3 pb = new TestProblem3(0.9);
+        double step = (pb.getFinalTime() - pb.getInitialTime()) * 0.0003;
+        LutherIntegrator integ = new LutherIntegrator(step);
+        integ.addStepHandler(new ContinuousOutputModel());
+        integ.integrate(pb,
+                        pb.getInitialTime(), pb.getInitialState(),
+                        pb.getFinalTime(), new double[pb.getDimension()]);
+
+        ByteArrayOutputStream bos = new ByteArrayOutputStream();
+        ObjectOutputStream    oos = new ObjectOutputStream(bos);
+        for (StepHandler handler : integ.getStepHandlers()) {
+            oos.writeObject(handler);
+        }
+
+        Assert.assertTrue(bos.size() > 1200000);
+        Assert.assertTrue(bos.size() < 1250000);
+
+        ByteArrayInputStream  bis = new ByteArrayInputStream(bos.toByteArray());
+        ObjectInputStream     ois = new ObjectInputStream(bis);
+        ContinuousOutputModel cm  = (ContinuousOutputModel) ois.readObject();
+
+        Random random = new Random(347588535632l);
+        double maxError = 0.0;
+        for (int i = 0; i < 1000; ++i) {
+            double r = random.nextDouble();
+            double time = r * pb.getInitialTime() + (1.0 - r) * pb.getFinalTime();
+            cm.setInterpolatedTime(time);
+            double[] interpolatedY = cm.getInterpolatedState ();
+            double[] theoreticalY  = pb.computeTheoreticalState(time);
+            double dx = interpolatedY[0] - theoreticalY[0];
+            double dy = interpolatedY[1] - theoreticalY[1];
+            double error = dx * dx + dy * dy;
+            if (error > maxError) {
+                maxError = error;
+            }
+        }
+
+        Assert.assertTrue(maxError < 2.2e-7);
+
+    }
+
+}

Propchange: commons/proper/math/trunk/src/test/java/org/apache/commons/math3/ode/nonstiff/LutherStepInterpolatorTest.java
------------------------------------------------------------------------------
    svn:eol-style = native

Propchange: commons/proper/math/trunk/src/test/java/org/apache/commons/math3/ode/nonstiff/LutherStepInterpolatorTest.java
------------------------------------------------------------------------------
    svn:keywords = "Author Date Id Revision"



Mime
View raw message