Commons Math example source code file (AdamsNordsieckTransformer.java)

This example Commons Math source code file (AdamsNordsieckTransformer.java) is included in the DevDaily.com "Java Source Code Warehouse" project. The intent of this project is to help you "Learn Java by Example" ^TM.

Java - Commons Math tags/keywords

adamsnordsiecktransformer, adamsnordsiecktransformer, array2drowfieldmatrix, array2drowrealmatrix, array2drowrealmatrix, bigfraction, bigfraction, cache, fielddecompositionsolver, fieldmatrix, fieldmatrix, hashmap, map, override, util

The Commons Math AdamsNordsieckTransformer.java source code

/*
 * Licensed to the Apache Software Foundation (ASF) under one or more
 * contributor license agreements.  See the NOTICE file distributed with
 * this work for additional information regarding copyright ownership.
 * The ASF licenses this file to You under the Apache License, Version 2.0
 * (the "License"); you may not use this file except in compliance with
 * the License.  You may obtain a copy of the License at
 *
 *      http://www.apache.org/licenses/LICENSE-2.0
 *
 * Unless required by applicable law or agreed to in writing, software
 * distributed under the License is distributed on an "AS IS" BASIS,
 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
 * See the License for the specific language governing permissions and
 * limitations under the License.
 */

package org.apache.commons.math.ode.nonstiff;

import java.util.Arrays;
import java.util.HashMap;
import java.util.Map;

import org.apache.commons.math.fraction.BigFraction;
import org.apache.commons.math.linear.Array2DRowFieldMatrix;
import org.apache.commons.math.linear.Array2DRowRealMatrix;
import org.apache.commons.math.linear.DefaultFieldMatrixChangingVisitor;
import org.apache.commons.math.linear.FieldDecompositionSolver;
import org.apache.commons.math.linear.FieldLUDecompositionImpl;
import org.apache.commons.math.linear.FieldMatrix;
import org.apache.commons.math.linear.MatrixUtils;

/** Transformer to Nordsieck vectors for Adams integrators.
 * <p>This class i used by {@link AdamsBashforthIntegrator Adams-Bashforth} and
 * {@link AdamsMoultonIntegrator Adams-Moulton} integrators to convert between
 * classical representation with several previous first derivatives and Nordsieck
 * representation with higher order scaled derivatives.</p>
 *
 * <p>We define scaled derivatives si(n) at step n as:
 * <pre>
 * s<sub>1(n) = h y'n for first derivative
 * s<sub>2(n) = h2/2 y''n for second derivative
 * s<sub>3(n) = h3/6 y'''n for third derivative
 * ...
 * s<sub>k(n) = hk/k! y(k)n for kth derivative
 * </pre>

 *
 * <p>With the previous definition, the classical representation of multistep methods
 * uses first derivatives only, i.e. it handles y<sub>n, s1(n) and
 * q<sub>n where qn is defined as:
 * <pre>
 *   q<sub>n = [ s1(n-1) s1(n-2) ... s1(n-(k-1)) ]T
 * </pre>
 * (we omit the k index in the notation for clarity).</p>
 *
 * <p>Another possible representation uses the Nordsieck vector with
 * higher degrees scaled derivatives all taken at the same step, i.e it handles y<sub>n,
 * s<sub>1(n) and rn) where rn is defined as:
 * <pre>
 * r<sub>n = [ s2(n), s3(n) ... sk(n) ]T
 * </pre>
 * (here again we omit the k index in the notation for clarity)
 * </p>
 *
 * <p>Taylor series formulas show that for any index offset i, s1(n-i) can be
 * computed from s<sub>1(n), s2(n) ... sk(n), the formula being exact
 * for degree k polynomials.
 * <pre>
 * s<sub>1(n-i) = s1(n) + ∑j j (-i)j-1 sj(n)
 * </pre>
 * The previous formula can be used with several values for i to compute the transform between
 * classical representation and Nordsieck vector at step end. The transform between r<sub>n
 * and q<sub>n resulting from the Taylor series formulas above is:
 * <pre>
 * q<sub>n = s1(n) u + P rn
 * </pre>
 * where u is the [ 1 1 ... 1 ]<sup>T vector and P is the (k-1)×(k-1) matrix built
 * with the j (-i)<sup>j-1 terms:
 * <pre>
 *        [  -2   3   -4    5  ... ]
 *        [  -4  12  -32   80  ... ]
 *   P =  [  -6  27 -108  405  ... ]
 *        [  -8  48 -256 1280  ... ]
 *        [          ...           ]
 * </pre>

 *
 * <p>Changing -i into +i in the formula above can be used to compute a similar transform between
 * classical representation and Nordsieck vector at step start. The resulting matrix is simply
 * the absolute value of matrix P.</p>
 *
 * <p>For {@link AdamsBashforthIntegrator Adams-Bashforth} method, the Nordsieck vector
 * at step n+1 is computed from the Nordsieck vector at step n as follows:
 * <ul>
 *   <li>yn+1 = yn + s1(n) + uT rn
 *   <li>s1(n+1) = h f(tn+1, yn+1)
 *   <li>rn+1 = (s1(n) - s1(n+1)) P-1 u + P-1 A P rn
 * </ul>
 * where A is a rows shifting matrix (the lower left part is an identity matrix):
 * <pre>
 *        [ 0 0   ...  0 0 | 0 ]
 *        [ ---------------+---]
 *        [ 1 0   ...  0 0 | 0 ]
 *    A = [ 0 1   ...  0 0 | 0 ]
 *        [       ...      | 0 ]
 *        [ 0 0   ...  1 0 | 0 ]
 *        [ 0 0   ...  0 1 | 0 ]
 * </pre>

 *
 * <p>For {@link AdamsMoultonIntegrator Adams-Moulton} method, the predicted Nordsieck vector
 * at step n+1 is computed from the Nordsieck vector at step n as follows:
 * <ul>
 *   <li>Yn+1 = yn + s1(n) + uT rn
 *   <li>S1(n+1) = h f(tn+1, Yn+1)
 *   <li>Rn+1 = (s1(n) - s1(n+1)) P-1 u + P-1 A P rn
 * </ul>
 * From this predicted vector, the corrected vector is computed as follows:
 * <ul>
 *   <li>yn+1 = yn + S1(n+1) + [ -1 +1 -1 +1 ... ±1 ] rn+1
 *   <li>s1(n+1) = h f(tn+1, yn+1)
 *   <li>rn+1 = Rn+1 + (s1(n+1) - S1(n+1)) P-1 u
 * </ul>
 * where the upper case Y<sub>n+1, S1(n+1) and Rn+1 represent the
 * predicted states whereas the lower case y<sub>n+1, sn+1 and rn+1
 * represent the corrected states.</p>
 *
 * <p>We observe that both methods use similar update formulas. In both cases a P-1u
 * vector and a P<sup>-1 A P matrix are used that do not depend on the state,
 * they only depend on k. This class handles these transformations.</p>
 *
 * @version $Revision: 810196 $ $Date: 2009-09-01 15:47:46 -0400 (Tue, 01 Sep 2009) $
 * @since 2.0
 */
public class AdamsNordsieckTransformer {

    /** Cache for already computed coefficients. */
    private static final Map<Integer, AdamsNordsieckTransformer> CACHE =
        new HashMap<Integer, AdamsNordsieckTransformer>();

    /** Initialization matrix for the higher order derivatives wrt y'', y''' ... */
    private final Array2DRowRealMatrix initialization;

    /** Update matrix for the higher order derivatives h<sup>2/2y'', h3/6 y''' ... */
    private final Array2DRowRealMatrix update;

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

    /** Simple constructor.
     * @param nSteps number of steps of the multistep method
     * (excluding the one being computed)
     */
    private AdamsNordsieckTransformer(final int nSteps) {

        // compute exact coefficients
        FieldMatrix<BigFraction> bigP = buildP(nSteps);
        FieldDecompositionSolver<BigFraction> pSolver =
            new FieldLUDecompositionImpl<BigFraction>(bigP).getSolver();

        BigFraction[] u = new BigFraction[nSteps];
        Arrays.fill(u, BigFraction.ONE);
        BigFraction[] bigC1 = pSolver.solve(u);

        // update coefficients are computed by combining transform from
        // Nordsieck to multistep, then shifting rows to represent step advance
        // then applying inverse transform
        BigFraction[][] shiftedP = bigP.getData();
        for (int i = shiftedP.length - 1; i > 0; --i) {
            // shift rows
            shiftedP[i] = shiftedP[i - 1];
        }
        shiftedP[0] = new BigFraction[nSteps];
        Arrays.fill(shiftedP[0], BigFraction.ZERO);
        FieldMatrix<BigFraction> bigMSupdate =
            pSolver.solve(new Array2DRowFieldMatrix<BigFraction>(shiftedP, false));

        // initialization coefficients, computed from a R matrix = abs(P)
        bigP.walkInOptimizedOrder(new DefaultFieldMatrixChangingVisitor<BigFraction>(BigFraction.ZERO) {
            /** {@inheritDoc} */
            @Override
            public BigFraction visit(int row, int column, BigFraction value) {
                return ((column & 0x1) == 0x1) ? value : value.negate();
            }
        });
        FieldMatrix<BigFraction> bigRInverse =
            new FieldLUDecompositionImpl<BigFraction>(bigP).getSolver().getInverse();

        // convert coefficients to double
        initialization = MatrixUtils.bigFractionMatrixToRealMatrix(bigRInverse);
        update         = MatrixUtils.bigFractionMatrixToRealMatrix(bigMSupdate);
        c1             = new double[nSteps];
        for (int i = 0; i < nSteps; ++i) {
            c1[i] = bigC1[i].doubleValue();
        }

    }

    /** Get the Nordsieck transformer for a given number of steps.
     * @param nSteps number of steps of the multistep method
     * (excluding the one being computed)
     * @return Nordsieck transformer for the specified number of steps
     */
    public static AdamsNordsieckTransformer getInstance(final int nSteps) {
        synchronized(CACHE) {
            AdamsNordsieckTransformer t = CACHE.get(nSteps);
            if (t == null) {
                t = new AdamsNordsieckTransformer(nSteps);
                CACHE.put(nSteps, t);
            }
            return t;
        }
    }

    /** Get the number of steps of the method
     * (excluding the one being computed).
     * @return number of steps of the method
     * (excluding the one being computed)
     */
    public int getNSteps() {
        return c1.length;
    }

    /** Build the P matrix.
     * <p>The P matrix general terms are shifted j (-i)j-1 terms:
     * <pre>
     *        [  -2   3   -4    5  ... ]
     *        [  -4  12  -32   80  ... ]
     *   P =  [  -6  27 -108  405  ... ]
     *        [  -8  48 -256 1280  ... ]
     *        [          ...           ]
     * </pre>

     * @param nSteps number of steps of the multistep method
     * (excluding the one being computed)
     * @return P matrix
     */
    private FieldMatrix<BigFraction> buildP(final int nSteps) {

        final BigFraction[][] pData = new BigFraction[nSteps][nSteps];

        for (int i = 0; i < pData.length; ++i) {
            // build the P matrix elements from Taylor series formulas
            final BigFraction[] pI = pData[i];
            final int factor = -(i + 1);
            int aj = factor;
            for (int j = 0; j < pI.length; ++j) {
                pI[j] = new BigFraction(aj * (j + 2));
                aj *= factor;
            }
        }

        return new Array2DRowFieldMatrix<BigFraction>(pData, false);

    }

    /** Initialize the high order scaled derivatives at step start.
     * @param first first scaled derivative at step start
     * @param multistep scaled derivatives after step start (hy'1, ..., hy'k-1)
     * will be modified
     * @return high order derivatives at step start
     */
    public Array2DRowRealMatrix initializeHighOrderDerivatives(final double[] first,
                                                     final double[][] multistep) {
        for (int i = 0; i < multistep.length; ++i) {
            final double[] msI = multistep[i];
            for (int j = 0; j < first.length; ++j) {
                msI[j] -= first[j];
            }
        }
        return initialization.multiply(new Array2DRowRealMatrix(multistep, false));
    }

    /** Update the high order scaled derivatives for Adams integrators (phase 1).
     * <p>The complete update of high order derivatives has a form similar to:
     * <pre>
     * r<sub>n+1 = (s1(n) - s1(n+1)) P-1 u + P-1 A P rn
     * </pre>
     * this method computes the P<sup>-1 A P rn part.

     * @param highOrder high order scaled derivatives
     * (h<sup>2/2 y'', ... hk/k! y(k))
     * @return updated high order derivatives
     * @see #updateHighOrderDerivativesPhase2(double[], double[], Array2DRowRealMatrix)
     */
    public Array2DRowRealMatrix updateHighOrderDerivativesPhase1(final Array2DRowRealMatrix highOrder) {
        return update.multiply(highOrder);
    }

    /** Update the high order scaled derivatives Adams integrators (phase 2).
     * <p>The complete update of high order derivatives has a form similar to:
     * <pre>
     * r<sub>n+1 = (s1(n) - s1(n+1)) P-1 u + P-1 A P rn
     * </pre>
     * this method computes the (s<sub>1(n) - s1(n+1)) P-1 u part.

     * <p>Phase 1 of the update must already have been performed.

     * @param start first order scaled derivatives at step start
     * @param end first order scaled derivatives at step end
     * @param highOrder high order scaled derivatives, will be modified
     * (h<sup>2/2 y'', ... hk/k! y(k))
     * @see #updateHighOrderDerivativesPhase1(Array2DRowRealMatrix)
     */
    public void updateHighOrderDerivativesPhase2(final double[] start,
                                                 final double[] end,
                                                 final Array2DRowRealMatrix highOrder) {
        final double[][] data = highOrder.getDataRef();
        for (int i = 0; i < data.length; ++i) {
            final double[] dataI = data[i];
            final double c1I = c1[i];
            for (int j = 0; j < dataI.length; ++j) {
                dataI[j] += c1I * (start[j] - end[j]);
            }
        }
    }

}