Repository: commons-math Updated Branches: refs/heads/field-ode b41ef1abd -> 23f3ca423
Field-based version of Luther method for solving ODE. Project: http://git-wip-us.apache.org/repos/asf/commons-math/repo Commit: http://git-wip-us.apache.org/repos/asf/commons-math/commit/23f3ca42 Tree: http://git-wip-us.apache.org/repos/asf/commons-math/tree/23f3ca42 Diff: http://git-wip-us.apache.org/repos/asf/commons-math/diff/23f3ca42 Branch: refs/heads/field-ode Commit: 23f3ca423ff7d79a814086b1d37b04d994e95cab Parents: b41ef1a Author: Luc Maisonobe <l...@apache.org> Authored: Sun Nov 15 14:03:03 2015 +0100 Committer: Luc Maisonobe <l...@apache.org> Committed: Sun Nov 15 14:03:03 2015 +0100 ---------------------------------------------------------------------- .../ode/nonstiff/LutherFieldIntegrator.java | 94 +++++++++ .../nonstiff/LutherFieldStepInterpolator.java | 208 +++++++++++++++++++ 2 files changed, 302 insertions(+) ---------------------------------------------------------------------- http://git-wip-us.apache.org/repos/asf/commons-math/blob/23f3ca42/src/main/java/org/apache/commons/math3/ode/nonstiff/LutherFieldIntegrator.java ---------------------------------------------------------------------- diff --git a/src/main/java/org/apache/commons/math3/ode/nonstiff/LutherFieldIntegrator.java b/src/main/java/org/apache/commons/math3/ode/nonstiff/LutherFieldIntegrator.java new file mode 100644 index 0000000..26156d1 --- /dev/null +++ b/src/main/java/org/apache/commons/math3/ode/nonstiff/LutherFieldIntegrator.java @@ -0,0 +1,94 @@ +/* + * 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.Field; +import org.apache.commons.math3.RealFieldElement; +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 EulerFieldIntegrator + * @see ClassicalRungeKuttaFieldIntegrator + * @see GillFieldIntegrator + * @see MidpointFieldIntegrator + * @see ThreeEighthesFieldIntegrator + * @param <T> the type of the field elements + * @since 3.6 + */ + +public class LutherFieldIntegrator<T extends RealFieldElement<T>> + extends RungeKuttaFieldIntegrator<T> { + + /** 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 field field to which the time and state vector elements belong + * @param step integration step + */ + public LutherFieldIntegrator(final Field<T> field, final T step) { + super(field, "Luther", STATIC_C, STATIC_A, STATIC_B, new LutherFieldStepInterpolator<T>(), step); + } + +} http://git-wip-us.apache.org/repos/asf/commons-math/blob/23f3ca42/src/main/java/org/apache/commons/math3/ode/nonstiff/LutherFieldStepInterpolator.java ---------------------------------------------------------------------- diff --git a/src/main/java/org/apache/commons/math3/ode/nonstiff/LutherFieldStepInterpolator.java b/src/main/java/org/apache/commons/math3/ode/nonstiff/LutherFieldStepInterpolator.java new file mode 100644 index 0000000..0647419 --- /dev/null +++ b/src/main/java/org/apache/commons/math3/ode/nonstiff/LutherFieldStepInterpolator.java @@ -0,0 +1,208 @@ +/* + * 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.RealFieldElement; +import org.apache.commons.math3.ode.FieldEquationsMapper; +import org.apache.commons.math3.ode.FieldODEStateAndDerivative; +import org.apache.commons.math3.util.FastMath; +import org.apache.commons.math3.util.MathArrays; + +/** + * This class represents an interpolator over the last step during an + * ODE integration for the 6th order Luther integrator. + * + * <p>This interpolator computes dense output inside the last + * step computed. The interpolation equation is consistent with the + * integration scheme.</p> + * + * @see LutherFieldIntegrator + * @param <T> the type of the field elements + * @since 3.6 + */ + +class LutherFieldStepInterpolator<T extends RealFieldElement<T>> + extends RungeKuttaFieldStepInterpolator<T> { + + /** 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.AbstractFieldStepInterpolator#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 + * RungeKuttaFieldIntegrator} class uses the prototyping design pattern + * to create the step interpolators by cloning an uninitialized model + * and later initializing the copy. + */ + LutherFieldStepInterpolator() { + } + + /** 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 + */ + LutherFieldStepInterpolator(final LutherFieldStepInterpolator<T> interpolator) { + super(interpolator); + } + + /** {@inheritDoc} */ + @Override + protected LutherFieldStepInterpolator<T> doCopy() { + return new LutherFieldStepInterpolator<T>(this); + } + + + /** {@inheritDoc} */ + @Override + protected FieldODEStateAndDerivative<T> computeInterpolatedStateAndDerivatives(final FieldEquationsMapper<T> mapper, + final T time, final T theta, + final T 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 T coeffDot1 = theta.multiply(theta.multiply(theta.multiply(theta.multiply( 21 ).add( -47 )).add( 36 )).add( -54 / 5.0)).add(1); + // not really needed as it is zero: final T coeffDot2 = theta.getField().getZero(); + final T coeffDot3 = theta.multiply(theta.multiply(theta.multiply(theta.multiply( 112 ).add(-608 / 3.0)).add( 320 / 3.0 )).add(-208 / 15.0)); + final T coeffDot4 = theta.multiply(theta.multiply(theta.multiply(theta.multiply( -567 / 5.0).add( 972 / 5.0)).add( -486 / 5.0 )).add( 324 / 25.0)); + final T coeffDot5 = theta.multiply(theta.multiply(theta.multiply(theta.multiply(( -49 - 49 * Q) / 5.0).add((392 + 287 * Q) / 15.0)).add((-637 - 357 * Q) / 30.0)).add((833 + 343 * Q) / 150.0)); + final T coeffDot6 = theta.multiply(theta.multiply(theta.multiply(theta.multiply(( -49 + 49 * Q) / 5.0).add((392 - 287 * Q) / 15.0)).add((-637 + 357 * Q) / 30.0)).add((833 - 343 * Q) / 150.0)); + final T coeffDot7 = theta.multiply(theta.multiply(theta.multiply( 3 ).add( -3 ).add( 3 / 5.0))); + final T[] interpolatedState = MathArrays.buildArray(theta.getField(), previousState.length); + final T[] interpolatedDerivatives = MathArrays.buildArray(theta.getField(), previousState.length); + + if ((previousState != null) && (theta.getReal() <= 0.5)) { + + final T coeff1 = theta.multiply(theta.multiply(theta.multiply(theta.multiply( 21 / 5.0).add( -47 / 4.0)).add( 12 )).add( -27 / 5.0)).add(1); + // not really needed as it is zero: final T coeff2 = theta.getField().getZero(); + final T coeff3 = theta.multiply(theta.multiply(theta.multiply(theta.multiply( 112 / 5.0).add(-152 / 3.0)).add( 320 / 9.0 )).add(-104 / 15.0)); + final T coeff4 = theta.multiply(theta.multiply(theta.multiply(theta.multiply(-567 / 25.0).add( 243 / 5.0)).add( -162 / 5.0 )).add( 162 / 25.0)); + final T coeff5 = theta.multiply(theta.multiply(theta.multiply(theta.multiply((-49 - 49 * Q) / 25.0).add((392 + 287 * Q) / 60.0)).add((-637 - 357 * Q) / 90.0)).add((833 + 343 * Q) / 300.0)); + final T coeff6 = theta.multiply(theta.multiply(theta.multiply(theta.multiply((-49 + 49 * Q) / 25.0).add((392 - 287 * Q) / 60.0)).add((-637 + 357 * Q) / 90.0)).add((833 - 343 * Q) / 300.0)); + final T coeff7 = theta.multiply(theta.multiply(theta.multiply( 3 / 4.0 ).add( -1 )).add( 3 / 10.0)); + for (int i = 0; i < interpolatedState.length; ++i) { + final T yDot1 = yDotK[0][i]; + // not really needed as associated coefficients are zero: final T yDot2 = yDotK[1][i]; + final T yDot3 = yDotK[2][i]; + final T yDot4 = yDotK[3][i]; + final T yDot5 = yDotK[4][i]; + final T yDot6 = yDotK[5][i]; + final T yDot7 = yDotK[6][i]; + interpolatedState[i] = previousState[i]. + add(theta.multiply(h). + multiply( coeff1.multiply(yDot1). + // not really needed as it is zero: add(coeff2.multiply(yDot2)). + add(coeff3.multiply(yDot3)). + add(coeff4.multiply(yDot4)). + add(coeff5.multiply(yDot5)). + add(coeff6.multiply(yDot6)). + add(coeff7.multiply(yDot7)))); + interpolatedDerivatives[i] = coeffDot1.multiply(yDot1). + // not really needed as it is zero: add(coeffDot2.multiply(yDot2)). + add(coeffDot3.multiply(yDot3)). + add(coeffDot4.multiply(yDot4)). + add(coeffDot5.multiply(yDot5)). + add(coeffDot6.multiply(yDot6)). + add(coeffDot7.multiply(yDot7)); + } + } else { + + final T coeff1 = theta.multiply(theta.multiply(theta.multiply(theta.multiply( -21 / 5.0).add( 151 / 20.0)).add( -89 / 20.0)).add( 19 / 20.0)).add( -1 / 20.0); + // not really needed as it is zero: final T coeff2 = theta.getField().getZero(); + final T coeff3 = theta.multiply(theta.multiply(theta.multiply(theta.multiply(-112 / 5.0).add( 424 / 15.0)).add( -328 / 45.0)).add( -16 / 45.0)).add(-16 / 45.0); + final T coeff4 = theta.multiply(theta.multiply(theta.multiply(theta.multiply( 567 / 25.0).add( -648 / 25.0)).add( 162 / 25.0))); + final T coeff5 = theta.multiply(theta.multiply(theta.multiply(theta.multiply(( 49 + 49 * Q) / 25.0).add((-1372 - 847 * Q) / 300.0)).add((2254 + 1029 * Q) / 900.0)).add( -49 / 180.0)).add(-49 / 180.0); + final T coeff6 = theta.multiply(theta.multiply(theta.multiply(theta.multiply(( 49 - 49 * Q) / 25.0).add((-1372 + 847 * Q) / 300.0)).add((2254 - 1029 * Q) / 900.0)).add( -49 / 180.0)).add(-49 / 180.0); + final T coeff7 = theta.multiply(theta.multiply(theta.multiply( -3 / 4.0 ).add( 1 / 4.0)).add( -1 / 20.0)).add( -1 / 20.0); + for (int i = 0; i < interpolatedState.length; ++i) { + final T yDot1 = yDotK[0][i]; + // not really needed as associated coefficients are zero: final T yDot2 = yDotK[1][i]; + final T yDot3 = yDotK[2][i]; + final T yDot4 = yDotK[3][i]; + final T yDot5 = yDotK[4][i]; + final T yDot6 = yDotK[5][i]; + final T yDot7 = yDotK[6][i]; + interpolatedState[i] = currentState[i]. + add(oneMinusThetaH. + multiply( coeff1.multiply(yDot1). + // not really needed as it is zero: add(coeff2.multiply(yDot2)). + add(coeff3.multiply(yDot3)). + add(coeff4.multiply(yDot4)). + add(coeff5.multiply(yDot5)). + add(coeff6.multiply(yDot6)). + add(coeff7.multiply(yDot7)))); + interpolatedDerivatives[i] = coeffDot1.multiply(yDot1). + // not really needed as it is zero: add(coeffDot2.multiply(yDot2)). + add(coeffDot3.multiply(yDot3)). + add(coeffDot4.multiply(yDot4)). + add(coeffDot5.multiply(yDot5)). + add(coeffDot6.multiply(yDot6)). + add(coeffDot7.multiply(yDot7)); + } + } + + return new FieldODEStateAndDerivative<T>(time, interpolatedState, yDotK[0]); + + } + +}
