From 7f162008a2eea123a43f241e5e214dadf6b9e88e Mon Sep 17 00:00:00 2001 From: Luc Maisonobe Date: Sun, 20 Apr 2014 13:25:11 +0000 Subject: [PATCH] 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. git-svn-id: https://svn.apache.org/repos/asf/commons/proper/math/trunk@1588753 13f79535-47bb-0310-9956-ffa450edef68 --- src/changes/changes.xml | 3 + .../math3/ode/nonstiff/LutherIntegrator.java | 90 +++++ .../ode/nonstiff/LutherStepInterpolator.java | 180 ++++++++++ .../commons/math3/ode/package-info.java | 1 + src/site/xdoc/userguide/ode.xml | 1 + .../ode/nonstiff/LutherIntegratorTest.java | 309 ++++++++++++++++++ .../nonstiff/LutherStepInterpolatorTest.java | 98 ++++++ 7 files changed, 682 insertions(+) create mode 100644 src/main/java/org/apache/commons/math3/ode/nonstiff/LutherIntegrator.java create mode 100644 src/main/java/org/apache/commons/math3/ode/nonstiff/LutherStepInterpolator.java create mode 100644 src/test/java/org/apache/commons/math3/ode/nonstiff/LutherIntegratorTest.java create mode 100644 src/test/java/org/apache/commons/math3/ode/nonstiff/LutherStepInterpolatorTest.java diff --git a/src/changes/changes.xml b/src/changes/changes.xml index 18e1dd20b..91c8703a6 100644 --- a/src/changes/changes.xml +++ b/src/changes/changes.xml @@ -56,6 +56,9 @@ If the output is not quite correct, check for invisible trailing spaces! Added new methods for testing floating-point equality between the real (resp. imaginary) parts of two complex numbers. + + Added an order 6 fixed-step ODE integrator designed by H. A. Luther in 1968. + Bracketing utility for univariate root solvers returns a tighter interval than before. It also allows choosing the search interval expansion rate, supporting both linear diff --git a/src/main/java/org/apache/commons/math3/ode/nonstiff/LutherIntegrator.java b/src/main/java/org/apache/commons/math3/ode/nonstiff/LutherIntegrator.java new file mode 100644 index 000000000..1e1cf7067 --- /dev/null +++ b/src/main/java/org/apache/commons/math3/ode/nonstiff/LutherIntegrator.java @@ -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. + + *

+ * This method is described in H. A. Luther 1968 paper + * An explicit Sixth-Order Runge-Kutta Formula. + *

+ + *

This method is an explicit Runge-Kutta method, its Butcher-array + * is the following one : + * 

+ *        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
+ *
+ * where q = √21

+ * + * @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); + } + +} diff --git a/src/main/java/org/apache/commons/math3/ode/nonstiff/LutherStepInterpolator.java b/src/main/java/org/apache/commons/math3/ode/nonstiff/LutherStepInterpolator.java new file mode 100644 index 000000000..5019aad2e --- /dev/null +++ b/src/main/java/org/apache/commons/math3/ode/nonstiff/LutherStepInterpolator.java @@ -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. + * + *

This interpolator computes dense output inside the last + * step computed. The interpolation equation is consistent with the + * integration scheme.

+ * + * @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 + // + + // the method is implemented in the rkcheck tool + // . + // 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; + } + } + + } + +} diff --git a/src/main/java/org/apache/commons/math3/ode/package-info.java b/src/main/java/org/apache/commons/math3/ode/package-info.java index 615cc7d8d..c40e764d0 100644 --- a/src/main/java/org/apache/commons/math3/ode/package-info.java +++ b/src/main/java/org/apache/commons/math3/ode/package-info.java @@ -136,6 +136,7 @@ * {@link org.apache.commons.math3.ode.nonstiff.ClassicalRungeKuttaIntegrator Classical Runge-Kutta}4 * {@link org.apache.commons.math3.ode.nonstiff.GillIntegrator Gill}4 * {@link org.apache.commons.math3.ode.nonstiff.ThreeEighthesIntegrator 3/8}4 + * {@link org.apache.commons.math3.ode.nonstiff.LutherIntegrator Luther}6 * *

* diff --git a/src/site/xdoc/userguide/ode.xml b/src/site/xdoc/userguide/ode.xml index bf8a19a02..599eb8ff9 100644 --- a/src/site/xdoc/userguide/ode.xml +++ b/src/site/xdoc/userguide/ode.xml @@ -265,6 +265,7 @@ public int eventOccurred(double t, double[] y, boolean increasing) { Classical Runge-Kutta4 Gill4 3/84 + Luther6

diff --git a/src/test/java/org/apache/commons/math3/ode/nonstiff/LutherIntegratorTest.java b/src/test/java/org/apache/commons/math3/ode/nonstiff/LutherIntegratorTest.java new file mode 100644 index 000000000..fb8cbcaa3 --- /dev/null +++ b/src/test/java/org/apache/commons/math3/ode/nonstiff/LutherIntegratorTest.java @@ -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]); + } + +} diff --git a/src/test/java/org/apache/commons/math3/ode/nonstiff/LutherStepInterpolatorTest.java b/src/test/java/org/apache/commons/math3/ode/nonstiff/LutherStepInterpolatorTest.java new file mode 100644 index 000000000..d29fe20c2 --- /dev/null +++ b/src/test/java/org/apache/commons/math3/ode/nonstiff/LutherStepInterpolatorTest.java @@ -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); + + } + +}