diff options
Diffstat (limited to 'src/main/java/org/apache/commons/math3/ode/nonstiff/DormandPrince54StepInterpolator.java')
| -rw-r--r-- | src/main/java/org/apache/commons/math3/ode/nonstiff/DormandPrince54StepInterpolator.java | 225 |
1 files changed, 225 insertions, 0 deletions
diff --git a/src/main/java/org/apache/commons/math3/ode/nonstiff/DormandPrince54StepInterpolator.java b/src/main/java/org/apache/commons/math3/ode/nonstiff/DormandPrince54StepInterpolator.java new file mode 100644 index 0000000..04aa970 --- /dev/null +++ b/src/main/java/org/apache/commons/math3/ode/nonstiff/DormandPrince54StepInterpolator.java @@ -0,0 +1,225 @@ +/* + * 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.AbstractIntegrator; +import org.apache.commons.math3.ode.EquationsMapper; +import org.apache.commons.math3.ode.sampling.StepInterpolator; + +/** + * This class represents an interpolator over the last step during an + * ODE integration for the 5(4) Dormand-Prince integrator. + * + * @see DormandPrince54Integrator + * + * @since 1.2 + */ + +class DormandPrince54StepInterpolator + extends RungeKuttaStepInterpolator { + + /** 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; + + /** Serializable version identifier. */ + private static final long serialVersionUID = 20111120L; + + /** First vector for interpolation. */ + private double[] v1; + + /** Second vector for interpolation. */ + private double[] v2; + + /** Third vector for interpolation. */ + private double[] v3; + + /** Fourth vector for interpolation. */ + private double[] 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. + */ + // CHECKSTYLE: stop RedundantModifier + // the public modifier here is needed for serialization + public DormandPrince54StepInterpolator() { + super(); + v1 = null; + v2 = null; + v3 = null; + v4 = null; + vectorsInitialized = false; + } + // CHECKSTYLE: resume RedundantModifier + + /** 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 + */ + DormandPrince54StepInterpolator(final DormandPrince54StepInterpolator 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 StepInterpolator doCopy() { + return new DormandPrince54StepInterpolator(this); + } + + + /** {@inheritDoc} */ + @Override + public void reinitialize(final AbstractIntegrator integrator, + final double[] y, final double[][] yDotK, final boolean forward, + final EquationsMapper primaryMapper, + final EquationsMapper[] secondaryMappers) { + super.reinitialize(integrator, y, yDotK, forward, primaryMapper, secondaryMappers); + v1 = null; + v2 = null; + v3 = null; + v4 = null; + vectorsInitialized = false; + } + + /** {@inheritDoc} */ + @Override + public void storeTime(final double t) { + super.storeTime(t); + vectorsInitialized = false; + } + + /** {@inheritDoc} */ + @Override + protected void computeInterpolatedStateAndDerivatives(final double theta, + final double oneMinusThetaH) { + + if (! vectorsInitialized) { + + if (v1 == null) { + v1 = new double[interpolatedState.length]; + v2 = new double[interpolatedState.length]; + v3 = new double[interpolatedState.length]; + v4 = new double[interpolatedState.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 < interpolatedState.length; ++i) { + final double yDot0 = yDotK[0][i]; + final double yDot2 = yDotK[2][i]; + final double yDot3 = yDotK[3][i]; + final double yDot4 = yDotK[4][i]; + final double yDot5 = yDotK[5][i]; + final double yDot6 = yDotK[6][i]; + v1[i] = A70 * yDot0 + A72 * yDot2 + A73 * yDot3 + A74 * yDot4 + A75 * yDot5; + v2[i] = yDot0 - v1[i]; + v3[i] = v1[i] - v2[i] - yDot6; + v4[i] = D0 * yDot0 + D2 * yDot2 + D3 * yDot3 + D4 * yDot4 + D5 * yDot5 + D6 * yDot6; + } + + vectorsInitialized = true; + + } + + // interpolate + final double eta = 1 - theta; + final double twoTheta = 2 * theta; + final double dot2 = 1 - twoTheta; + final double dot3 = theta * (2 - 3 * theta); + final double dot4 = twoTheta * (1 + theta * (twoTheta - 3)); + if ((previousState != null) && (theta <= 0.5)) { + for (int i = 0; i < interpolatedState.length; ++i) { + interpolatedState[i] = + previousState[i] + theta * h * (v1[i] + eta * (v2[i] + theta * (v3[i] + eta * v4[i]))); + interpolatedDerivatives[i] = v1[i] + dot2 * v2[i] + dot3 * v3[i] + dot4 * v4[i]; + } + } else { + for (int i = 0; i < interpolatedState.length; ++i) { + interpolatedState[i] = + currentState[i] - oneMinusThetaH * (v1[i] - theta * (v2[i] + theta * (v3[i] + eta * v4[i]))); + interpolatedDerivatives[i] = v1[i] + dot2 * v2[i] + dot3 * v3[i] + dot4 * v4[i]; + } + } + + } + +} |