DormandPrince853StepInterpolator.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 java.io.IOException;
  21. import java.io.ObjectInput;
  22. import java.io.ObjectOutput;

  24. import org.apache.commons.math4.legacy.exception.MaxCountExceededException;
  25. import org.apache.commons.math4.legacy.ode.AbstractIntegrator;
  26. import org.apache.commons.math4.legacy.ode.EquationsMapper;
  27. import org.apache.commons.math4.legacy.ode.sampling.StepInterpolator;

  29. /**
  30.  * This class represents an interpolator over the last step during an
  31.  * ODE integration for the 8(5,3) Dormand-Prince integrator.
  32.  *
  33.  * @see DormandPrince853Integrator
  34.  *
  35.  * @since 1.2
  36.  */

  38. class DormandPrince853StepInterpolator
  39.   extends RungeKuttaStepInterpolator {

  41.     /** Serializable version identifier. */
  42.     private static final long serialVersionUID = 20111120L;

  44.     /** Propagation weights, element 1. */
  45.     private static final double B_01 =         104257.0 / 1920240.0;

  47.     // elements 2 to 5 are zero, so they are neither stored nor used

  49.     /** Propagation weights, element 6. */
  50.     private static final double B_06 =        3399327.0 / 763840.0;

  52.     /** Propagation weights, element 7. */
  53.     private static final double B_07 =       66578432.0 / 35198415.0;

  55.     /** Propagation weights, element 8. */
  56.     private static final double B_08 =    -1674902723.0 / 288716400.0;

  58.     /** Propagation weights, element 9. */
  59.     private static final double B_09 = 54980371265625.0 / 176692375811392.0;

  61.     /** Propagation weights, element 10. */
  62.     private static final double B_10 =        -734375.0 / 4826304.0;

  64.     /** Propagation weights, element 11. */
  65.     private static final double B_11 =      171414593.0 / 851261400.0;

  67.     /** Propagation weights, element 12. */
  68.     private static final double B_12 =         137909.0 / 3084480.0;

  70.     /** Time step for stage 14 (interpolation only). */
  71.     private static final double C14    = 1.0 / 10.0;

  73.     /** Internal weights for stage 14, element 1. */
  74.     private static final double K14_01 =       13481885573.0 / 240030000000.0      - B_01;

  76.     // elements 2 to 5 are zero, so they are neither stored nor used

  78.     /** Internal weights for stage 14, element 6. */
  79.     private static final double K14_06 =                 0.0                       - B_06;

  81.     /** Internal weights for stage 14, element 7. */
  82.     private static final double K14_07 =      139418837528.0 / 549975234375.0      - B_07;

  84.     /** Internal weights for stage 14, element 8. */
  85.     private static final double K14_08 =   -11108320068443.0 / 45111937500000.0    - B_08;

  87.     /** Internal weights for stage 14, element 9. */
  88.     private static final double K14_09 = -1769651421925959.0 / 14249385146080000.0 - B_09;

  90.     /** Internal weights for stage 14, element 10. */
  91.     private static final double K14_10 =          57799439.0 / 377055000.0         - B_10;

  93.     /** Internal weights for stage 14, element 11. */
  94.     private static final double K14_11 =      793322643029.0 / 96734250000000.0    - B_11;

  96.     /** Internal weights for stage 14, element 12. */
  97.     private static final double K14_12 =        1458939311.0 / 192780000000.0      - B_12;

  99.     /** Internal weights for stage 14, element 13. */
  100.     private static final double K14_13 =             -4149.0 / 500000.0;

  102.     /** Time step for stage 15 (interpolation only). */
  103.     private static final double C15    = 1.0 / 5.0;


  106.     /** Internal weights for stage 15, element 1. */
  107.     private static final double K15_01 =     1595561272731.0 / 50120273500000.0    - B_01;

  109.     // elements 2 to 5 are zero, so they are neither stored nor used

  111.     /** Internal weights for stage 15, element 6. */
  112.     private static final double K15_06 =      975183916491.0 / 34457688031250.0    - B_06;

  114.     /** Internal weights for stage 15, element 7. */
  115.     private static final double K15_07 =    38492013932672.0 / 718912673015625.0   - B_07;

  117.     /** Internal weights for stage 15, element 8. */
  118.     private static final double K15_08 = -1114881286517557.0 / 20298710767500000.0 - B_08;

  120.     /** Internal weights for stage 15, element 9. */
  121.     private static final double K15_09 =                 0.0                       - B_09;

  123.     /** Internal weights for stage 15, element 10. */
  124.     private static final double K15_10 =                 0.0                       - B_10;

  126.     /** Internal weights for stage 15, element 11. */
  127.     private static final double K15_11 =    -2538710946863.0 / 23431227861250000.0 - B_11;

  129.     /** Internal weights for stage 15, element 12. */
  130.     private static final double K15_12 =        8824659001.0 / 23066716781250.0    - B_12;

  132.     /** Internal weights for stage 15, element 13. */
  133.     private static final double K15_13 =      -11518334563.0 / 33831184612500.0;

  135.     /** Internal weights for stage 15, element 14. */
  136.     private static final double K15_14 =        1912306948.0 / 13532473845.0;

  138.     /** Time step for stage 16 (interpolation only). */
  139.     private static final double C16    = 7.0 / 9.0;


  142.     /** Internal weights for stage 16, element 1. */
  143.     private static final double K16_01 =      -13613986967.0 / 31741908048.0       - B_01;

  145.     // elements 2 to 5 are zero, so they are neither stored nor used

  147.     /** Internal weights for stage 16, element 6. */
  148.     private static final double K16_06 =       -4755612631.0 / 1012344804.0        - B_06;

  150.     /** Internal weights for stage 16, element 7. */
  151.     private static final double K16_07 =    42939257944576.0 / 5588559685701.0     - B_07;

  153.     /** Internal weights for stage 16, element 8. */
  154.     private static final double K16_08 =    77881972900277.0 / 19140370552944.0    - B_08;

  156.     /** Internal weights for stage 16, element 9. */
  157.     private static final double K16_09 =    22719829234375.0 / 63689648654052.0    - B_09;

  159.     /** Internal weights for stage 16, element 10. */
  160.     private static final double K16_10 =                 0.0                       - B_10;

  162.     /** Internal weights for stage 16, element 11. */
  163.     private static final double K16_11 =                 0.0                       - B_11;

  165.     /** Internal weights for stage 16, element 12. */
  166.     private static final double K16_12 =                 0.0                       - B_12;

  168.     /** Internal weights for stage 16, element 13. */
  169.     private static final double K16_13 =       -1199007803.0 / 857031517296.0;

  171.     /** Internal weights for stage 16, element 14. */
  172.     private static final double K16_14 =      157882067000.0 / 53564469831.0;

  174.     /** Internal weights for stage 16, element 15. */
  175.     private static final double K16_15 =     -290468882375.0 / 31741908048.0;

  177.     /** Interpolation weights.
  178.      * (beware that only the non-null values are in the table)
  179.      */
  180.     private static final double[][] D = {

  182.       {        -17751989329.0 / 2106076560.0,               4272954039.0 / 7539864640.0,
  183.               -118476319744.0 / 38604839385.0,            755123450731.0 / 316657731600.0,
  184.         3692384461234828125.0 / 1744130441634250432.0,     -4612609375.0 / 5293382976.0,
  185.               2091772278379.0 / 933644586600.0,             2136624137.0 / 3382989120.0,
  186.                     -126493.0 / 1421424.0,                    98350000.0 / 5419179.0,
  187.                   -18878125.0 / 2053168.0,                 -1944542619.0 / 438351368.0},

  189.       {         32941697297.0 / 3159114840.0,             456696183123.0 / 1884966160.0,
  190.              19132610714624.0 / 115814518155.0,       -177904688592943.0 / 474986597400.0,
  191.        -4821139941836765625.0 / 218016305204281304.0,      30702015625.0 / 3970037232.0,
  192.             -85916079474274.0 / 2800933759800.0,           -5919468007.0 / 634310460.0,
  193.                     2479159.0 / 157936.0,                    -18750000.0 / 602131.0,
  194.                   -19203125.0 / 2053168.0,                 15700361463.0 / 438351368.0},

  196.       {         12627015655.0 / 631822968.0,              -72955222965.0 / 188496616.0,
  197.             -13145744952320.0 / 69488710893.0,          30084216194513.0 / 56998391688.0,
  198.         -296858761006640625.0 / 25648977082856624.0,         569140625.0 / 82709109.0,
  199.                -18684190637.0 / 18672891732.0,                69644045.0 / 89549712.0,
  200.                   -11847025.0 / 4264272.0,                  -978650000.0 / 16257537.0,
  201.                   519371875.0 / 6159504.0,                  5256837225.0 / 438351368.0},

  203.       {          -450944925.0 / 17550638.0,               -14532122925.0 / 94248308.0,
  204.               -595876966400.0 / 2573655959.0,             188748653015.0 / 527762886.0,
  205.         2545485458115234375.0 / 27252038150535163.0,       -1376953125.0 / 36759604.0,
  206.                 53995596795.0 / 518691437.0,                 210311225.0 / 7047894.0,
  207.                    -1718875.0 / 39484.0,                      58000000.0 / 602131.0,
  208.                    -1546875.0 / 39484.0,                   -1262172375.0 / 8429834.0}
  209.     };

  211.     /** Last evaluations. */
  212.     private double[][] yDotKLast;

  214.     /** Vectors for interpolation. */
  215.     private double[][] v;

  217.     /** Initialization indicator for the interpolation vectors. */
  218.     private boolean vectorsInitialized;

  220.   /** Simple constructor.
  221.    * This constructor builds an instance that is not usable yet, the
  222.    * {@link #reinitialize} method should be called before using the
  223.    * instance in order to initialize the internal arrays. This
  224.    * constructor is used only in order to delay the initialization in
  225.    * some cases. The {@link EmbeddedRungeKuttaIntegrator} uses the
  226.    * prototyping design pattern to create the step interpolators by
  227.    * cloning an uninitialized model and latter initializing the copy.
  228.    */
  229.   // CHECKSTYLE: stop RedundantModifier
  230.   // the public modifier here is needed for serialization
  231.   public DormandPrince853StepInterpolator() {
  232.     super();
  233.     yDotKLast = null;
  234.     v         = null;
  235.     vectorsInitialized = false;
  236.   }
  237.   // CHECKSTYLE: resume RedundantModifier

  239.   /** Copy constructor.
  240.    * @param interpolator interpolator to copy from. The copy is a deep
  241.    * copy: its arrays are separated from the original arrays of the
  242.    * instance
  243.    */
  244.   DormandPrince853StepInterpolator(final DormandPrince853StepInterpolator interpolator) {

  246.     super(interpolator);

  248.     if (interpolator.currentState == null) {

  250.       yDotKLast = null;
  251.       v         = null;
  252.       vectorsInitialized = false;
  253.     } else {

  255.       final int dimension = interpolator.currentState.length;

  257.       yDotKLast    = new double[3][];
  258.       for (int k = 0; k < yDotKLast.length; ++k) {
  259.         yDotKLast[k] = new double[dimension];
  260.         System.arraycopy(interpolator.yDotKLast[k], 0, yDotKLast[k], 0,
  261.                          dimension);
  262.       }

  264.       v = new double[7][];
  265.       for (int k = 0; k < v.length; ++k) {
  266.         v[k] = new double[dimension];
  267.         System.arraycopy(interpolator.v[k], 0, v[k], 0, dimension);
  268.       }

  270.       vectorsInitialized = interpolator.vectorsInitialized;
  271.     }
  272.   }

  274.   /** {@inheritDoc} */
  275.   @Override
  276.   protected StepInterpolator doCopy() {
  277.     return new DormandPrince853StepInterpolator(this);
  278.   }

  280.   /** {@inheritDoc} */
  281.   @Override
  282.   public void reinitialize(final AbstractIntegrator integrator,
  283.                            final double[] y, final double[][] yDotK, final boolean forward,
  284.                            final EquationsMapper primaryMapper,
  285.                            final EquationsMapper[] secondaryMappers) {

  287.     super.reinitialize(integrator, y, yDotK, forward, primaryMapper, secondaryMappers);

  289.     final int dimension = currentState.length;

  291.     yDotKLast = new double[3][];
  292.     for (int k = 0; k < yDotKLast.length; ++k) {
  293.       yDotKLast[k] = new double[dimension];
  294.     }

  296.     v = new double[7][];
  297.     for (int k = 0; k < v.length; ++k) {
  298.       v[k]  = new double[dimension];
  299.     }

  301.     vectorsInitialized = false;
  302.   }

  304.   /** {@inheritDoc} */
  305.   @Override
  306.   public void storeTime(final double t) {
  307.     super.storeTime(t);
  308.     vectorsInitialized = false;
  309.   }

  311.   /** {@inheritDoc} */
  312.   @Override
  313.   protected void computeInterpolatedStateAndDerivatives(final double theta,
  314.                                           final double oneMinusThetaH)
  315.       throws MaxCountExceededException {

  317.     if (! vectorsInitialized) {

  319.       if (v == null) {
  320.         v = new double[7][];
  321.         for (int k = 0; k < 7; ++k) {
  322.           v[k] = new double[interpolatedState.length];
  323.         }
  324.       }

  326.       // perform the last evaluations if they have not been done yet
  327.       finalizeStep();

  329.       // compute the interpolation vectors for this time step
  330.       for (int i = 0; i < interpolatedState.length; ++i) {
  331.           final double yDot1  = yDotK[0][i];
  332.           final double yDot6  = yDotK[5][i];
  333.           final double yDot7  = yDotK[6][i];
  334.           final double yDot8  = yDotK[7][i];
  335.           final double yDot9  = yDotK[8][i];
  336.           final double yDot10 = yDotK[9][i];
  337.           final double yDot11 = yDotK[10][i];
  338.           final double yDot12 = yDotK[11][i];
  339.           final double yDot13 = yDotK[12][i];
  340.           final double yDot14 = yDotKLast[0][i];
  341.           final double yDot15 = yDotKLast[1][i];
  342.           final double yDot16 = yDotKLast[2][i];
  343.           v[0][i] = B_01 * yDot1  + B_06 * yDot6 + B_07 * yDot7 +
  344.                     B_08 * yDot8  + B_09 * yDot9 + B_10 * yDot10 +
  345.                     B_11 * yDot11 + B_12 * yDot12;
  346.           v[1][i] = yDot1 - v[0][i];
  347.           v[2][i] = v[0][i] - v[1][i] - yDotK[12][i];
  348.           for (int k = 0; k < D.length; ++k) {
  349.               v[k+3][i] = D[k][0] * yDot1  + D[k][1]  * yDot6  + D[k][2]  * yDot7  +
  350.                           D[k][3] * yDot8  + D[k][4]  * yDot9  + D[k][5]  * yDot10 +
  351.                           D[k][6] * yDot11 + D[k][7]  * yDot12 + D[k][8]  * yDot13 +
  352.                           D[k][9] * yDot14 + D[k][10] * yDot15 + D[k][11] * yDot16;
  353.           }
  354.       }

  356.       vectorsInitialized = true;
  357.     }

  359.     final double eta      = 1 - theta;
  360.     final double twoTheta = 2 * theta;
  361.     final double theta2   = theta * theta;
  362.     final double dot1 = 1 - twoTheta;
  363.     final double dot2 = theta * (2 - 3 * theta);
  364.     final double dot3 = twoTheta * (1 + theta * (twoTheta -3));
  365.     final double dot4 = theta2 * (3 + theta * (5 * theta - 8));
  366.     final double dot5 = theta2 * (3 + theta * (-12 + theta * (15 - 6 * theta)));
  367.     final double dot6 = theta2 * theta * (4 + theta * (-15 + theta * (18 - 7 * theta)));

  369.     if (previousState != null && theta <= 0.5) {
  370.         for (int i = 0; i < interpolatedState.length; ++i) {
  371.             interpolatedState[i] = previousState[i] +
  372.                     theta * h * (v[0][i] +
  373.                             eta * (v[1][i] +
  374.                                     theta * (v[2][i] +
  375.                                             eta * (v[3][i] +
  376.                                                     theta * (v[4][i] +
  377.                                                             eta * (v[5][i] +
  378.                                                                     theta * (v[6][i])))))));
  379.             interpolatedDerivatives[i] =  v[0][i] + dot1 * v[1][i] + dot2 * v[2][i] +
  380.                     dot3 * v[3][i] + dot4 * v[4][i] +
  381.                     dot5 * v[5][i] + dot6 * v[6][i];
  382.         }
  383.     } else {
  384.         for (int i = 0; i < interpolatedState.length; ++i) {
  385.             interpolatedState[i] = currentState[i] -
  386.                     oneMinusThetaH * (v[0][i] -
  387.                             theta * (v[1][i] +
  388.                                     theta * (v[2][i] +
  389.                                             eta * (v[3][i] +
  390.                                                     theta * (v[4][i] +
  391.                                                             eta * (v[5][i] +
  392.                                                                     theta * (v[6][i])))))));
  393.             interpolatedDerivatives[i] =  v[0][i] + dot1 * v[1][i] + dot2 * v[2][i] +
  394.                     dot3 * v[3][i] + dot4 * v[4][i] +
  395.                     dot5 * v[5][i] + dot6 * v[6][i];
  396.         }
  397.     }
  398.   }

  400.   /** {@inheritDoc} */
  401.   @Override
  402.   protected void doFinalize() throws MaxCountExceededException {

  404.       if (currentState == null) {
  405.           // we are finalizing an uninitialized instance
  406.           return;
  407.       }

  409.       double s;
  410.       final double[] yTmp = new double[currentState.length];
  411.       final double pT = getGlobalPreviousTime();

  413.       // k14
  414.       for (int j = 0; j < currentState.length; ++j) {
  415.           s = K14_01 * yDotK[0][j]  + K14_06 * yDotK[5][j]  + K14_07 * yDotK[6][j] +
  416.                   K14_08 * yDotK[7][j]  + K14_09 * yDotK[8][j]  + K14_10 * yDotK[9][j] +
  417.                   K14_11 * yDotK[10][j] + K14_12 * yDotK[11][j] + K14_13 * yDotK[12][j];
  418.           yTmp[j] = currentState[j] + h * s;
  419.       }
  420.       integrator.computeDerivatives(pT + C14 * h, yTmp, yDotKLast[0]);

  422.       // k15
  423.       for (int j = 0; j < currentState.length; ++j) {
  424.           s = K15_01 * yDotK[0][j]  + K15_06 * yDotK[5][j]  + K15_07 * yDotK[6][j] +
  425.                   K15_08 * yDotK[7][j]  + K15_09 * yDotK[8][j]  + K15_10 * yDotK[9][j] +
  426.                   K15_11 * yDotK[10][j] + K15_12 * yDotK[11][j] + K15_13 * yDotK[12][j] +
  427.                   K15_14 * yDotKLast[0][j];
  428.           yTmp[j] = currentState[j] + h * s;
  429.       }
  430.       integrator.computeDerivatives(pT + C15 * h, yTmp, yDotKLast[1]);

  432.       // k16
  433.       for (int j = 0; j < currentState.length; ++j) {
  434.           s = K16_01 * yDotK[0][j]  + K16_06 * yDotK[5][j]  + K16_07 * yDotK[6][j] +
  435.                   K16_08 * yDotK[7][j]  + K16_09 * yDotK[8][j]  + K16_10 * yDotK[9][j] +
  436.                   K16_11 * yDotK[10][j] + K16_12 * yDotK[11][j] + K16_13 * yDotK[12][j] +
  437.                   K16_14 * yDotKLast[0][j] +  K16_15 * yDotKLast[1][j];
  438.           yTmp[j] = currentState[j] + h * s;
  439.       }
  440.       integrator.computeDerivatives(pT + C16 * h, yTmp, yDotKLast[2]);
  441.   }

  443.   /** {@inheritDoc} */
  444.   @Override
  445.   public void writeExternal(final ObjectOutput out)
  446.     throws IOException {

  448.     try {
  449.         // save the local attributes
  450.         finalizeStep();
  451.     } catch (MaxCountExceededException mcee) {
  452.         final IOException ioe = new IOException(mcee.getLocalizedMessage());
  453.         ioe.initCause(mcee);
  454.         throw ioe;
  455.     }

  457.     final int dimension = (currentState == null) ? -1 : currentState.length;
  458.     out.writeInt(dimension);
  459.     for (int i = 0; i < dimension; ++i) {
  460.       out.writeDouble(yDotKLast[0][i]);
  461.       out.writeDouble(yDotKLast[1][i]);
  462.       out.writeDouble(yDotKLast[2][i]);
  463.     }

  465.     // save the state of the base class
  466.     super.writeExternal(out);
  467.   }

  469.   /** {@inheritDoc} */
  470.   @Override
  471.   public void readExternal(final ObjectInput in)
  472.     throws IOException, ClassNotFoundException {

  474.     // read the local attributes
  475.     yDotKLast = new double[3][];
  476.     final int dimension = in.readInt();
  477.     yDotKLast[0] = (dimension < 0) ? null : new double[dimension];
  478.     yDotKLast[1] = (dimension < 0) ? null : new double[dimension];
  479.     yDotKLast[2] = (dimension < 0) ? null : new double[dimension];

  481.     for (int i = 0; i < dimension; ++i) {
  482.       yDotKLast[0][i] = in.readDouble();
  483.       yDotKLast[1][i] = in.readDouble();
  484.       yDotKLast[2][i] = in.readDouble();
  485.     }

  487.     // read the base state
  488.     super.readExternal(in);
  489.   }
  490. }
Created with JaCoCo 0.8.10.202304240956Code Coverage Report for Apache Commons Math 4.0-SNAPSHOT