    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.math.ode.nonstiff;
    019
    020    import java.util.Arrays;
    021    import java.util.HashMap;
    022    import java.util.Map;
    023
    024    import org.apache.commons.math.fraction.BigFraction;
    025    import org.apache.commons.math.linear.Array2DRowFieldMatrix;
    026    import org.apache.commons.math.linear.Array2DRowRealMatrix;
    027    import org.apache.commons.math.linear.DefaultFieldMatrixChangingVisitor;
    028    import org.apache.commons.math.linear.FieldDecompositionSolver;
    029    import org.apache.commons.math.linear.FieldLUDecompositionImpl;
    030    import org.apache.commons.math.linear.FieldMatrix;
    031    import org.apache.commons.math.linear.MatrixUtils;
    032
    033    /** Transformer to Nordsieck vectors for Adams integrators.
    034     * <p>This class i used by {@link AdamsBashforthIntegrator Adams-Bashforth} and
    035     * {@link AdamsMoultonIntegrator Adams-Moulton} integrators to convert between
    036     * classical representation with several previous first derivatives and Nordsieck
    037     * representation with higher order scaled derivatives.</p>
    038     *
    039     * <p>We define scaled derivatives s<sub>i</sub>(n) at step n as:
    040     * <pre>
    041     * s<sub>1</sub>(n) = h y'<sub>n</sub> for first derivative
    042     * s<sub>2</sub>(n) = h<sup>2</sup>/2 y''<sub>n</sub> for second derivative
    043     * s<sub>3</sub>(n) = h<sup>3</sup>/6 y'''<sub>n</sub> for third derivative
    044     * ...
    045     * s<sub>k</sub>(n) = h<sup>k</sup>/k! y(k)<sub>n</sub> for k<sup>th</sup> derivative
    046     * </pre></p>
    047     *
    048     * <p>With the previous definition, the classical representation of multistep methods
    049     * uses first derivatives only, i.e. it handles y<sub>n</sub>, s<sub>1</sub>(n) and
    050     * q<sub>n</sub> where q<sub>n</sub> is defined as:
    051     * <pre>
    052     *   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>
    053     * </pre>
    054     * (we omit the k index in the notation for clarity).</p>
    055     *
    056     * <p>Another possible representation uses the Nordsieck vector with
    057     * higher degrees scaled derivatives all taken at the same step, i.e it handles y<sub>n</sub>,
    058     * s<sub>1</sub>(n) and r<sub>n</sub>) where r<sub>n</sub> is defined as:
    059     * <pre>
    060     * r<sub>n</sub> = [ s<sub>2</sub>(n), s<sub>3</sub>(n) ... s<sub>k</sub>(n) ]<sup>T</sup>
    061     * </pre>
    062     * (here again we omit the k index in the notation for clarity)
    063     * </p>
    064     *
    065     * <p>Taylor series formulas show that for any index offset i, s<sub>1</sub>(n-i) can be
    066     * computed from s<sub>1</sub>(n), s<sub>2</sub>(n) ... s<sub>k</sub>(n), the formula being exact
    067     * for degree k polynomials.
    068     * <pre>
    069     * 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)
    070     * </pre>
    071     * The previous formula can be used with several values for i to compute the transform between
    072     * classical representation and Nordsieck vector at step end. The transform between r<sub>n</sub>
    073     * and q<sub>n</sub> resulting from the Taylor series formulas above is:
    074     * <pre>
    075     * q<sub>n</sub> = s<sub>1</sub>(n) u + P r<sub>n</sub>
    076     * </pre>
    077     * where u is the [ 1 1 ... 1 ]<sup>T</sup> vector and P is the (k-1)&times;(k-1) matrix built
    078     * with the j (-i)<sup>j-1</sup> terms:
    079     * <pre>
    080     *        [  -2   3   -4    5  ... ]
    081     *        [  -4  12  -32   80  ... ]
    082     *   P =  [  -6  27 -108  405  ... ]
    083     *        [  -8  48 -256 1280  ... ]
    084     *        [          ...           ]
    085     * </pre></p>
    086     *
    087     * <p>Changing -i into +i in the formula above can be used to compute a similar transform between
    088     * classical representation and Nordsieck vector at step start. The resulting matrix is simply
    089     * the absolute value of matrix P.</p>
    090     *
    091     * <p>For {@link AdamsBashforthIntegrator Adams-Bashforth} method, the Nordsieck vector
    092     * at step n+1 is computed from the Nordsieck vector at step n as follows:
    093     * <ul>
    094     *   <li>y<sub>n+1</sub> = y<sub>n</sub> + s<sub>1</sub>(n) + u<sup>T</sup> r<sub>n</sub></li>
    095     *   <li>s<sub>1</sub>(n+1) = h f(t<sub>n+1</sub>, y<sub>n+1</sub>)</li>
    096     *   <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>
    097     * </ul>
    098     * where A is a rows shifting matrix (the lower left part is an identity matrix):
    099     * <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: 810196 $ $Date: 2009-09-01 21:47:46 +0200 (mar. 01 sept. 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    }