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
This commit is contained in:
Luc Maisonobe
parent
4b724bdc2e
commit
7f162008a2
|
|
@ -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
|
Added new methods for testing floating-point equality between the real
|
||||||
(resp. imaginary) parts of two complex numbers.
|
(resp. imaginary) parts of two complex numbers.
|
||||||
</action>
|
</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" >
|
<action dev="luc" type="update" >
|
||||||
Bracketing utility for univariate root solvers returns a tighter interval than before.
|
Bracketing utility for univariate root solvers returns a tighter interval than before.
|
||||||
It also allows choosing the search interval expansion rate, supporting both linear
|
It also allows choosing the search interval expansion rate, supporting both linear
|
||||||
|
|
|
||||||
|
|
@ -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 = √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);
|
||||||
|
}
|
||||||
|
|
||||||
|
}
|
||||||
|
|
@ -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;
|
||||||
|
}
|
||||||
|
}
|
||||||
|
|
||||||
|
}
|
||||||
|
|
||||||
|
}
|
||||||
|
|
@ -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.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.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.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>
|
* </table>
|
||||||
* </p>
|
* </p>
|
||||||
*
|
*
|
||||||
|
|
|
||||||
|
|
@ -265,6 +265,7 @@ public int eventOccurred(double t, double[] y, boolean increasing) {
|
||||||
<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/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/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/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>
|
</table>
|
||||||
</p>
|
</p>
|
||||||
<p>
|
<p>
|
||||||
|
|
|
||||||
|
|
@ -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]);
|
||||||
|
}
|
||||||
|
|
||||||
|
}
|
||||||
|
|
@ -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);
|
||||||
|
|
||||||
|
}
|
||||||
|
|
||||||
|
}
|
||||||
Loading…
Reference in New Issue