AdamsNordsieckFieldTransformer.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.  *      https://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. /*
  19.  * This is not the original file distributed by the Apache Software Foundation
  20.  * It has been modified by the Hipparchus project
  21.  */

  23. package org.hipparchus.ode.nonstiff;

  25. import java.util.Arrays;
  26. import java.util.HashMap;
  27. import java.util.Map;

  29. import org.hipparchus.Field;
  30. import org.hipparchus.CalculusFieldElement;
  31. import org.hipparchus.linear.Array2DRowFieldMatrix;
  32. import org.hipparchus.linear.ArrayFieldVector;
  33. import org.hipparchus.linear.FieldDecompositionSolver;
  34. import org.hipparchus.linear.FieldLUDecomposition;
  35. import org.hipparchus.linear.FieldMatrix;
  36. import org.hipparchus.util.MathArrays;

  38. /** Transformer to Nordsieck vectors for Adams integrators.
  39.  * <p>This class is used by {@link AdamsBashforthIntegrator Adams-Bashforth} and
  40.  * {@link AdamsMoultonIntegrator Adams-Moulton} integrators to convert between
  41.  * classical representation with several previous first derivatives and Nordsieck
  42.  * representation with higher order scaled derivatives.</p>
  43.  *
  44.  * <p>We define scaled derivatives s<sub>i</sub>(n) at step n as:
  45.  * \[
  46.  *   \left\{\begin{align}
  47.  *   s_1(n) &amp;= h y'_n \text{ for first derivative}\\
  48.  *   s_2(n) &amp;= \frac{h^2}{2} y_n'' \text{ for second derivative}\\
  49.  *   s_3(n) &amp;= \frac{h^3}{6} y_n''' \text{ for third derivative}\\
  50.  *   &amp;\cdots\\
  51.  *   s_k(n) &amp;= \frac{h^k}{k!} y_n^{(k)} \text{ for } k^\mathrm{th} \text{ derivative}
  52.  *   \end{align}\right.
  53.  * \]</p>
  54.  *
  55.  * <p>With the previous definition, the classical representation of multistep methods
  56.  * uses first derivatives only, i.e. it handles y<sub>n</sub>, s<sub>1</sub>(n) and
  57.  * q<sub>n</sub> where q<sub>n</sub> is defined as:
  58.  * \[
  59.  *   q_n = [ s_1(n-1) s_1(n-2) \ldots s_1(n-(k-1)) ]^T
  60.  * \]
  61.  * (we omit the k index in the notation for clarity).</p>
  62.  *
  63.  * <p>Another possible representation uses the Nordsieck vector with
  64.  * higher degrees scaled derivatives all taken at the same step, i.e it handles y<sub>n</sub>,
  65.  * s<sub>1</sub>(n) and r<sub>n</sub>) where r<sub>n</sub> is defined as:
  66.  * \[
  67.  * r_n = [ s_2(n), s_3(n) \ldots s_k(n) ]^T
  68.  * \]
  69.  * (here again we omit the k index in the notation for clarity)
  70.  * </p>
  71.  *
  72.  * <p>Taylor series formulas show that for any index offset i, s<sub>1</sub>(n-i) can be
  73.  * computed from s<sub>1</sub>(n), s<sub>2</sub>(n) ... s<sub>k</sub>(n), the formula being exact
  74.  * for degree k polynomials.
  75.  * \[
  76.  * s_1(n-i) = s_1(n) + \sum_{j\gt 0} (j+1) (-i)^j s_{j+1}(n)
  77.  * \]
  78.  * The previous formula can be used with several values for i to compute the transform between
  79.  * classical representation and Nordsieck vector at step end. The transform between r<sub>n</sub>
  80.  * and q<sub>n</sub> resulting from the Taylor series formulas above is:
  81.  * \[
  82.  * q_n = s_1(n) u + P r_n
  83.  * \]
  84.  * where u is the [ 1 1 ... 1 ]<sup>T</sup> vector and P is the (k-1)&times;(k-1) matrix built
  85.  * with the \((j+1) (-i)^j\) terms with i being the row number starting from 1 and j being
  86.  * the column number starting from 1:
  87.  * \[
  88.  *   P=\begin{bmatrix}
  89.  *   -2  &amp;  3 &amp;   -4 &amp;    5 &amp; \ldots \\
  90.  *   -4  &amp; 12 &amp;  -32 &amp;   80 &amp; \ldots \\
  91.  *   -6  &amp; 27 &amp; -108 &amp;  405 &amp; \ldots \\
  92.  *   -8  &amp; 48 &amp; -256 &amp; 1280 &amp; \ldots \\
  93.  *       &amp;    &amp;  \ldots\\
  94.  *    \end{bmatrix}
  95.  * \]
  96.  *
  97.  * <p>Changing -i into +i in the formula above can be used to compute a similar transform between
  98.  * classical representation and Nordsieck vector at step start. The resulting matrix is simply
  99.  * the absolute value of matrix P.</p>
  100.  *
  101.  * <p>For {@link AdamsBashforthIntegrator Adams-Bashforth} method, the Nordsieck vector
  102.  * at step n+1 is computed from the Nordsieck vector at step n as follows:
  103.  * <ul>
  104.  *   <li>y<sub>n+1</sub> = y<sub>n</sub> + s<sub>1</sub>(n) + u<sup>T</sup> r<sub>n</sub></li>
  105.  *   <li>s<sub>1</sub>(n+1) = h f(t<sub>n+1</sub>, y<sub>n+1</sub>)</li>
  106.  *   <li>r<sub>n+1</sub> = (s<sub>1</sub>(n) - s<sub>1</sub>(n+1)) P<sup>-1</sup> u + P<sup>-1</sup> A P r<sub>n</sub></li>
  107.  * </ul>
  108.  * <p>where A is a rows shifting matrix (the lower left part is an identity matrix):</p>
  109.  * <pre>
  110.  *        [ 0 0   ...  0 0 | 0 ]
  111.  *        [ ---------------+---]
  112.  *        [ 1 0   ...  0 0 | 0 ]
  113.  *    A = [ 0 1   ...  0 0 | 0 ]
  114.  *        [       ...      | 0 ]
  115.  *        [ 0 0   ...  1 0 | 0 ]
  116.  *        [ 0 0   ...  0 1 | 0 ]
  117.  * </pre>
  118.  *
  119.  * <p>For {@link AdamsMoultonIntegrator Adams-Moulton} method, the predicted Nordsieck vector
  120.  * at step n+1 is computed from the Nordsieck vector at step n as follows:
  121.  * <ul>
  122.  *   <li>Y<sub>n+1</sub> = y<sub>n</sub> + s<sub>1</sub>(n) + u<sup>T</sup> r<sub>n</sub></li>
  123.  *   <li>S<sub>1</sub>(n+1) = h f(t<sub>n+1</sub>, Y<sub>n+1</sub>)</li>
  124.  *   <li>R<sub>n+1</sub> = (s<sub>1</sub>(n) - s<sub>1</sub>(n+1)) P<sup>-1</sup> u + P<sup>-1</sup> A P r<sub>n</sub></li>
  125.  * </ul>
  126.  * From this predicted vector, the corrected vector is computed as follows:
  127.  * <ul>
  128.  *   <li>y<sub>n+1</sub> = y<sub>n</sub> + S<sub>1</sub>(n+1) + [ -1 +1 -1 +1 ... &plusmn;1 ] r<sub>n+1</sub></li>
  129.  *   <li>s<sub>1</sub>(n+1) = h f(t<sub>n+1</sub>, y<sub>n+1</sub>)</li>
  130.  *   <li>r<sub>n+1</sub> = R<sub>n+1</sub> + (s<sub>1</sub>(n+1) - S<sub>1</sub>(n+1)) P<sup>-1</sup> u</li>
  131.  * </ul>
  132.  * <p>where the upper case Y<sub>n+1</sub>, S<sub>1</sub>(n+1) and R<sub>n+1</sub> represent the
  133.  * predicted states whereas the lower case y<sub>n+1</sub>, s<sub>n+1</sub> and r<sub>n+1</sub>
  134.  * represent the corrected states.</p>
  135.  *
  136.  * <p>We observe that both methods use similar update formulas. In both cases a P<sup>-1</sup>u
  137.  * vector and a P<sup>-1</sup> A P matrix are used that do not depend on the state,
  138.  * they only depend on k. This class handles these transformations.</p>
  139.  *
  140.  * @param <T> the type of the field elements
  141.  */
  142. public class AdamsNordsieckFieldTransformer<T extends CalculusFieldElement<T>> {

  144.     /** Cache for already computed coefficients. */
  145.     private static final Map<Integer,
  146.                          Map<Field<? extends CalculusFieldElement<?>>,
  147.                                    AdamsNordsieckFieldTransformer<? extends CalculusFieldElement<?>>>> CACHE = new HashMap<>();

  149.     /** Field to which the time and state vector elements belong. */
  150.     private final Field<T> field;

  152.     /** Update matrix for the higher order derivatives h<sup>2</sup>/2 y'', h<sup>3</sup>/6 y''' ... */
  153.     private final Array2DRowFieldMatrix<T> update;

  155.     /** Update coefficients of the higher order derivatives wrt y'. */
  156.     private final T[] c1;

  158.     /** Simple constructor.
  159.      * @param field field to which the time and state vector elements belong
  160.      * @param n number of steps of the multistep method
  161.      * (excluding the one being computed)
  162.      */
  163.     private AdamsNordsieckFieldTransformer(final Field<T> field, final int n) {

  165.         this.field = field;
  166.         final int rows = n - 1;

  168.         // compute coefficients
  169.         FieldMatrix<T> bigP = buildP(rows);
  170.         FieldDecompositionSolver<T> pSolver =
  171.                 new FieldLUDecomposition<>(bigP).getSolver();

  173.         T[] u = MathArrays.buildArray(field, rows);
  174.         Arrays.fill(u, field.getOne());
  175.         c1 = pSolver.solve(new ArrayFieldVector<>(u, false)).toArray();

  177.         // update coefficients are computed by combining transform from
  178.         // Nordsieck to multistep, then shifting rows to represent step advance
  179.         // then applying inverse transform
  180.         T[][] shiftedP = bigP.getData();
  181.         for (int i = shiftedP.length - 1; i > 0; --i) {
  182.             // shift rows
  183.             shiftedP[i] = shiftedP[i - 1];
  184.         }
  185.         shiftedP[0] = MathArrays.buildArray(field, rows);
  186.         Arrays.fill(shiftedP[0], field.getZero());
  187.         update = new Array2DRowFieldMatrix<>(pSolver.solve(new Array2DRowFieldMatrix<>(shiftedP, false)).getData());

  189.     }

  191.     /** Get the Nordsieck transformer for a given field and number of steps.
  192.      * @param field field to which the time and state vector elements belong
  193.      * @param nSteps number of steps of the multistep method
  194.      * (excluding the one being computed)
  195.      * @return Nordsieck transformer for the specified field and number of steps
  196.      * @param <T> the type of the field elements
  197.      */
  198.     public static <T extends CalculusFieldElement<T>> AdamsNordsieckFieldTransformer<T>
  199.     getInstance(final Field<T> field, final int nSteps) { // NOPMD - PMD false positive
  200.         synchronized(CACHE) {
  201.             Map<Field<? extends CalculusFieldElement<?>>,
  202.                     AdamsNordsieckFieldTransformer<? extends CalculusFieldElement<?>>> map = CACHE.computeIfAbsent(nSteps, k -> new HashMap<>());
  203.             @SuppressWarnings("unchecked")
  204.             AdamsNordsieckFieldTransformer<T> t = (AdamsNordsieckFieldTransformer<T>) map.get(field);
  205.             if (t == null) {
  206.                 t = new AdamsNordsieckFieldTransformer<>(field, nSteps);
  207.                 map.put(field, t);
  208.             }
  209.             return t;

  211.         }
  212.     }

  214.     /** Build the P matrix.
  215.      * <p>The P matrix general terms are shifted \((j+1) (-i)^j\) terms
  216.      * with i being the row number starting from 1 and j being the column
  217.      * number starting from 1:
  218.      * <pre>
  219.      *        [  -2   3   -4    5  ... ]
  220.      *        [  -4  12  -32   80  ... ]
  221.      *   P =  [  -6  27 -108  405  ... ]
  222.      *        [  -8  48 -256 1280  ... ]
  223.      *        [          ...           ]
  224.      * </pre></p>
  225.      * @param rows number of rows of the matrix
  226.      * @return P matrix
  227.      */
  228.     private FieldMatrix<T> buildP(final int rows) {

  230.         final T[][] pData = MathArrays.buildArray(field, rows, rows);

  232.         for (int i = 1; i <= pData.length; ++i) {
  233.             // build the P matrix elements from Taylor series formulas
  234.             final T[] pI = pData[i - 1];
  235.             final int factor = -i;
  236.             T aj = field.getZero().add(factor);
  237.             for (int j = 1; j <= pI.length; ++j) {
  238.                 pI[j - 1] = aj.multiply(j + 1);
  239.                 aj = aj.multiply(factor);
  240.             }
  241.         }

  243.         return new Array2DRowFieldMatrix<>(pData, false);

  245.     }

  247.     /** Initialize the high order scaled derivatives at step start.
  248.      * @param h step size to use for scaling
  249.      * @param t first steps times
  250.      * @param y first steps states
  251.      * @param yDot first steps derivatives
  252.      * @return Nordieck vector at start of first step (h<sup>2</sup>/2 y''<sub>n</sub>,
  253.      * h<sup>3</sup>/6 y'''<sub>n</sub> ... h<sup>k</sup>/k! y<sup>(k)</sup><sub>n</sub>)
  254.      */

  256.     public Array2DRowFieldMatrix<T> initializeHighOrderDerivatives(final T h, final T[] t,
  257.                                                                    final T[][] y,
  258.                                                                    final T[][] yDot) {

  260.         // using Taylor series with di = ti - t0, we get:
  261.         //  y(ti)  - y(t0)  - di y'(t0) =   di^2 / h^2 s2 + ... +   di^k     / h^k sk + O(h^k)
  262.         //  y'(ti) - y'(t0)             = 2 di   / h^2 s2 + ... + k di^(k-1) / h^k sk + O(h^(k-1))
  263.         // we write these relations for i = 1 to i= 1+n/2 as a set of n + 2 linear
  264.         // equations depending on the Nordsieck vector [s2 ... sk rk], so s2 to sk correspond
  265.         // to the appropriately truncated Taylor expansion, and rk is the Taylor remainder.
  266.         // The goal is to have s2 to sk as accurate as possible considering the fact the sum is
  267.         // truncated and we don't want the error terms to be included in s2 ... sk, so we need
  268.         // to solve also for the remainder
  269.         final T[][] a     = MathArrays.buildArray(field, c1.length + 1, c1.length + 1);
  270.         final T[][] b     = MathArrays.buildArray(field, c1.length + 1, y[0].length);
  271.         final T[]   y0    = y[0];
  272.         final T[]   yDot0 = yDot[0];
  273.         for (int i = 1; i < y.length; ++i) {

  275.             final T di    = t[i].subtract(t[0]);
  276.             final T ratio = di.divide(h);
  277.             T dikM1Ohk    = h.reciprocal();

  279.             // linear coefficients of equations
  280.             // y(ti) - y(t0) - di y'(t0) and y'(ti) - y'(t0)
  281.             final T[] aI    = a[2 * i - 2];
  282.             final T[] aDotI = (2 * i - 1) < a.length ? a[2 * i - 1] : null;
  283.             for (int j = 0; j < aI.length; ++j) {
  284.                 dikM1Ohk = dikM1Ohk.multiply(ratio);
  285.                 aI[j]    = di.multiply(dikM1Ohk);
  286.                 if (aDotI != null) {
  287.                     aDotI[j]  = dikM1Ohk.multiply(j + 2);
  288.                 }
  289.             }

  291.             // expected value of the previous equations
  292.             final T[] yI    = y[i];
  293.             final T[] yDotI = yDot[i];
  294.             final T[] bI    = b[2 * i - 2];
  295.             final T[] bDotI = (2 * i - 1) < b.length ? b[2 * i - 1] : null;
  296.             for (int j = 0; j < yI.length; ++j) {
  297.                 bI[j]    = yI[j].subtract(y0[j]).subtract(di.multiply(yDot0[j]));
  298.                 if (bDotI != null) {
  299.                     bDotI[j] = yDotI[j].subtract(yDot0[j]);
  300.                 }
  301.             }

  303.         }

  305.         // solve the linear system to get the best estimate of the Nordsieck vector [s2 ... sk],
  306.         // with the additional terms s(k+1) and c grabbing the parts after the truncated Taylor expansion
  307.         final FieldLUDecomposition<T> decomposition = new FieldLUDecomposition<>(new Array2DRowFieldMatrix<>(a, false));
  308.         final FieldMatrix<T> x = decomposition.getSolver().solve(new Array2DRowFieldMatrix<>(b, false));

  310.         // extract just the Nordsieck vector [s2 ... sk]
  311.         final Array2DRowFieldMatrix<T> truncatedX =
  312.                         new Array2DRowFieldMatrix<>(field, x.getRowDimension() - 1, x.getColumnDimension());
  313.         for (int i = 0; i < truncatedX.getRowDimension(); ++i) {
  314.             for (int j = 0; j < truncatedX.getColumnDimension(); ++j) {
  315.                 truncatedX.setEntry(i, j, x.getEntry(i, j));
  316.             }
  317.         }
  318.         return truncatedX;

  320.     }

  322.     /** Update the high order scaled derivatives for Adams integrators (phase 1).
  323.      * <p>The complete update of high order derivatives has a form similar to:
  324.      * \[
  325.      * r_{n+1} = (s_1(n) - s_1(n+1)) P^{-1} u + P^{-1} A P r_n
  326.      * \]
  327.      * this method computes the P<sup>-1</sup> A P r<sub>n</sub> part.</p>
  328.      * @param highOrder high order scaled derivatives
  329.      * (h<sup>2</sup>/2 y'', ... h<sup>k</sup>/k! y(k))
  330.      * @return updated high order derivatives
  331.      * @see #updateHighOrderDerivativesPhase2(CalculusFieldElement[], CalculusFieldElement[], Array2DRowFieldMatrix)
  332.      */
  333.     public Array2DRowFieldMatrix<T> updateHighOrderDerivativesPhase1(final Array2DRowFieldMatrix<T> highOrder) {
  334.         return update.multiply(highOrder);
  335.     }

  337.     /** Update the high order scaled derivatives Adams integrators (phase 2).
  338.      * <p>The complete update of high order derivatives has a form similar to:
  339.      * \[
  340.      * r_{n+1} = (s_1(n) - s_1(n+1)) P^{-1} u + P^{-1} A P r_n
  341.      * \]
  342.      * this method computes the (s<sub>1</sub>(n) - s<sub>1</sub>(n+1)) P<sup>-1</sup> u part.</p>
  343.      * <p>Phase 1 of the update must already have been performed.</p>
  344.      * @param start first order scaled derivatives at step start
  345.      * @param end first order scaled derivatives at step end
  346.      * @param highOrder high order scaled derivatives, will be modified
  347.      * (h<sup>2</sup>/2 y'', ... h<sup>k</sup>/k! y(k))
  348.      * @see #updateHighOrderDerivativesPhase1(Array2DRowFieldMatrix)
  349.      */
  350.     public void updateHighOrderDerivativesPhase2(final T[] start,
  351.                                                  final T[] end,
  352.                                                  final Array2DRowFieldMatrix<T> highOrder) {
  353.         final T[][] data = highOrder.getDataRef();
  354.         for (int i = 0; i < data.length; ++i) {
  355.             final T[] dataI = data[i];
  356.             final T c1I = c1[i];
  357.             for (int j = 0; j < dataI.length; ++j) {
  358.                 dataI[j] = dataI[j].add(c1I.multiply(start[j].subtract(end[j])));
  359.             }
  360.         }
  361.     }

  363. }
Created with JaCoCo 0.8.12.202403310830