001    /*
002     * Licensed to the Apache Software Foundation (ASF) under one or more
003     * contributor license agreements.  See the NOTICE file distributed with
004     * this work for additional information regarding copyright ownership.
005     * The ASF licenses this file to You under the Apache License, Version 2.0
006     * (the "License"); you may not use this file except in compliance with
007     * the License.  You may obtain a copy of the License at
008     *
009     *      http://www.apache.org/licenses/LICENSE-2.0
010     *
011     * Unless required by applicable law or agreed to in writing, software
012     * distributed under the License is distributed on an "AS IS" BASIS,
013     * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
014     * See the License for the specific language governing permissions and
015     * limitations under the License.
016     */
017    
018    package org.apache.commons.math3.ode.nonstiff;
019    
020    import java.util.Arrays;
021    import java.util.HashMap;
022    import java.util.Map;
023    
024    import org.apache.commons.math3.fraction.BigFraction;
025    import org.apache.commons.math3.linear.Array2DRowFieldMatrix;
026    import org.apache.commons.math3.linear.Array2DRowRealMatrix;
027    import org.apache.commons.math3.linear.ArrayFieldVector;
028    import org.apache.commons.math3.linear.FieldDecompositionSolver;
029    import org.apache.commons.math3.linear.FieldLUDecomposition;
030    import org.apache.commons.math3.linear.FieldMatrix;
031    import org.apache.commons.math3.linear.MatrixUtils;
032    import org.apache.commons.math3.linear.QRDecomposition;
033    import org.apache.commons.math3.linear.RealMatrix;
034    
035    /** Transformer to Nordsieck vectors for Adams integrators.
036     * <p>This class is used by {@link AdamsBashforthIntegrator Adams-Bashforth} and
037     * {@link AdamsMoultonIntegrator Adams-Moulton} integrators to convert between
038     * classical representation with several previous first derivatives and Nordsieck
039     * representation with higher order scaled derivatives.</p>
040     *
041     * <p>We define scaled derivatives s<sub>i</sub>(n) at step n as:
042     * <pre>
043     * s<sub>1</sub>(n) = h y'<sub>n</sub> for first derivative
044     * s<sub>2</sub>(n) = h<sup>2</sup>/2 y''<sub>n</sub> for second derivative
045     * s<sub>3</sub>(n) = h<sup>3</sup>/6 y'''<sub>n</sub> for third derivative
046     * ...
047     * s<sub>k</sub>(n) = h<sup>k</sup>/k! y<sup>(k)</sup><sub>n</sub> for k<sup>th</sup> derivative
048     * </pre></p>
049     *
050     * <p>With the previous definition, the classical representation of multistep methods
051     * uses first derivatives only, i.e. it handles y<sub>n</sub>, s<sub>1</sub>(n) and
052     * q<sub>n</sub> where q<sub>n</sub> is defined as:
053     * <pre>
054     *   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>
055     * </pre>
056     * (we omit the k index in the notation for clarity).</p>
057     *
058     * <p>Another possible representation uses the Nordsieck vector with
059     * higher degrees scaled derivatives all taken at the same step, i.e it handles y<sub>n</sub>,
060     * s<sub>1</sub>(n) and r<sub>n</sub>) where r<sub>n</sub> is defined as:
061     * <pre>
062     * r<sub>n</sub> = [ s<sub>2</sub>(n), s<sub>3</sub>(n) ... s<sub>k</sub>(n) ]<sup>T</sup>
063     * </pre>
064     * (here again we omit the k index in the notation for clarity)
065     * </p>
066     *
067     * <p>Taylor series formulas show that for any index offset i, s<sub>1</sub>(n-i) can be
068     * computed from s<sub>1</sub>(n), s<sub>2</sub>(n) ... s<sub>k</sub>(n), the formula being exact
069     * for degree k polynomials.
070     * <pre>
071     * s<sub>1</sub>(n-i) = s<sub>1</sub>(n) + &sum;<sub>j&gt;1</sub> j (-i)<sup>j-1</sup> s<sub>j</sub>(n)
072     * </pre>
073     * The previous formula can be used with several values for i to compute the transform between
074     * classical representation and Nordsieck vector at step end. The transform between r<sub>n</sub>
075     * and q<sub>n</sub> resulting from the Taylor series formulas above is:
076     * <pre>
077     * q<sub>n</sub> = s<sub>1</sub>(n) u + P r<sub>n</sub>
078     * </pre>
079     * where u is the [ 1 1 ... 1 ]<sup>T</sup> vector and P is the (k-1)&times;(k-1) matrix built
080     * with the j (-i)<sup>j-1</sup> terms:
081     * <pre>
082     *        [  -2   3   -4    5  ... ]
083     *        [  -4  12  -32   80  ... ]
084     *   P =  [  -6  27 -108  405  ... ]
085     *        [  -8  48 -256 1280  ... ]
086     *        [          ...           ]
087     * </pre></p>
088     *
089     * <p>Changing -i into +i in the formula above can be used to compute a similar transform between
090     * classical representation and Nordsieck vector at step start. The resulting matrix is simply
091     * the absolute value of matrix P.</p>
092     *
093     * <p>For {@link AdamsBashforthIntegrator Adams-Bashforth} method, the Nordsieck vector
094     * at step n+1 is computed from the Nordsieck vector at step n as follows:
095     * <ul>
096     *   <li>y<sub>n+1</sub> = y<sub>n</sub> + s<sub>1</sub>(n) + u<sup>T</sup> r<sub>n</sub></li>
097     *   <li>s<sub>1</sub>(n+1) = h f(t<sub>n+1</sub>, y<sub>n+1</sub>)</li>
098     *   <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>
099     * </ul>
100     * where A is a rows shifting matrix (the lower left part is an identity matrix):
101     * <pre>
102     *        [ 0 0   ...  0 0 | 0 ]
103     *        [ ---------------+---]
104     *        [ 1 0   ...  0 0 | 0 ]
105     *    A = [ 0 1   ...  0 0 | 0 ]
106     *        [       ...      | 0 ]
107     *        [ 0 0   ...  1 0 | 0 ]
108     *        [ 0 0   ...  0 1 | 0 ]
109     * </pre></p>
110     *
111     * <p>For {@link AdamsMoultonIntegrator Adams-Moulton} method, the predicted Nordsieck vector
112     * at step n+1 is computed from the Nordsieck vector at step n as follows:
113     * <ul>
114     *   <li>Y<sub>n+1</sub> = y<sub>n</sub> + s<sub>1</sub>(n) + u<sup>T</sup> r<sub>n</sub></li>
115     *   <li>S<sub>1</sub>(n+1) = h f(t<sub>n+1</sub>, Y<sub>n+1</sub>)</li>
116     *   <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>
117     * </ul>
118     * From this predicted vector, the corrected vector is computed as follows:
119     * <ul>
120     *   <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>
121     *   <li>s<sub>1</sub>(n+1) = h f(t<sub>n+1</sub>, y<sub>n+1</sub>)</li>
122     *   <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>
123     * </ul>
124     * where the upper case Y<sub>n+1</sub>, S<sub>1</sub>(n+1) and R<sub>n+1</sub> represent the
125     * predicted states whereas the lower case y<sub>n+1</sub>, s<sub>n+1</sub> and r<sub>n+1</sub>
126     * represent the corrected states.</p>
127     *
128     * <p>We observe that both methods use similar update formulas. In both cases a P<sup>-1</sup>u
129     * vector and a P<sup>-1</sup> A P matrix are used that do not depend on the state,
130     * they only depend on k. This class handles these transformations.</p>
131     *
132     * @version $Id: AdamsNordsieckTransformer.java 1416643 2012-12-03 19:37:14Z tn $
133     * @since 2.0
134     */
135    public class AdamsNordsieckTransformer {
136    
137        /** Cache for already computed coefficients. */
138        private static final Map<Integer, AdamsNordsieckTransformer> CACHE =
139            new HashMap<Integer, AdamsNordsieckTransformer>();
140    
141        /** Update matrix for the higher order derivatives h<sup>2</sup>/2y'', h<sup>3</sup>/6 y''' ... */
142        private final Array2DRowRealMatrix update;
143    
144        /** Update coefficients of the higher order derivatives wrt y'. */
145        private final double[] c1;
146    
147        /** Simple constructor.
148         * @param nSteps number of steps of the multistep method
149         * (excluding the one being computed)
150         */
151        private AdamsNordsieckTransformer(final int nSteps) {
152    
153            // compute exact coefficients
154            FieldMatrix<BigFraction> bigP = buildP(nSteps);
155            FieldDecompositionSolver<BigFraction> pSolver =
156                new FieldLUDecomposition<BigFraction>(bigP).getSolver();
157    
158            BigFraction[] u = new BigFraction[nSteps];
159            Arrays.fill(u, BigFraction.ONE);
160            BigFraction[] bigC1 = pSolver
161                .solve(new ArrayFieldVector<BigFraction>(u, false)).toArray();
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            // convert coefficients to double
177            update         = MatrixUtils.bigFractionMatrixToRealMatrix(bigMSupdate);
178            c1             = new double[nSteps];
179            for (int i = 0; i < nSteps; ++i) {
180                c1[i] = bigC1[i].doubleValue();
181            }
182    
183        }
184    
185        /** Get the Nordsieck transformer for a given number of steps.
186         * @param nSteps number of steps of the multistep method
187         * (excluding the one being computed)
188         * @return Nordsieck transformer for the specified number of steps
189         */
190        public static AdamsNordsieckTransformer getInstance(final int nSteps) {
191            synchronized(CACHE) {
192                AdamsNordsieckTransformer t = CACHE.get(nSteps);
193                if (t == null) {
194                    t = new AdamsNordsieckTransformer(nSteps);
195                    CACHE.put(nSteps, t);
196                }
197                return t;
198            }
199        }
200    
201        /** Get the number of steps of the method
202         * (excluding the one being computed).
203         * @return number of steps of the method
204         * (excluding the one being computed)
205         */
206        public int getNSteps() {
207            return c1.length;
208        }
209    
210        /** Build the P matrix.
211         * <p>The P matrix general terms are shifted j (-i)<sup>j-1</sup> terms:
212         * <pre>
213         *        [  -2   3   -4    5  ... ]
214         *        [  -4  12  -32   80  ... ]
215         *   P =  [  -6  27 -108  405  ... ]
216         *        [  -8  48 -256 1280  ... ]
217         *        [          ...           ]
218         * </pre></p>
219         * @param nSteps number of steps of the multistep method
220         * (excluding the one being computed)
221         * @return P matrix
222         */
223        private FieldMatrix<BigFraction> buildP(final int nSteps) {
224    
225            final BigFraction[][] pData = new BigFraction[nSteps][nSteps];
226    
227            for (int i = 0; i < pData.length; ++i) {
228                // build the P matrix elements from Taylor series formulas
229                final BigFraction[] pI = pData[i];
230                final int factor = -(i + 1);
231                int aj = factor;
232                for (int j = 0; j < pI.length; ++j) {
233                    pI[j] = new BigFraction(aj * (j + 2));
234                    aj *= factor;
235                }
236            }
237    
238            return new Array2DRowFieldMatrix<BigFraction>(pData, false);
239    
240        }
241    
242        /** Initialize the high order scaled derivatives at step start.
243         * @param h step size to use for scaling
244         * @param t first steps times
245         * @param y first steps states
246         * @param yDot first steps derivatives
247         * @return Nordieck vector at first step (h<sup>2</sup>/2 y''<sub>n</sub>,
248         * h<sup>3</sup>/6 y'''<sub>n</sub> ... h<sup>k</sup>/k! y<sup>(k)</sup><sub>n</sub>)
249         */
250        public Array2DRowRealMatrix initializeHighOrderDerivatives(final double h, final double[] t,
251                                                                   final double[][] y,
252                                                                   final double[][] yDot) {
253    
254            // using Taylor series with di = ti - t0, we get:
255            //  y(ti)  - y(t0)  - di y'(t0) =   di^2 / h^2 s2 + ... +   di^k     / h^k sk + O(h^(k+1))
256            //  y'(ti) - y'(t0)             = 2 di   / h^2 s2 + ... + k di^(k-1) / h^k sk + O(h^k)
257            // we write these relations for i = 1 to i= n-1 as a set of 2(n-1) linear
258            // equations depending on the Nordsieck vector [s2 ... sk]
259            final double[][] a     = new double[2 * (y.length - 1)][c1.length];
260            final double[][] b     = new double[2 * (y.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 = a[2 * i - 1];
273                for (int j = 0; j < aI.length; ++j) {
274                    dikM1Ohk *= ratio;
275                    aI[j]     = di      * dikM1Ohk;
276                    aDotI[j]  = (j + 2) * dikM1Ohk;
277                }
278    
279                // expected value of the previous equations
280                final double[] yI    = y[i];
281                final double[] yDotI = yDot[i];
282                final double[] bI    = b[2 * i - 2];
283                final double[] bDotI = b[2 * i - 1];
284                for (int j = 0; j < yI.length; ++j) {
285                    bI[j]    = yI[j] - y0[j] - di * yDot0[j];
286                    bDotI[j] = yDotI[j] - yDot0[j];
287                }
288    
289            }
290    
291            // solve the rectangular system in the least square sense
292            // to get the best estimate of the Nordsieck vector [s2 ... sk]
293            QRDecomposition decomposition;
294            decomposition = new QRDecomposition(new Array2DRowRealMatrix(a, false));
295            RealMatrix x = decomposition.getSolver().solve(new Array2DRowRealMatrix(b, false));
296            return new Array2DRowRealMatrix(x.getData(), false);
297        }
298    
299        /** Update the high order scaled derivatives for Adams integrators (phase 1).
300         * <p>The complete update of high order derivatives has a form similar to:
301         * <pre>
302         * 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>
303         * </pre>
304         * this method computes the P<sup>-1</sup> A P r<sub>n</sub> part.</p>
305         * @param highOrder high order scaled derivatives
306         * (h<sup>2</sup>/2 y'', ... h<sup>k</sup>/k! y(k))
307         * @return updated high order derivatives
308         * @see #updateHighOrderDerivativesPhase2(double[], double[], Array2DRowRealMatrix)
309         */
310        public Array2DRowRealMatrix updateHighOrderDerivativesPhase1(final Array2DRowRealMatrix highOrder) {
311            return update.multiply(highOrder);
312        }
313    
314        /** Update the high order scaled derivatives Adams integrators (phase 2).
315         * <p>The complete update of high order derivatives has a form similar to:
316         * <pre>
317         * 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>
318         * </pre>
319         * this method computes the (s<sub>1</sub>(n) - s<sub>1</sub>(n+1)) P<sup>-1</sup> u part.</p>
320         * <p>Phase 1 of the update must already have been performed.</p>
321         * @param start first order scaled derivatives at step start
322         * @param end first order scaled derivatives at step end
323         * @param highOrder high order scaled derivatives, will be modified
324         * (h<sup>2</sup>/2 y'', ... h<sup>k</sup>/k! y(k))
325         * @see #updateHighOrderDerivativesPhase1(Array2DRowRealMatrix)
326         */
327        public void updateHighOrderDerivativesPhase2(final double[] start,
328                                                     final double[] end,
329                                                     final Array2DRowRealMatrix highOrder) {
330            final double[][] data = highOrder.getDataRef();
331            for (int i = 0; i < data.length; ++i) {
332                final double[] dataI = data[i];
333                final double c1I = c1[i];
334                for (int j = 0; j < dataI.length; ++j) {
335                    dataI[j] += c1I * (start[j] - end[j]);
336                }
337            }
338        }
339    
340    }