View Javadoc

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   */
17  
18  package org.apache.commons.math.ode.nonstiff;
19  
20  import java.util.Arrays;
21  import java.util.HashMap;
22  import java.util.Map;
23  
24  import org.apache.commons.math.fraction.BigFraction;
25  import org.apache.commons.math.linear.Array2DRowFieldMatrix;
26  import org.apache.commons.math.linear.Array2DRowRealMatrix;
27  import org.apache.commons.math.linear.DefaultFieldMatrixChangingVisitor;
28  import org.apache.commons.math.linear.FieldDecompositionSolver;
29  import org.apache.commons.math.linear.FieldLUDecompositionImpl;
30  import org.apache.commons.math.linear.FieldMatrix;
31  import org.apache.commons.math.linear.MatrixUtils;
32  
33  /** Transformer to Nordsieck vectors for Adams integrators.
34   * <p>This class i used by {@link AdamsBashforthIntegrator Adams-Bashforth} and
35   * {@link AdamsMoultonIntegrator Adams-Moulton} integrators to convert between
36   * classical representation with several previous first derivatives and Nordsieck
37   * representation with higher order scaled derivatives.</p>
38   *
39   * <p>We define scaled derivatives s<sub>i</sub>(n) at step n as:
40   * <pre>
41   * s<sub>1</sub>(n) = h y'<sub>n</sub> for first derivative
42   * s<sub>2</sub>(n) = h<sup>2</sup>/2 y''<sub>n</sub> for second derivative
43   * s<sub>3</sub>(n) = h<sup>3</sup>/6 y'''<sub>n</sub> for third derivative
44   * ...
45   * s<sub>k</sub>(n) = h<sup>k</sup>/k! y(k)<sub>n</sub> for k<sup>th</sup> derivative
46   * </pre></p>
47   *
48   * <p>With the previous definition, the classical representation of multistep methods
49   * uses first derivatives only, i.e. it handles y<sub>n</sub>, s<sub>1</sub>(n) and
50   * q<sub>n</sub> where q<sub>n</sub> is defined as:
51   * <pre>
52   *   q<sub>n</sub> = [ s<sub>1</sub>(n-1) s<sub>1</sub>(n-2) ... s<sub>1</sub>(n-(k-1)) ]<sup>T</sup>
53   * </pre>
54   * (we omit the k index in the notation for clarity).</p>
55   *
56   * <p>Another possible representation uses the Nordsieck vector with
57   * higher degrees scaled derivatives all taken at the same step, i.e it handles y<sub>n</sub>,
58   * s<sub>1</sub>(n) and r<sub>n</sub>) where r<sub>n</sub> is defined as:
59   * <pre>
60   * r<sub>n</sub> = [ s<sub>2</sub>(n), s<sub>3</sub>(n) ... s<sub>k</sub>(n) ]<sup>T</sup>
61   * </pre>
62   * (here again we omit the k index in the notation for clarity)
63   * </p>
64   *
65   * <p>Taylor series formulas show that for any index offset i, s<sub>1</sub>(n-i) can be
66   * computed from s<sub>1</sub>(n), s<sub>2</sub>(n) ... s<sub>k</sub>(n), the formula being exact
67   * for degree k polynomials.
68   * <pre>
69   * s<sub>1</sub>(n-i) = s<sub>1</sub>(n) + &sum;<sub>j</sub> j (-i)<sup>j-1</sup> s<sub>j</sub>(n)
70   * </pre>
71   * The previous formula can be used with several values for i to compute the transform between
72   * classical representation and Nordsieck vector at step end. The transform between r<sub>n</sub>
73   * and q<sub>n</sub> resulting from the Taylor series formulas above is:
74   * <pre>
75   * q<sub>n</sub> = s<sub>1</sub>(n) u + P r<sub>n</sub>
76   * </pre>
77   * where u is the [ 1 1 ... 1 ]<sup>T</sup> vector and P is the (k-1)&times;(k-1) matrix built
78   * with the j (-i)<sup>j-1</sup> terms:
79   * <pre>
80   *        [  -2   3   -4    5  ... ]
81   *        [  -4  12  -32   80  ... ]
82   *   P =  [  -6  27 -108  405  ... ]
83   *        [  -8  48 -256 1280  ... ]
84   *        [          ...           ]
85   * </pre></p>
86   *
87   * <p>Changing -i into +i in the formula above can be used to compute a similar transform between
88   * classical representation and Nordsieck vector at step start. The resulting matrix is simply
89   * the absolute value of matrix P.</p>
90   *
91   * <p>For {@link AdamsBashforthIntegrator Adams-Bashforth} method, the Nordsieck vector
92   * at step n+1 is computed from the Nordsieck vector at step n as follows:
93   * <ul>
94   *   <li>y<sub>n+1</sub> = y<sub>n</sub> + s<sub>1</sub>(n) + u<sup>T</sup> r<sub>n</sub></li>
95   *   <li>s<sub>1</sub>(n+1) = h f(t<sub>n+1</sub>, y<sub>n+1</sub>)</li>
96   *   <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>
97   * </ul>
98   * where A is a rows shifting matrix (the lower left part is an identity matrix):
99   * <pre>
100  *        [ 0 0   ...  0 0 | 0 ]
101  *        [ ---------------+---]
102  *        [ 1 0   ...  0 0 | 0 ]
103  *    A = [ 0 1   ...  0 0 | 0 ]
104  *        [       ...      | 0 ]
105  *        [ 0 0   ...  1 0 | 0 ]
106  *        [ 0 0   ...  0 1 | 0 ]
107  * </pre></p>
108  *
109  * <p>For {@link AdamsMoultonIntegrator Adams-Moulton} method, the predicted Nordsieck vector
110  * at step n+1 is computed from the Nordsieck vector at step n as follows:
111  * <ul>
112  *   <li>Y<sub>n+1</sub> = y<sub>n</sub> + s<sub>1</sub>(n) + u<sup>T</sup> r<sub>n</sub></li>
113  *   <li>S<sub>1</sub>(n+1) = h f(t<sub>n+1</sub>, Y<sub>n+1</sub>)</li>
114  *   <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>
115  * </ul>
116  * From this predicted vector, the corrected vector is computed as follows:
117  * <ul>
118  *   <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>
119  *   <li>s<sub>1</sub>(n+1) = h f(t<sub>n+1</sub>, y<sub>n+1</sub>)</li>
120  *   <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>
121  * </ul>
122  * where the upper case Y<sub>n+1</sub>, S<sub>1</sub>(n+1) and R<sub>n+1</sub> represent the
123  * predicted states whereas the lower case y<sub>n+1</sub>, s<sub>n+1</sub> and r<sub>n+1</sub>
124  * represent the corrected states.</p>
125  *
126  * <p>We observe that both methods use similar update formulas. In both cases a P<sup>-1</sup>u
127  * vector and a P<sup>-1</sup> A P matrix are used that do not depend on the state,
128  * they only depend on k. This class handles these transformations.</p>
129  *
130  * @version $Revision: 790374 $ $Date: 2009-07-01 16:57:20 -0400 (Wed, 01 Jul 2009) $
131  * @since 2.0
132  */
133 public class AdamsNordsieckTransformer {
134 
135     /** Cache for already computed coefficients. */
136     private static final Map<Integer, AdamsNordsieckTransformer> cache =
137         new HashMap<Integer, AdamsNordsieckTransformer>();
138 
139     /** Initialization matrix for the higher order derivatives wrt y'', y''' ... */
140     private final Array2DRowRealMatrix initialization;
141 
142     /** Update matrix for the higher order derivatives h<sup>2</sup>/2y'', h<sup>3</sup>/6 y''' ... */
143     private final Array2DRowRealMatrix update;
144 
145     /** Update coefficients of the higher order derivatives wrt y'. */
146     private final double[] c1;
147 
148     /** Simple constructor.
149      * @param nSteps number of steps of the multistep method
150      * (excluding the one being computed)
151      */
152     private AdamsNordsieckTransformer(final int nSteps) {
153 
154         // compute exact coefficients
155         FieldMatrix<BigFraction> bigP = buildP(nSteps);
156         FieldDecompositionSolver<BigFraction> pSolver =
157             new FieldLUDecompositionImpl<BigFraction>(bigP).getSolver();
158 
159         BigFraction[] u = new BigFraction[nSteps];
160         Arrays.fill(u, BigFraction.ONE);
161         BigFraction[] bigC1 = pSolver.solve(u);
162 
163         // update coefficients are computed by combining transform from
164         // Nordsieck to multistep, then shifting rows to represent step advance
165         // then applying inverse transform
166         BigFraction[][] shiftedP = bigP.getData();
167         for (int i = shiftedP.length - 1; i > 0; --i) {
168             // shift rows
169             shiftedP[i] = shiftedP[i - 1];
170         }
171         shiftedP[0] = new BigFraction[nSteps];
172         Arrays.fill(shiftedP[0], BigFraction.ZERO);
173         FieldMatrix<BigFraction> bigMSupdate =
174             pSolver.solve(new Array2DRowFieldMatrix<BigFraction>(shiftedP, false));
175 
176         // initialization coefficients, computed from a R matrix = abs(P)
177         bigP.walkInOptimizedOrder(new DefaultFieldMatrixChangingVisitor<BigFraction>(BigFraction.ZERO) {
178             /** {@inheritDoc} */
179             @Override
180             public BigFraction visit(int row, int column, BigFraction value) {
181                 return ((column & 0x1) == 0x1) ? value : value.negate();
182             }
183         });
184         FieldMatrix<BigFraction> bigRInverse =
185             new FieldLUDecompositionImpl<BigFraction>(bigP).getSolver().getInverse();
186 
187         // convert coefficients to double
188         initialization = MatrixUtils.bigFractionMatrixToRealMatrix(bigRInverse);
189         update         = MatrixUtils.bigFractionMatrixToRealMatrix(bigMSupdate);
190         c1             = new double[nSteps];
191         for (int i = 0; i < nSteps; ++i) {
192             c1[i] = bigC1[i].doubleValue();
193         }
194 
195     }
196 
197     /** Get the Nordsieck transformer for a given number of steps.
198      * @param nSteps number of steps of the multistep method
199      * (excluding the one being computed)
200      * @return Nordsieck transformer for the specified number of steps
201      */
202     public static AdamsNordsieckTransformer getInstance(final int nSteps) {
203         synchronized(cache) {
204             AdamsNordsieckTransformer t = cache.get(nSteps);
205             if (t == null) {
206                 t = new AdamsNordsieckTransformer(nSteps);
207                 cache.put(nSteps, t);
208             }
209             return t;
210         }
211     }
212 
213     /** Get the number of steps of the method
214      * (excluding the one being computed).
215      * @return number of steps of the method
216      * (excluding the one being computed)
217      */
218     public int getNSteps() {
219         return c1.length;
220     }
221 
222     /** Build the P matrix.
223      * <p>The P matrix general terms are shifted j (-i)<sup>j-1</sup> terms:
224      * <pre>
225      *        [  -2   3   -4    5  ... ]
226      *        [  -4  12  -32   80  ... ]
227      *   P =  [  -6  27 -108  405  ... ]
228      *        [  -8  48 -256 1280  ... ]
229      *        [          ...           ]
230      * </pre></p>
231      * @param nSteps number of steps of the multistep method
232      * (excluding the one being computed)
233      * @return P matrix
234      */
235     private FieldMatrix<BigFraction> buildP(final int nSteps) {
236 
237         final BigFraction[][] pData = new BigFraction[nSteps][nSteps];
238 
239         for (int i = 0; i < pData.length; ++i) {
240             // build the P matrix elements from Taylor series formulas
241             final BigFraction[] pI = pData[i];
242             final int factor = -(i + 1);
243             int aj = factor;
244             for (int j = 0; j < pI.length; ++j) {
245                 pI[j] = new BigFraction(aj * (j + 2));
246                 aj *= factor;
247             }
248         }
249 
250         return new Array2DRowFieldMatrix<BigFraction>(pData, false);
251 
252     }
253 
254     /** Initialize the high order scaled derivatives at step start.
255      * @param first first scaled derivative at step start
256      * @param multistep scaled derivatives after step start (hy'1, ..., hy'k-1)
257      * will be modified
258      * @return high order derivatives at step start
259      */
260     public Array2DRowRealMatrix initializeHighOrderDerivatives(final double[] first,
261                                                      final double[][] multistep) {
262         for (int i = 0; i < multistep.length; ++i) {
263             final double[] msI = multistep[i];
264             for (int j = 0; j < first.length; ++j) {
265                 msI[j] -= first[j];
266             }
267         }
268         return initialization.multiply(new Array2DRowRealMatrix(multistep, false));
269     }
270 
271     /** Update the high order scaled derivatives for Adams integrators (phase 1).
272      * <p>The complete update of high order derivatives has a form similar to:
273      * <pre>
274      * 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>
275      * </pre>
276      * this method computes the P<sup>-1</sup> A P r<sub>n</sub> part.</p>
277      * @param highOrder high order scaled derivatives
278      * (h<sup>2</sup>/2 y'', ... h<sup>k</sup>/k! y(k))
279      * @return updated high order derivatives
280      * @see #updateHighOrderDerivativesPhase2(double[], double[], Array2DRowRealMatrix)
281      */
282     public Array2DRowRealMatrix updateHighOrderDerivativesPhase1(final Array2DRowRealMatrix highOrder) {
283         return update.multiply(highOrder);
284     }
285 
286     /** Update the high order scaled derivatives Adams integrators (phase 2).
287      * <p>The complete update of high order derivatives has a form similar to:
288      * <pre>
289      * 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>
290      * </pre>
291      * this method computes the (s<sub>1</sub>(n) - s<sub>1</sub>(n+1)) P<sup>-1</sup> u part.</p>
292      * <p>Phase 1 of the update must already have been performed.</p>
293      * @param start first order scaled derivatives at step start
294      * @param end first order scaled derivatives at step end
295      * @param highOrder high order scaled derivatives, will be modified
296      * (h<sup>2</sup>/2 y'', ... h<sup>k</sup>/k! y(k))
297      * @see #updateHighOrderDerivativesPhase1(Array2DRowRealMatrix)
298      */
299     public void updateHighOrderDerivativesPhase2(final double[] start,
300                                                  final double[] end,
301                                                  final Array2DRowRealMatrix highOrder) {
302         final double[][] data = highOrder.getDataRef();
303         for (int i = 0; i < data.length; ++i) {
304             final double[] dataI = data[i];
305             final double c1I = c1[i];
306             for (int j = 0; j < dataI.length; ++j) {
307                 dataI[j] += c1I * (start[j] - end[j]);
308             }
309         }
310     }
311 
312 }