View Javadoc
1   /*
2    * Licensed to the Hipparchus project 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 Hipparchus project 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   */
17  
18  package org.hipparchus.ode.nonstiff;
19  
20  import java.util.Arrays;
21  import java.util.HashMap;
22  import java.util.Map;
23  
24  import org.hipparchus.fraction.BigFraction;
25  import org.hipparchus.linear.Array2DRowFieldMatrix;
26  import org.hipparchus.linear.Array2DRowRealMatrix;
27  import org.hipparchus.linear.ArrayFieldVector;
28  import org.hipparchus.linear.FieldDecompositionSolver;
29  import org.hipparchus.linear.FieldLUDecomposition;
30  import org.hipparchus.linear.FieldMatrix;
31  import org.hipparchus.linear.MatrixUtils;
32  import org.hipparchus.linear.QRDecomposition;
33  import org.hipparchus.linear.RealMatrix;
34  
35  /** Transformer to Nordsieck vectors for Adams integrators.
36   * <p>This class is used by {@link AdamsBashforthIntegrator Adams-Bashforth} and
37   * {@link AdamsMoultonIntegrator Adams-Moulton} integrators to convert between
38   * classical representation with several previous first derivatives and Nordsieck
39   * representation with higher order scaled derivatives.</p>
40   *
41   * <p>We define scaled derivatives s<sub>i</sub>(n) at step n as:
42   * \[
43   *   \left\{\begin{align}
44   *   s_1(n) &amp;= h y'_n \text{ for first derivative}\\
45   *   s_2(n) &amp;= \frac{h^2}{2} y_n'' \text{ for second derivative}\\
46   *   s_3(n) &amp;= \frac{h^3}{6} y_n''' \text{ for third derivative}\\
47   *   &amp;\cdots\\
48   *   s_k(n) &amp;= \frac{h^k}{k!} y_n^{(k)} \text{ for } k^\mathrm{th} \text{ derivative}
49   *   \end{align}\right.
50   * \]</p>
51   *
52   * <p>With the previous definition, the classical representation of multistep methods
53   * uses first derivatives only, i.e. it handles y<sub>n</sub>, s<sub>1</sub>(n) and
54   * q<sub>n</sub> where q<sub>n</sub> is defined as:
55   * \[
56   *   q_n = [ s_1(n-1) s_1(n-2) \ldots s_1(n-(k-1)) ]^T
57   * \]
58   * (we omit the k index in the notation for clarity).</p>
59   *
60   * <p>Another possible representation uses the Nordsieck vector with
61   * higher degrees scaled derivatives all taken at the same step, i.e it handles y<sub>n</sub>,
62   * s<sub>1</sub>(n) and r<sub>n</sub>) where r<sub>n</sub> is defined as:
63   * \[
64   * r_n = [ s_2(n), s_3(n) \ldots s_k(n) ]^T
65   * \]
66   * (here again we omit the k index in the notation for clarity)
67   * </p>
68   *
69   * <p>Taylor series formulas show that for any index offset i, s<sub>1</sub>(n-i) can be
70   * computed from s<sub>1</sub>(n), s<sub>2</sub>(n) ... s<sub>k</sub>(n), the formula being exact
71   * for degree k polynomials.
72   * \[
73   * s_1(n-i) = s_1(n) + \sum_{j\gt 0} (j+1) (-i)^j s_{j+1}(n)
74   * \]
75   * The previous formula can be used with several values for i to compute the transform between
76   * classical representation and Nordsieck vector at step end. The transform between r<sub>n</sub>
77   * and q<sub>n</sub> resulting from the Taylor series formulas above is:
78   * \[
79   * q_n = s_1(n) u + P r_n
80   * \]
81   * where u is the [ 1 1 ... 1 ]<sup>T</sup> vector and P is the (k-1)&times;(k-1) matrix built
82   * with the \((j+1) (-i)^j\) terms with i being the row number starting from 1 and j being
83   * the column number starting from 1:
84   * \[
85   *   P=\begin{bmatrix}
86   *   -2  &amp;  3 &amp;   -4 &amp;    5 &amp; \ldots \\
87   *   -4  &amp; 12 &amp;  -32 &amp;   80 &amp; \ldots \\
88   *   -6  &amp; 27 &amp; -108 &amp;  405 &amp; \ldots \\
89   *   -8  &amp; 48 &amp; -256 &amp; 1280 &amp; \ldots \\
90   *       &amp;    &amp;  \ldots\\
91   *    \end{bmatrix}
92   * \]
93   *
94   * <p>Changing -i into +i in the formula above can be used to compute a similar transform between
95   * classical representation and Nordsieck vector at step start. The resulting matrix is simply
96   * the absolute value of matrix P.</p>
97   *
98   * <p>For {@link AdamsBashforthIntegrator Adams-Bashforth} method, the Nordsieck vector
99   * at step n+1 is computed from the Nordsieck vector at step n as follows:
100  * <ul>
101  *   <li>y<sub>n+1</sub> = y<sub>n</sub> + s<sub>1</sub>(n) + u<sup>T</sup> r<sub>n</sub></li>
102  *   <li>s<sub>1</sub>(n+1) = h f(t<sub>n+1</sub>, y<sub>n+1</sub>)</li>
103  *   <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>
104  * </ul>
105  * <p>where A is a rows shifting matrix (the lower left part is an identity matrix):</p>
106  * <pre>
107  *        [ 0 0   ...  0 0 | 0 ]
108  *        [ ---------------+---]
109  *        [ 1 0   ...  0 0 | 0 ]
110  *    A = [ 0 1   ...  0 0 | 0 ]
111  *        [       ...      | 0 ]
112  *        [ 0 0   ...  1 0 | 0 ]
113  *        [ 0 0   ...  0 1 | 0 ]
114  * </pre>
115  *
116  * <p>For {@link AdamsMoultonIntegrator Adams-Moulton} method, the predicted Nordsieck vector
117  * at step n+1 is computed from the Nordsieck vector at step n as follows:
118  * <ul>
119  *   <li>Y<sub>n+1</sub> = y<sub>n</sub> + s<sub>1</sub>(n) + u<sup>T</sup> r<sub>n</sub></li>
120  *   <li>S<sub>1</sub>(n+1) = h f(t<sub>n+1</sub>, Y<sub>n+1</sub>)</li>
121  *   <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>
122  * </ul>
123  * From this predicted vector, the corrected vector is computed as follows:
124  * <ul>
125  *   <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>
126  *   <li>s<sub>1</sub>(n+1) = h f(t<sub>n+1</sub>, y<sub>n+1</sub>)</li>
127  *   <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>
128  * </ul>
129  * <p>where the upper case Y<sub>n+1</sub>, S<sub>1</sub>(n+1) and R<sub>n+1</sub> represent the
130  * predicted states whereas the lower case y<sub>n+1</sub>, s<sub>n+1</sub> and r<sub>n+1</sub>
131  * represent the corrected states.</p>
132  *
133  * <p>We observe that both methods use similar update formulas. In both cases a P<sup>-1</sup>u
134  * vector and a P<sup>-1</sup> A P matrix are used that do not depend on the state,
135  * they only depend on k. This class handles these transformations.</p>
136  *
137  */
138 public class AdamsNordsieckTransformer {
139 
140     /** Cache for already computed coefficients. */
141     private static final Map<Integer, AdamsNordsieckTransformer> CACHE = new HashMap<>();
142 
143     /** Update matrix for the higher order derivatives h<sup>2</sup>/2 y'', h<sup>3</sup>/6 y''' ... */
144     private final Array2DRowRealMatrix update;
145 
146     /** Update coefficients of the higher order derivatives wrt y'. */
147     private final double[] c1;
148 
149     /** Simple constructor.
150      * @param n number of steps of the multistep method
151      * (excluding the one being computed)
152      */
153     private AdamsNordsieckTransformer(final int n) {
154 
155         final int rows = n - 1;
156 
157         // compute exact coefficients
158         FieldMatrix<BigFraction> bigP = buildP(rows);
159         FieldDecompositionSolver<BigFraction> pSolver =
160             new FieldLUDecomposition<BigFraction>(bigP).getSolver();
161 
162         BigFraction[] u = new BigFraction[rows];
163         Arrays.fill(u, BigFraction.ONE);
164         BigFraction[] bigC1 = pSolver.solve(new ArrayFieldVector<BigFraction>(u, false)).toArray();
165 
166         // update coefficients are computed by combining transform from
167         // Nordsieck to multistep, then shifting rows to represent step advance
168         // then applying inverse transform
169         BigFraction[][] shiftedP = bigP.getData();
170         for (int i = shiftedP.length - 1; i > 0; --i) {
171             // shift rows
172             shiftedP[i] = shiftedP[i - 1];
173         }
174         shiftedP[0] = new BigFraction[rows];
175         Arrays.fill(shiftedP[0], BigFraction.ZERO);
176         FieldMatrix<BigFraction> bigMSupdate =
177             pSolver.solve(new Array2DRowFieldMatrix<BigFraction>(shiftedP, false));
178 
179         // convert coefficients to double
180         update         = MatrixUtils.bigFractionMatrixToRealMatrix(bigMSupdate);
181         c1             = new double[rows];
182         for (int i = 0; i < rows; ++i) {
183             c1[i] = bigC1[i].doubleValue();
184         }
185 
186     }
187 
188     /** Get the Nordsieck transformer for a given number of steps.
189      * @param nSteps number of steps of the multistep method
190      * (excluding the one being computed)
191      * @return Nordsieck transformer for the specified number of steps
192      */
193     public static AdamsNordsieckTransformer getInstance(final int nSteps) { // NOPMD - PMD false positive
194         synchronized(CACHE) {
195             AdamsNordsieckTransformer t = CACHE.get(nSteps);
196             if (t == null) {
197                 t = new AdamsNordsieckTransformer(nSteps);
198                 CACHE.put(nSteps, t);
199             }
200             return t;
201         }
202     }
203 
204     /** Build the P matrix.
205      * <p>The P matrix general terms are shifted \((j+1) (-i)^j\) terms
206      * with i being the row number starting from 1 and j being the column
207      * number starting from 1:
208      * <pre>
209      *        [  -2   3   -4    5  ... ]
210      *        [  -4  12  -32   80  ... ]
211      *   P =  [  -6  27 -108  405  ... ]
212      *        [  -8  48 -256 1280  ... ]
213      *        [          ...           ]
214      * </pre></p>
215      * @param rows number of rows of the matrix
216      * @return P matrix
217      */
218     private FieldMatrix<BigFraction> buildP(final int rows) {
219 
220         final BigFraction[][] pData = new BigFraction[rows][rows];
221 
222         for (int i = 1; i <= pData.length; ++i) {
223             // build the P matrix elements from Taylor series formulas
224             final BigFraction[] pI = pData[i - 1];
225             final int factor = -i;
226             int aj = factor;
227             for (int j = 1; j <= pI.length; ++j) {
228                 pI[j - 1] = new BigFraction(aj * (j + 1));
229                 aj *= factor;
230             }
231         }
232 
233         return new Array2DRowFieldMatrix<BigFraction>(pData, false);
234 
235     }
236 
237     /** Initialize the high order scaled derivatives at step start.
238      * @param h step size to use for scaling
239      * @param t first steps times
240      * @param y first steps states
241      * @param yDot first steps derivatives
242      * @return Nordieck vector at start of first step (h<sup>2</sup>/2 y''<sub>n</sub>,
243      * h<sup>3</sup>/6 y'''<sub>n</sub> ... h<sup>k</sup>/k! y<sup>(k)</sup><sub>n</sub>)
244      */
245 
246     public Array2DRowRealMatrix initializeHighOrderDerivatives(final double h, final double[] t,
247                                                                final double[][] y,
248                                                                final double[][] yDot) {
249 
250         // using Taylor series with di = ti - t0, we get:
251         //  y(ti)  - y(t0)  - di y'(t0) =   di^2 / h^2 s2 + ... +   di^k     / h^k sk + O(h^k)
252         //  y'(ti) - y'(t0)             = 2 di   / h^2 s2 + ... + k di^(k-1) / h^k sk + O(h^(k-1))
253         // we write these relations for i = 1 to i= 1+n/2 as a set of n + 2 linear
254         // equations depending on the Nordsieck vector [s2 ... sk rk], so s2 to sk correspond
255         // to the appropriately truncated Taylor expansion, and rk is the Taylor remainder.
256         // The goal is to have s2 to sk as accurate as possible considering the fact the sum is
257         // truncated and we don't want the error terms to be included in s2 ... sk, so we need
258         // to solve also for the remainder
259         final double[][] a     = new double[c1.length + 1][c1.length + 1];
260         final double[][] b     = new double[c1.length + 1][y[0].length];
261         final double[]   y0    = y[0];
262         final double[]   yDot0 = yDot[0];
263         for (int i = 1; i < y.length; ++i) {
264 
265             final double di    = t[i] - t[0];
266             final double ratio = di / h;
267             double dikM1Ohk    =  1 / h;
268 
269             // linear coefficients of equations
270             // y(ti) - y(t0) - di y'(t0) and y'(ti) - y'(t0)
271             final double[] aI    = a[2 * i - 2];
272             final double[] aDotI = (2 * i - 1) < a.length ? a[2 * i - 1] : null;
273             for (int j = 0; j < aI.length; ++j) {
274                 dikM1Ohk *= ratio;
275                 aI[j]     = di      * dikM1Ohk;
276                 if (aDotI != null) {
277                     aDotI[j]  = (j + 2) * dikM1Ohk;
278                 }
279             }
280 
281             // expected value of the previous equations
282             final double[] yI    = y[i];
283             final double[] yDotI = yDot[i];
284             final double[] bI    = b[2 * i - 2];
285             final double[] bDotI = (2 * i - 1) < b.length ? b[2 * i - 1] : null;
286             for (int j = 0; j < yI.length; ++j) {
287                 bI[j]    = yI[j] - y0[j] - di * yDot0[j];
288                 if (bDotI != null) {
289                     bDotI[j] = yDotI[j] - yDot0[j];
290                 }
291             }
292 
293         }
294 
295         // solve the linear system to get the best estimate of the Nordsieck vector [s2 ... sk],
296         // with the additional terms s(k+1) and c grabbing the parts after the truncated Taylor expansion
297         final QRDecomposition decomposition = new QRDecomposition(new Array2DRowRealMatrix(a, false));
298         final RealMatrix x = decomposition.getSolver().solve(new Array2DRowRealMatrix(b, false));
299 
300         // extract just the Nordsieck vector [s2 ... sk]
301         final Array2DRowRealMatrix truncatedX = new Array2DRowRealMatrix(x.getRowDimension() - 1, x.getColumnDimension());
302         for (int i = 0; i < truncatedX.getRowDimension(); ++i) {
303             for (int j = 0; j < truncatedX.getColumnDimension(); ++j) {
304                 truncatedX.setEntry(i, j, x.getEntry(i, j));
305             }
306         }
307         return truncatedX;
308 
309     }
310 
311     /** Update the high order scaled derivatives for Adams integrators (phase 1).
312      * <p>The complete update of high order derivatives has a form similar to:
313      * \[
314      * r_{n+1} = (s_1(n) - s_1(n+1)) P^{-1} u + P^{-1} A P r_n
315      * \]
316      * this method computes the P<sup>-1</sup> A P r<sub>n</sub> part.</p>
317      * @param highOrder high order scaled derivatives
318      * (h<sup>2</sup>/2 y'', ... h<sup>k</sup>/k! y(k))
319      * @return updated high order derivatives
320      * @see #updateHighOrderDerivativesPhase2(double[], double[], Array2DRowRealMatrix)
321      */
322     public Array2DRowRealMatrix updateHighOrderDerivativesPhase1(final Array2DRowRealMatrix highOrder) {
323         return update.multiply(highOrder);
324     }
325 
326     /** Update the high order scaled derivatives Adams integrators (phase 2).
327      * <p>The complete update of high order derivatives has a form similar to:
328      * \[
329      * r_{n+1} = (s_1(n) - s_1(n+1)) P^{-1} u + P^{-1} A P r_n
330      * \]
331      * this method computes the (s<sub>1</sub>(n) - s<sub>1</sub>(n+1)) P<sup>-1</sup> u part.</p>
332      * <p>Phase 1 of the update must already have been performed.</p>
333      * @param start first order scaled derivatives at step start
334      * @param end first order scaled derivatives at step end
335      * @param highOrder high order scaled derivatives, will be modified
336      * (h<sup>2</sup>/2 y'', ... h<sup>k</sup>/k! y(k))
337      * @see #updateHighOrderDerivativesPhase1(Array2DRowRealMatrix)
338      */
339     public void updateHighOrderDerivativesPhase2(final double[] start,
340                                                  final double[] end,
341                                                  final Array2DRowRealMatrix highOrder) {
342         final double[][] data = highOrder.getDataRef();
343         for (int i = 0; i < data.length; ++i) {
344             final double[] dataI = data[i];
345             final double c1I = c1[i];
346             for (int j = 0; j < dataI.length; ++j) {
347                 dataI[j] += c1I * (start[j] - end[j]);
348             }
349         }
350     }
351 
352 }