Field-based version of Dormand-Prince 5(4) method for solving ODE.
This commit is contained in:
Luc Maisonobe
parent
301b0a8110
commit
e7a46ac6ca
|
|
@ -0,0 +1,172 @@
|
|||
/*
|
||||
* 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;
|
||||
import org.apache.commons.math4.util.MathUtils;
|
||||
|
||||
|
||||
/**
|
||||
* This class implements the 5(4) Dormand-Prince integrator for Ordinary
|
||||
* Differential Equations.
|
||||
|
||||
* <p>This integrator is an embedded Runge-Kutta integrator
|
||||
* of order 5(4) used in local extrapolation mode (i.e. the solution
|
||||
* is computed using the high order formula) with stepsize control
|
||||
* (and automatic step initialization) and continuous output. This
|
||||
* method uses 7 functions evaluations per step. However, since this
|
||||
* is an <i>fsal</i>, the last evaluation of one step is the same as
|
||||
* the first evaluation of the next step and hence can be avoided. So
|
||||
* the cost is really 6 functions evaluations per step.</p>
|
||||
*
|
||||
* <p>This method has been published (whithout the continuous output
|
||||
* that was added by Shampine in 1986) in the following article :
|
||||
* <pre>
|
||||
* A family of embedded Runge-Kutta formulae
|
||||
* J. R. Dormand and P. J. Prince
|
||||
* Journal of Computational and Applied Mathematics
|
||||
* volume 6, no 1, 1980, pp. 19-26
|
||||
* </pre></p>
|
||||
*
|
||||
* @param <T> the type of the field elements
|
||||
* @since 3.6
|
||||
*/
|
||||
|
||||
public class DormandPrince54FieldIntegrator<T extends RealFieldElement<T>>
|
||||
extends EmbeddedRungeKuttaFieldIntegrator<T> {
|
||||
|
||||
/** Integrator method name. */
|
||||
private static final String METHOD_NAME = "Dormand-Prince 5(4)";
|
||||
|
||||
/** Time steps Butcher array. */
|
||||
private static final double[] STATIC_C = {
|
||||
1.0/5.0, 3.0/10.0, 4.0/5.0, 8.0/9.0, 1.0, 1.0
|
||||
};
|
||||
|
||||
/** Internal weights Butcher array. */
|
||||
private static final double[][] STATIC_A = {
|
||||
{1.0/5.0},
|
||||
{3.0/40.0, 9.0/40.0},
|
||||
{44.0/45.0, -56.0/15.0, 32.0/9.0},
|
||||
{19372.0/6561.0, -25360.0/2187.0, 64448.0/6561.0, -212.0/729.0},
|
||||
{9017.0/3168.0, -355.0/33.0, 46732.0/5247.0, 49.0/176.0, -5103.0/18656.0},
|
||||
{35.0/384.0, 0.0, 500.0/1113.0, 125.0/192.0, -2187.0/6784.0, 11.0/84.0}
|
||||
};
|
||||
|
||||
/** Propagation weights Butcher array. */
|
||||
private static final double[] STATIC_B = {
|
||||
35.0/384.0, 0.0, 500.0/1113.0, 125.0/192.0, -2187.0/6784.0, 11.0/84.0, 0.0
|
||||
};
|
||||
|
||||
/** Error array, element 1. */
|
||||
private static final double E1 = 71.0 / 57600.0;
|
||||
|
||||
// element 2 is zero, so it is neither stored nor used
|
||||
|
||||
/** Error array, element 3. */
|
||||
private static final double E3 = -71.0 / 16695.0;
|
||||
|
||||
/** Error array, element 4. */
|
||||
private static final double E4 = 71.0 / 1920.0;
|
||||
|
||||
/** Error array, element 5. */
|
||||
private static final double E5 = -17253.0 / 339200.0;
|
||||
|
||||
/** Error array, element 6. */
|
||||
private static final double E6 = 22.0 / 525.0;
|
||||
|
||||
/** Error array, element 7. */
|
||||
private static final double E7 = -1.0 / 40.0;
|
||||
|
||||
/** Simple constructor.
|
||||
* Build a fifth order Dormand-Prince integrator with the given step bounds
|
||||
* @param field field to which the time and state vector elements belong
|
||||
* @param minStep minimal step (sign is irrelevant, regardless of
|
||||
* integration direction, forward or backward), the last step can
|
||||
* be smaller than this
|
||||
* @param maxStep maximal step (sign is irrelevant, regardless of
|
||||
* integration direction, forward or backward), the last step can
|
||||
* be smaller than this
|
||||
* @param scalAbsoluteTolerance allowed absolute error
|
||||
* @param scalRelativeTolerance allowed relative error
|
||||
*/
|
||||
public DormandPrince54FieldIntegrator(final Field<T> field,
|
||||
final double minStep, final double maxStep,
|
||||
final double scalAbsoluteTolerance,
|
||||
final double scalRelativeTolerance) {
|
||||
super(field, METHOD_NAME, true, STATIC_C, STATIC_A, STATIC_B,
|
||||
new DormandPrince54FieldStepInterpolator<T>(),
|
||||
minStep, maxStep, scalAbsoluteTolerance, scalRelativeTolerance);
|
||||
}
|
||||
|
||||
/** Simple constructor.
|
||||
* Build a fifth order Dormand-Prince integrator with the given step bounds
|
||||
* @param field field to which the time and state vector elements belong
|
||||
* @param minStep minimal step (sign is irrelevant, regardless of
|
||||
* integration direction, forward or backward), the last step can
|
||||
* be smaller than this
|
||||
* @param maxStep maximal step (sign is irrelevant, regardless of
|
||||
* integration direction, forward or backward), the last step can
|
||||
* be smaller than this
|
||||
* @param vecAbsoluteTolerance allowed absolute error
|
||||
* @param vecRelativeTolerance allowed relative error
|
||||
*/
|
||||
public DormandPrince54FieldIntegrator(final Field<T> field,
|
||||
final double minStep, final double maxStep,
|
||||
final double[] vecAbsoluteTolerance,
|
||||
final double[] vecRelativeTolerance) {
|
||||
super(field, METHOD_NAME, true, STATIC_C, STATIC_A, STATIC_B,
|
||||
new DormandPrince54FieldStepInterpolator<T>(),
|
||||
minStep, maxStep, vecAbsoluteTolerance, vecRelativeTolerance);
|
||||
}
|
||||
|
||||
/** {@inheritDoc} */
|
||||
@Override
|
||||
public int getOrder() {
|
||||
return 5;
|
||||
}
|
||||
|
||||
/** {@inheritDoc} */
|
||||
@Override
|
||||
protected T estimateError(final T[][] yDotK, final T[] y0, final T[] y1, final T h) {
|
||||
|
||||
T error = getField().getZero();
|
||||
|
||||
for (int j = 0; j < mainSetDimension; ++j) {
|
||||
final T errSum = yDotK[0][j].multiply(E1).
|
||||
add(yDotK[2][j].multiply(E3)).
|
||||
add(yDotK[3][j].multiply(E4)).
|
||||
add(yDotK[4][j].multiply(E5)).
|
||||
add(yDotK[5][j].multiply(E6)).
|
||||
add(yDotK[6][j].multiply(E7));
|
||||
|
||||
final T yScale = MathUtils.max(y0[j].abs(), y1[j].abs());
|
||||
final T tol = (vecAbsoluteTolerance == null) ?
|
||||
yScale.multiply(scalRelativeTolerance).add(scalAbsoluteTolerance) :
|
||||
yScale.multiply(vecRelativeTolerance[j]).add(vecAbsoluteTolerance[j]);
|
||||
final T ratio = h.multiply(errSum).divide(tol);
|
||||
error = error.add(ratio.multiply(ratio));
|
||||
|
||||
}
|
||||
|
||||
return error.divide(mainSetDimension).sqrt();
|
||||
|
||||
}
|
||||
|
||||
}
|
||||
|
|
@ -0,0 +1,245 @@
|
|||
/*
|
||||
* 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 represents an interpolator over the last step during an
|
||||
* ODE integration for the 5(4) Dormand-Prince integrator.
|
||||
*
|
||||
* @see DormandPrince54Integrator
|
||||
*
|
||||
* @param <T> the type of the field elements
|
||||
* @since 3.6
|
||||
*/
|
||||
|
||||
class DormandPrince54FieldStepInterpolator<T extends RealFieldElement<T>>
|
||||
extends RungeKuttaFieldStepInterpolator<T> {
|
||||
|
||||
/** Last row of the Butcher-array internal weights, element 0. */
|
||||
private static final double A70 = 35.0 / 384.0;
|
||||
|
||||
// element 1 is zero, so it is neither stored nor used
|
||||
|
||||
/** Last row of the Butcher-array internal weights, element 2. */
|
||||
private static final double A72 = 500.0 / 1113.0;
|
||||
|
||||
/** Last row of the Butcher-array internal weights, element 3. */
|
||||
private static final double A73 = 125.0 / 192.0;
|
||||
|
||||
/** Last row of the Butcher-array internal weights, element 4. */
|
||||
private static final double A74 = -2187.0 / 6784.0;
|
||||
|
||||
/** Last row of the Butcher-array internal weights, element 5. */
|
||||
private static final double A75 = 11.0 / 84.0;
|
||||
|
||||
/** Shampine (1986) Dense output, element 0. */
|
||||
private static final double D0 = -12715105075.0 / 11282082432.0;
|
||||
|
||||
// element 1 is zero, so it is neither stored nor used
|
||||
|
||||
/** Shampine (1986) Dense output, element 2. */
|
||||
private static final double D2 = 87487479700.0 / 32700410799.0;
|
||||
|
||||
/** Shampine (1986) Dense output, element 3. */
|
||||
private static final double D3 = -10690763975.0 / 1880347072.0;
|
||||
|
||||
/** Shampine (1986) Dense output, element 4. */
|
||||
private static final double D4 = 701980252875.0 / 199316789632.0;
|
||||
|
||||
/** Shampine (1986) Dense output, element 5. */
|
||||
private static final double D5 = -1453857185.0 / 822651844.0;
|
||||
|
||||
/** Shampine (1986) Dense output, element 6. */
|
||||
private static final double D6 = 69997945.0 / 29380423.0;
|
||||
|
||||
/** First vector for interpolation. */
|
||||
private T[] v1;
|
||||
|
||||
/** Second vector for interpolation. */
|
||||
private T[] v2;
|
||||
|
||||
/** Third vector for interpolation. */
|
||||
private T[] v3;
|
||||
|
||||
/** Fourth vector for interpolation. */
|
||||
private T[] v4;
|
||||
|
||||
/** Initialization indicator for the interpolation vectors. */
|
||||
private boolean vectorsInitialized;
|
||||
|
||||
/** Simple constructor.
|
||||
* This constructor builds an instance that is not usable yet, the
|
||||
* {@link #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 EmbeddedRungeKuttaIntegrator} uses the
|
||||
* prototyping design pattern to create the step interpolators by
|
||||
* cloning an uninitialized model and latter initializing the copy.
|
||||
*/
|
||||
DormandPrince54FieldStepInterpolator() {
|
||||
super();
|
||||
v1 = null;
|
||||
v2 = null;
|
||||
v3 = null;
|
||||
v4 = null;
|
||||
vectorsInitialized = false;
|
||||
}
|
||||
|
||||
/** 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
|
||||
*/
|
||||
DormandPrince54FieldStepInterpolator(final DormandPrince54FieldStepInterpolator<T> interpolator) {
|
||||
|
||||
super(interpolator);
|
||||
|
||||
if (interpolator.v1 == null) {
|
||||
|
||||
v1 = null;
|
||||
v2 = null;
|
||||
v3 = null;
|
||||
v4 = null;
|
||||
vectorsInitialized = false;
|
||||
|
||||
} else {
|
||||
|
||||
v1 = interpolator.v1.clone();
|
||||
v2 = interpolator.v2.clone();
|
||||
v3 = interpolator.v3.clone();
|
||||
v4 = interpolator.v4.clone();
|
||||
vectorsInitialized = interpolator.vectorsInitialized;
|
||||
|
||||
}
|
||||
|
||||
}
|
||||
|
||||
/** {@inheritDoc} */
|
||||
@Override
|
||||
protected DormandPrince54FieldStepInterpolator<T> doCopy() {
|
||||
return new DormandPrince54FieldStepInterpolator<T>(this);
|
||||
}
|
||||
|
||||
|
||||
/** {@inheritDoc} */
|
||||
@Override
|
||||
protected void reinitialize(final T[] y, final boolean isForward, final FieldEquationsMapper<T> equationsMapper) {
|
||||
super.reinitialize(y, isForward, equationsMapper);
|
||||
v1 = null;
|
||||
v2 = null;
|
||||
v3 = null;
|
||||
v4 = null;
|
||||
vectorsInitialized = false;
|
||||
}
|
||||
|
||||
/** {@inheritDoc} */
|
||||
@Override
|
||||
public void storeState(final FieldODEStateAndDerivative<T> state) {
|
||||
super.storeState(state);
|
||||
vectorsInitialized = false;
|
||||
}
|
||||
|
||||
/** {@inheritDoc} */
|
||||
@Override
|
||||
protected FieldODEStateAndDerivative<T> computeInterpolatedStateAndDerivatives(final FieldEquationsMapper<T> mapper,
|
||||
final T time, final T theta,
|
||||
final T oneMinusThetaH) {
|
||||
|
||||
if (! vectorsInitialized) {
|
||||
|
||||
if (v1 == null) {
|
||||
v1 = MathArrays.buildArray(time.getField(), previousState.length);
|
||||
v2 = MathArrays.buildArray(time.getField(), previousState.length);
|
||||
v3 = MathArrays.buildArray(time.getField(), previousState.length);
|
||||
v4 = MathArrays.buildArray(time.getField(), previousState.length);
|
||||
}
|
||||
|
||||
// no step finalization is needed for this interpolator
|
||||
|
||||
// we need to compute the interpolation vectors for this time step
|
||||
for (int i = 0; i < previousState.length; ++i) {
|
||||
final T yDot0 = yDotK[0][i];
|
||||
final T yDot2 = yDotK[2][i];
|
||||
final T yDot3 = yDotK[3][i];
|
||||
final T yDot4 = yDotK[4][i];
|
||||
final T yDot5 = yDotK[5][i];
|
||||
final T yDot6 = yDotK[6][i];
|
||||
v1[i] = yDot0.multiply(A70).
|
||||
add(yDot2.multiply(A72)).
|
||||
add(yDot3.multiply(A73)).
|
||||
add(yDot4.multiply(A74)).
|
||||
add(yDot5.multiply(A75));
|
||||
v2[i] = yDot0.subtract(v1[i]);
|
||||
v3[i] = v1[i].subtract(v2[i]).subtract(yDot6);
|
||||
v4[i] = yDot0.multiply(D0).
|
||||
add(yDot2.multiply(D2)).
|
||||
add(yDot3.multiply(D3)).
|
||||
add(yDot4.multiply(D4)).
|
||||
add(yDot5.multiply(D5)).
|
||||
add(yDot6.multiply(D6));
|
||||
}
|
||||
|
||||
vectorsInitialized = true;
|
||||
|
||||
}
|
||||
|
||||
// interpolate
|
||||
final T one = theta.getField().getOne();
|
||||
final T eta = one.subtract(theta);
|
||||
final T twoTheta = theta.multiply(2);
|
||||
final T dot2 = one.subtract(twoTheta);
|
||||
final T dot3 = theta.multiply(theta.multiply(-3).add(2));
|
||||
final T dot4 = twoTheta.multiply(theta.multiply(twoTheta.subtract(3)).add(1));
|
||||
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)) {
|
||||
for (int i = 0; i < interpolatedState.length; ++i) {
|
||||
interpolatedState[i] = previousState[i].
|
||||
add(theta.multiply(h.multiply(v1[i].
|
||||
add(eta.multiply(v2[i].
|
||||
add(theta.multiply(v3[i].
|
||||
add(eta.multiply(v4[i])))))))));
|
||||
interpolatedDerivatives[i] = v1[i].
|
||||
add(dot2.multiply(v2[i])).
|
||||
add(dot3.multiply(v3[i])).
|
||||
add(dot4.multiply(v4[i]));
|
||||
}
|
||||
} else {
|
||||
for (int i = 0; i < interpolatedState.length; ++i) {
|
||||
interpolatedState[i] = currentState[i].
|
||||
subtract(oneMinusThetaH.multiply(v1[i].
|
||||
subtract(theta.multiply(v2[i].
|
||||
add(theta.multiply(v3[i].
|
||||
add(eta.multiply(v4[i]))))))));
|
||||
interpolatedDerivatives[i] = v1[i].
|
||||
add(dot2.multiply(v2[i])).
|
||||
add(dot3.multiply(v3[i])).
|
||||
add(dot4.multiply(v4[i]));
|
||||
}
|
||||
}
|
||||
|
||||
return new FieldODEStateAndDerivative<T>(time, interpolatedState, yDotK[0]);
|
||||
|
||||
}
|
||||
|
||||
}
|
||||
Loading…
Reference in New Issue