From d67f7e2af93a1f715ce14ee373b07bfae9389370 Mon Sep 17 00:00:00 2001 From: Luc Maisonobe Date: Wed, 6 Jan 2016 12:24:16 +0100 Subject: [PATCH] Field-based version of classical Runge-Kutta method for solving ODE. --- .../ClassicalRungeKuttaFieldIntegrator.java | 79 ++++++++++ ...ssicalRungeKuttaFieldStepInterpolator.java | 141 ++++++++++++++++++ 2 files changed, 220 insertions(+) create mode 100644 src/main/java/org/apache/commons/math4/ode/nonstiff/ClassicalRungeKuttaFieldIntegrator.java create mode 100644 src/main/java/org/apache/commons/math4/ode/nonstiff/ClassicalRungeKuttaFieldStepInterpolator.java diff --git a/src/main/java/org/apache/commons/math4/ode/nonstiff/ClassicalRungeKuttaFieldIntegrator.java b/src/main/java/org/apache/commons/math4/ode/nonstiff/ClassicalRungeKuttaFieldIntegrator.java new file mode 100644 index 000000000..b69f5e36b --- /dev/null +++ b/src/main/java/org/apache/commons/math4/ode/nonstiff/ClassicalRungeKuttaFieldIntegrator.java @@ -0,0 +1,79 @@ +/* + * 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.math4.ode.nonstiff; + +import org.apache.commons.math4.Field; +import org.apache.commons.math4.RealFieldElement; + +/** + * This class implements the classical fourth order Runge-Kutta + * integrator for Ordinary Differential Equations (it is the most + * often used Runge-Kutta method). + * + *

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

+ *    0  |  0    0    0    0
+ *   1/2 | 1/2   0    0    0
+ *   1/2 |  0   1/2   0    0
+ *    1  |  0    0    1    0
+ *       |--------------------
+ *       | 1/6  1/3  1/3  1/6
+ *
+ *

+ * + * @see EulerFieldIntegrator + * @see GillFieldIntegrator + * @see MidpointFieldIntegrator + * @see ThreeEighthesFieldIntegrator + * @see LutherFieldIntegrator + * @param the type of the field elements + * @since 3.6 + */ + +public class ClassicalRungeKuttaFieldIntegrator> + extends RungeKuttaFieldIntegrator { + + /** Time steps Butcher array. */ + private static final double[] STATIC_C = { + 1.0 / 2.0, 1.0 / 2.0, 1.0 + }; + + /** Internal weights Butcher array. */ + private static final double[][] STATIC_A = { + { 1.0 / 2.0 }, + { 0.0, 1.0 / 2.0 }, + { 0.0, 0.0, 1.0 } + }; + + /** Propagation weights Butcher array. */ + private static final double[] STATIC_B = { + 1.0 / 6.0, 1.0 / 3.0, 1.0 / 3.0, 1.0 / 6.0 + }; + + /** Simple constructor. + * Build a fourth-order Runge-Kutta integrator with the given step. + * @param field field to which the time and state vector elements belong + * @param step integration step + */ + public ClassicalRungeKuttaFieldIntegrator(final Field field, final T step) { + super(field, "classical Runge-Kutta", STATIC_C, STATIC_A, STATIC_B, + new ClassicalRungeKuttaFieldStepInterpolator(), step); + } + +} diff --git a/src/main/java/org/apache/commons/math4/ode/nonstiff/ClassicalRungeKuttaFieldStepInterpolator.java b/src/main/java/org/apache/commons/math4/ode/nonstiff/ClassicalRungeKuttaFieldStepInterpolator.java new file mode 100644 index 000000000..98c0d9e98 --- /dev/null +++ b/src/main/java/org/apache/commons/math4/ode/nonstiff/ClassicalRungeKuttaFieldStepInterpolator.java @@ -0,0 +1,141 @@ +/* + * 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.math4.ode.nonstiff; + +import org.apache.commons.math4.RealFieldElement; +import org.apache.commons.math4.ode.FieldEquationsMapper; +import org.apache.commons.math4.ode.FieldODEStateAndDerivative; +import org.apache.commons.math4.util.MathArrays; + +/** + * This class implements a step interpolator for the classical fourth + * order Runge-Kutta integrator. + * + *

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

    + *
  • Using reference point at step start:
    + * y(tn + θ h) = y (tn) + * + θ (h/6) [ (6 - 9 θ + 4 θ2) y'1 + * + ( 6 θ - 4 θ2) (y'2 + y'3) + * + ( -3 θ + 4 θ2) y'4 + * ] + *
    • + *
  • Using reference point at step end:
    + * y(tn + θ h) = y (tn + h) + * + (1 - θ) (h/6) [ (-4 θ^2 + 5 θ - 1) y'1 + * +(4 θ^2 - 2 θ - 2) (y'2 + y'3) + * -(4 θ^2 + θ + 1) y'4 + * ] + *
    • + *
+ *

+ * + * where θ belongs to [0 ; 1] and where y'1 to y'4 are the four + * evaluations of the derivatives already computed during the + * step.

+ * + * @see ClassicalRungeKuttaFieldIntegrator + * @param the type of the field elements + * @since 3.6 + */ + +class ClassicalRungeKuttaFieldStepInterpolator> + extends RungeKuttaFieldStepInterpolator { + + /** Simple constructor. + * This constructor builds an instance that is not usable yet, the + * {@link RungeKuttaFieldStepInterpolator#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 latter initializing + * the copy. + */ + ClassicalRungeKuttaFieldStepInterpolator() { + } + + /** 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 + */ + ClassicalRungeKuttaFieldStepInterpolator(final ClassicalRungeKuttaFieldStepInterpolator interpolator) { + super(interpolator); + } + + /** {@inheritDoc} */ + @Override + protected ClassicalRungeKuttaFieldStepInterpolator doCopy() { + return new ClassicalRungeKuttaFieldStepInterpolator(this); + } + + /** {@inheritDoc} */ + @Override + protected FieldODEStateAndDerivative computeInterpolatedStateAndDerivatives(final FieldEquationsMapper mapper, + final T time, final T theta, + final T oneMinusThetaH) { + + final T one = time.getField().getOne(); + final T oneMinusTheta = one.subtract(theta); + final T oneMinus2Theta = one.subtract(theta.multiply(2)); + final T coeffDot1 = oneMinusTheta.multiply(oneMinus2Theta); + final T coeffDot23 = theta.multiply(oneMinusTheta).multiply(2); + final T coeffDot4 = theta.multiply(oneMinus2Theta).negate(); + 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 fourTheta2 = theta.multiply(theta).multiply(4); + final T s = theta.multiply(h).divide(6.0); + final T coeff1 = s.multiply(fourTheta2.subtract(theta.multiply(9)).add(6)); + final T coeff23 = s.multiply(theta.multiply(6).subtract(fourTheta2)); + final T coeff4 = s.multiply(fourTheta2.subtract(theta.multiply(3))); + for (int i = 0; i < interpolatedState.length; ++i) { + final T yDot1 = yDotK[0][i]; + final T yDot23 = yDotK[1][i].add(yDotK[2][i]); + final T yDot4 = yDotK[3][i]; + interpolatedState[i] = + previousState[i].add(coeff1.multiply(yDot1)).add(coeff23.multiply(yDot23)).add(coeff4.multiply(yDot4)); + interpolatedDerivatives[i] = + coeffDot1.multiply(yDot1).add(coeffDot23.multiply(yDot23)).add(coeffDot4.multiply(yDot4)); + } + } else { + final T fourTheta = theta.multiply(4); + final T s = oneMinusThetaH.divide(6); + final T coeff1 = s.multiply(theta.multiply(fourTheta.negate().add(5)).subtract(1)); + final T coeff23 = s.multiply(theta.multiply(fourTheta.subtract(2)).subtract(2)); + final T coeff4 = s.multiply(theta.multiply(fourTheta.negate().subtract(1)).subtract(1)); + for (int i = 0; i < interpolatedState.length; ++i) { + final T yDot1 = yDotK[0][i]; + final T yDot23 = yDotK[1][i].add(yDotK[2][i]); + final T yDot4 = yDotK[3][i]; + interpolatedState[i] = + currentState[i].add(coeff1.multiply(yDot1)).add(coeff23.multiply(yDot23)).add(coeff4.multiply(yDot4)); + interpolatedDerivatives[i] = + coeffDot1.multiply(yDot1).add(coeffDot23.multiply(yDot23)).add(coeffDot4.multiply(yDot4)); + } + } + + return new FieldODEStateAndDerivative(time, interpolatedState, yDotK[0]); + + } + +}