DormandPrince54StepInterpolator.java

  1. /*
  2.  * Licensed to the Apache Software Foundation (ASF) under one or more
  3.  * contributor license agreements.  See the NOTICE file distributed with
  4.  * this work for additional information regarding copyright ownership.
  5.  * The ASF licenses this file to You under the Apache License, Version 2.0
  6.  * (the "License"); you may not use this file except in compliance with
  7.  * the License.  You may obtain a copy of the License at
  8.  *
  9.  *      http://www.apache.org/licenses/LICENSE-2.0
  10.  *
  11.  * Unless required by applicable law or agreed to in writing, software
  12.  * distributed under the License is distributed on an "AS IS" BASIS,
  13.  * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
  14.  * See the License for the specific language governing permissions and
  15.  * limitations under the License.
  16.  */

  18. package org.apache.commons.math4.legacy.ode.nonstiff;

  20. import org.apache.commons.math4.legacy.ode.AbstractIntegrator;
  21. import org.apache.commons.math4.legacy.ode.EquationsMapper;
  22. import org.apache.commons.math4.legacy.ode.sampling.StepInterpolator;

  24. /**
  25.  * This class represents an interpolator over the last step during an
  26.  * ODE integration for the 5(4) Dormand-Prince integrator.
  27.  *
  28.  * @see DormandPrince54Integrator
  29.  *
  30.  * @since 1.2
  31.  */

  33. class DormandPrince54StepInterpolator
  34.   extends RungeKuttaStepInterpolator {

  36.     /** Last row of the Butcher-array internal weights, element 0. */
  37.     private static final double A70 =    35.0 /  384.0;

  39.     // element 1 is zero, so it is neither stored nor used

  41.     /** Last row of the Butcher-array internal weights, element 2. */
  42.     private static final double A72 =   500.0 / 1113.0;

  44.     /** Last row of the Butcher-array internal weights, element 3. */
  45.     private static final double A73 =   125.0 /  192.0;

  47.     /** Last row of the Butcher-array internal weights, element 4. */
  48.     private static final double A74 = -2187.0 / 6784.0;

  50.     /** Last row of the Butcher-array internal weights, element 5. */
  51.     private static final double A75 =    11.0 /   84.0;

  53.     /** Shampine (1986) Dense output, element 0. */
  54.     private static final double D0 =  -12715105075.0 /  11282082432.0;

  56.     // element 1 is zero, so it is neither stored nor used

  58.     /** Shampine (1986) Dense output, element 2. */
  59.     private static final double D2 =   87487479700.0 /  32700410799.0;

  61.     /** Shampine (1986) Dense output, element 3. */
  62.     private static final double D3 =  -10690763975.0 /   1880347072.0;

  64.     /** Shampine (1986) Dense output, element 4. */
  65.     private static final double D4 =  701980252875.0 / 199316789632.0;

  67.     /** Shampine (1986) Dense output, element 5. */
  68.     private static final double D5 =   -1453857185.0 /    822651844.0;

  70.     /** Shampine (1986) Dense output, element 6. */
  71.     private static final double D6 =      69997945.0 /     29380423.0;

  73.     /** Serializable version identifier. */
  74.     private static final long serialVersionUID = 20111120L;

  76.     /** First vector for interpolation. */
  77.     private double[] v1;

  79.     /** Second vector for interpolation. */
  80.     private double[] v2;

  82.     /** Third vector for interpolation. */
  83.     private double[] v3;

  85.     /** Fourth vector for interpolation. */
  86.     private double[] v4;

  88.     /** Initialization indicator for the interpolation vectors. */
  89.     private boolean vectorsInitialized;

  91.   /** Simple constructor.
  92.    * This constructor builds an instance that is not usable yet, the
  93.    * {@link #reinitialize} method should be called before using the
  94.    * instance in order to initialize the internal arrays. This
  95.    * constructor is used only in order to delay the initialization in
  96.    * some cases. The {@link EmbeddedRungeKuttaIntegrator} uses the
  97.    * prototyping design pattern to create the step interpolators by
  98.    * cloning an uninitialized model and latter initializing the copy.
  99.    */
  100.   // CHECKSTYLE: stop RedundantModifier
  101.   // the public modifier here is needed for serialization
  102.   public DormandPrince54StepInterpolator() {
  103.     super();
  104.     v1 = null;
  105.     v2 = null;
  106.     v3 = null;
  107.     v4 = null;
  108.     vectorsInitialized = false;
  109.   }
  110.   // CHECKSTYLE: resume RedundantModifier

  112.   /** Copy constructor.
  113.    * @param interpolator interpolator to copy from. The copy is a deep
  114.    * copy: its arrays are separated from the original arrays of the
  115.    * instance
  116.    */
  117.   DormandPrince54StepInterpolator(final DormandPrince54StepInterpolator interpolator) {

  119.     super(interpolator);

  121.     if (interpolator.v1 == null) {

  123.       v1 = null;
  124.       v2 = null;
  125.       v3 = null;
  126.       v4 = null;
  127.       vectorsInitialized = false;
  128.     } else {

  130.       v1 = interpolator.v1.clone();
  131.       v2 = interpolator.v2.clone();
  132.       v3 = interpolator.v3.clone();
  133.       v4 = interpolator.v4.clone();
  134.       vectorsInitialized = interpolator.vectorsInitialized;
  135.     }
  136.   }

  138.   /** {@inheritDoc} */
  139.   @Override
  140.   protected StepInterpolator doCopy() {
  141.     return new DormandPrince54StepInterpolator(this);
  142.   }


  145.   /** {@inheritDoc} */
  146.   @Override
  147.   public void reinitialize(final AbstractIntegrator integrator,
  148.                            final double[] y, final double[][] yDotK, final boolean forward,
  149.                            final EquationsMapper primaryMapper,
  150.                            final EquationsMapper[] secondaryMappers) {
  151.     super.reinitialize(integrator, y, yDotK, forward, primaryMapper, secondaryMappers);
  152.     v1 = null;
  153.     v2 = null;
  154.     v3 = null;
  155.     v4 = null;
  156.     vectorsInitialized = false;
  157.   }

  159.   /** {@inheritDoc} */
  160.   @Override
  161.   public void storeTime(final double t) {
  162.     super.storeTime(t);
  163.     vectorsInitialized = false;
  164.   }

  166.   /** {@inheritDoc} */
  167.   @Override
  168.   protected void computeInterpolatedStateAndDerivatives(final double theta,
  169.                                           final double oneMinusThetaH) {

  171.     if (! vectorsInitialized) {

  173.       if (v1 == null) {
  174.         v1 = new double[interpolatedState.length];
  175.         v2 = new double[interpolatedState.length];
  176.         v3 = new double[interpolatedState.length];
  177.         v4 = new double[interpolatedState.length];
  178.       }

  180.       // no step finalization is needed for this interpolator

  182.       // we need to compute the interpolation vectors for this time step
  183.       for (int i = 0; i < interpolatedState.length; ++i) {
  184.           final double yDot0 = yDotK[0][i];
  185.           final double yDot2 = yDotK[2][i];
  186.           final double yDot3 = yDotK[3][i];
  187.           final double yDot4 = yDotK[4][i];
  188.           final double yDot5 = yDotK[5][i];
  189.           final double yDot6 = yDotK[6][i];
  190.           v1[i] = A70 * yDot0 + A72 * yDot2 + A73 * yDot3 + A74 * yDot4 + A75 * yDot5;
  191.           v2[i] = yDot0 - v1[i];
  192.           v3[i] = v1[i] - v2[i] - yDot6;
  193.           v4[i] = D0 * yDot0 + D2 * yDot2 + D3 * yDot3 + D4 * yDot4 + D5 * yDot5 + D6 * yDot6;
  194.       }

  196.       vectorsInitialized = true;
  197.     }

  199.     // interpolate
  200.     final double eta = 1 - theta;
  201.     final double twoTheta = 2 * theta;
  202.     final double dot2 = 1 - twoTheta;
  203.     final double dot3 = theta * (2 - 3 * theta);
  204.     final double dot4 = twoTheta * (1 + theta * (twoTheta - 3));
  205.     if (previousState != null && theta <= 0.5) {
  206.         for (int i = 0; i < interpolatedState.length; ++i) {
  207.             interpolatedState[i] =
  208.                     previousState[i] + theta * h * (v1[i] + eta * (v2[i] + theta * (v3[i] + eta * v4[i])));
  209.             interpolatedDerivatives[i] = v1[i] + dot2 * v2[i] + dot3 * v3[i] + dot4 * v4[i];
  210.         }
  211.     } else {
  212.         for (int i = 0; i < interpolatedState.length; ++i) {
  213.             interpolatedState[i] =
  214.                     currentState[i] - oneMinusThetaH * (v1[i] - theta * (v2[i] + theta * (v3[i] + eta * v4[i])));
  215.             interpolatedDerivatives[i] = v1[i] + dot2 * v2[i] + dot3 * v3[i] + dot4 * v4[i];
  216.         }
  217.     }
  218.   }
  219. }
Created with JaCoCo 0.8.10.202304240956Code Coverage Report for Apache Commons Math 4.0-SNAPSHOT