Commons Math example source code file (MultistepIntegrator.java)

This example Commons Math source code file (MultistepIntegrator.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

array2drowrealmatrix, countingdifferentialequations, derivativeexception, derivativeexception, dormandprince853integrator, dormandprince853integrator, initializationcompletedmarkerexception, initializationcompletedmarkerexception, integratorexception, multistepintegrator, multistepintegrator, nordsieckinitializer, realmatrix, string

The Commons Math MultistepIntegrator.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;

import org.apache.commons.math.MathRuntimeException;
import org.apache.commons.math.linear.Array2DRowRealMatrix;
import org.apache.commons.math.linear.RealMatrix;
import org.apache.commons.math.ode.nonstiff.AdaptiveStepsizeIntegrator;
import org.apache.commons.math.ode.nonstiff.DormandPrince853Integrator;
import org.apache.commons.math.ode.sampling.StepHandler;
import org.apache.commons.math.ode.sampling.StepInterpolator;

/**
 * This class is the base class for multistep integrators for Ordinary
 * Differential Equations.
 * <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>Rather than storing several previous steps separately, this implementation uses
 * the Nordsieck vector with higher degrees scaled derivatives all taken at the same
 * step (y<sub>n, s1(n) and rn) where rn is defined as:
 * <pre>
 * r<sub>n = [ s2(n), s3(n) ... sk(n) ]T
 * </pre>
 * (we omit the k index in the notation for clarity)</p>
 * <p>
 * Multistep integrators with Nordsieck representation are highly sensitive to
 * large step changes because when the step is multiplied by a factor a, the
 * k<sup>th component of the Nordsieck vector is multiplied by ak
 * and the last components are the least accurate ones. The default max growth
 * factor is therefore set to a quite low value: 2<sup>1/order.
 * </p>
 *
 * @see org.apache.commons.math.ode.nonstiff.AdamsBashforthIntegrator
 * @see org.apache.commons.math.ode.nonstiff.AdamsMoultonIntegrator
 * @version $Revision: 811827 $ $Date: 2009-09-06 11:32:50 -0400 (Sun, 06 Sep 2009) $
 * @since 2.0
 */
public abstract class MultistepIntegrator extends AdaptiveStepsizeIntegrator {

    /** First scaled derivative (h y'). */
    protected double[] scaled;

    /** Nordsieck matrix of the higher scaled derivatives.
     * <p>(h2/2 y'', h3/6 y''' ..., hk/k! y(k))

     */
    protected Array2DRowRealMatrix nordsieck;

    /** Starter integrator. */
    private FirstOrderIntegrator starter;

    /** Number of steps of the multistep method (excluding the one being computed). */
    private final int nSteps;

    /** Stepsize control exponent. */
    private double exp;

    /** Safety factor for stepsize control. */
    private double safety;

    /** Minimal reduction factor for stepsize control. */
    private double minReduction;

    /** Maximal growth factor for stepsize control. */
    private double maxGrowth;

    /**
     * Build a multistep integrator with the given stepsize bounds.
     * <p>The default starter integrator is set to the {@link
     * DormandPrince853Integrator Dormand-Prince 8(5,3)} integrator with
     * some defaults settings.</p>
     * <p>
     * The default max growth factor is set to a quite low value: 2<sup>1/order.
     * </p>
     * @param name name of the method
     * @param nSteps number of steps of the multistep method
     * (excluding the one being computed)
     * @param order order of the method
     * @param minStep minimal step (must be positive even for backward
     * integration), the last step can be smaller than this
     * @param maxStep maximal step (must be positive even for backward
     * integration)
     * @param scalAbsoluteTolerance allowed absolute error
     * @param scalRelativeTolerance allowed relative error
     */
    protected MultistepIntegrator(final String name, final int nSteps,
                                  final int order,
                                  final double minStep, final double maxStep,
                                  final double scalAbsoluteTolerance,
                                  final double scalRelativeTolerance) {

        super(name, minStep, maxStep, scalAbsoluteTolerance, scalRelativeTolerance);

        if (nSteps <= 0) {
            throw MathRuntimeException.createIllegalArgumentException(
                  "{0} method needs at least one previous point",
                  name);
        }

        starter = new DormandPrince853Integrator(minStep, maxStep,
                                                 scalAbsoluteTolerance,
                                                 scalRelativeTolerance);
        this.nSteps = nSteps;

        exp = -1.0 / order;

        // set the default values of the algorithm control parameters
        setSafety(0.9);
        setMinReduction(0.2);
        setMaxGrowth(Math.pow(2.0, -exp));

    }

    /**
     * Build a multistep integrator with the given stepsize bounds.
     * <p>The default starter integrator is set to the {@link
     * DormandPrince853Integrator Dormand-Prince 8(5,3)} integrator with
     * some defaults settings.</p>
     * <p>
     * The default max growth factor is set to a quite low value: 2<sup>1/order.
     * </p>
     * @param name name of the method
     * @param nSteps number of steps of the multistep method
     * (excluding the one being computed)
     * @param order order of the method
     * @param minStep minimal step (must be positive even for backward
     * integration), the last step can be smaller than this
     * @param maxStep maximal step (must be positive even for backward
     * integration)
     * @param vecAbsoluteTolerance allowed absolute error
     * @param vecRelativeTolerance allowed relative error
     */
    protected MultistepIntegrator(final String name, final int nSteps,
                                  final int order,
                                  final double minStep, final double maxStep,
                                  final double[] vecAbsoluteTolerance,
                                  final double[] vecRelativeTolerance) {
        super(name, minStep, maxStep, vecAbsoluteTolerance, vecRelativeTolerance);
        starter = new DormandPrince853Integrator(minStep, maxStep,
                                                 vecAbsoluteTolerance,
                                                 vecRelativeTolerance);
        this.nSteps = nSteps;

        exp = -1.0 / order;

        // set the default values of the algorithm control parameters
        setSafety(0.9);
        setMinReduction(0.2);
        setMaxGrowth(Math.pow(2.0, -exp));

    }

    /**
     * Get the starter integrator.
     * @return starter integrator
     */
    public ODEIntegrator getStarterIntegrator() {
        return starter;
    }

    /**
     * Set the starter integrator.
     * <p>The various step and event handlers for this starter integrator
     * will be managed automatically by the multi-step integrator. Any
     * user configuration for these elements will be cleared before use.</p>
     * @param starterIntegrator starter integrator
     */
    public void setStarterIntegrator(FirstOrderIntegrator starterIntegrator) {
        this.starter = starterIntegrator;
    }

    /** Start the integration.
     * <p>This method computes one step using the underlying starter integrator,
     * and initializes the Nordsieck vector at step start. The starter integrator
     * purpose is only to establish initial conditions, it does not really change
     * time by itself. The top level multistep integrator remains in charge of
     * handling time propagation and events handling as it will starts its own
     * computation right from the beginning. In a sense, the starter integrator
     * can be seen as a dummy one and so it will never trigger any user event nor
     * call any user step handler.</p>
     * @param t0 initial time
     * @param y0 initial value of the state vector at t0
     * @param t target time for the integration
     * (can be set to a value smaller than <code>t0 for backward integration)
     * @throws IntegratorException if the integrator cannot perform integration
     * @throws DerivativeException this exception is propagated to the caller if
     * the underlying user function triggers one
     */
    protected void start(final double t0, final double[] y0, final double t)
        throws DerivativeException, IntegratorException {

        // make sure NO user event nor user step handler is triggered,
        // this is the task of the top level integrator, not the task
        // of the starter integrator
        starter.clearEventHandlers();
        starter.clearStepHandlers();

        // set up one specific step handler to extract initial Nordsieck vector
        starter.addStepHandler(new NordsieckInitializer(y0.length));

        // start integration, expecting a InitializationCompletedMarkerException
        try {
            starter.integrate(new CountingDifferentialEquations(y0.length),
                              t0, y0, t, new double[y0.length]);
        } catch (DerivativeException de) {
            if (!(de instanceof InitializationCompletedMarkerException)) {
                // this is not the expected nominal interruption of the start integrator
                throw de;
            }
        }

        // remove the specific step handler
        starter.clearStepHandlers();

    }

    /** 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 scaled derivatives at step start
     */
    protected abstract Array2DRowRealMatrix initializeHighOrderDerivatives(final double[] first,
                                                                           final double[][] multistep);

    /** Get the minimal reduction factor for stepsize control.
     * @return minimal reduction factor
     */
    public double getMinReduction() {
        return minReduction;
    }

    /** Set the minimal reduction factor for stepsize control.
     * @param minReduction minimal reduction factor
     */
    public void setMinReduction(final double minReduction) {
        this.minReduction = minReduction;
    }

    /** Get the maximal growth factor for stepsize control.
     * @return maximal growth factor
     */
    public double getMaxGrowth() {
        return maxGrowth;
    }

    /** Set the maximal growth factor for stepsize control.
     * @param maxGrowth maximal growth factor
     */
    public void setMaxGrowth(final double maxGrowth) {
        this.maxGrowth = maxGrowth;
    }

    /** Get the safety factor for stepsize control.
     * @return safety factor
     */
    public double getSafety() {
      return safety;
    }

    /** Set the safety factor for stepsize control.
     * @param safety safety factor
     */
    public void setSafety(final double safety) {
      this.safety = safety;
    }

    /** Compute step grow/shrink factor according to normalized error.
     * @param error normalized error of the current step
     * @return grow/shrink factor for next step
     */
    protected double computeStepGrowShrinkFactor(final double error) {
        return Math.min(maxGrowth, Math.max(minReduction, safety * Math.pow(error, exp)));
    }

    /** Transformer used to convert the first step to Nordsieck representation. */
    public static interface NordsieckTransformer {
        /** 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
         */
        RealMatrix initializeHighOrderDerivatives(double[] first, double[][] multistep);
    }

    /** Specialized step handler storing the first step. */
    private class NordsieckInitializer implements StepHandler {

        /** Problem dimension. */
        private final int n;

        /** Simple constructor.
         * @param n problem dimension
         */
        public NordsieckInitializer(final int n) {
            this.n = n;
        }

        /** {@inheritDoc} */
        public void handleStep(StepInterpolator interpolator, boolean isLast)
            throws DerivativeException {

            final double prev = interpolator.getPreviousTime();
            final double curr = interpolator.getCurrentTime();
            stepStart = prev;
            stepSize  = (curr - prev) / (nSteps + 1);

            // compute the first scaled derivative
            interpolator.setInterpolatedTime(prev);
            scaled = interpolator.getInterpolatedDerivatives().clone();
            for (int j = 0; j < n; ++j) {
                scaled[j] *= stepSize;
            }

            // compute the high order scaled derivatives
            final double[][] multistep = new double[nSteps][];
            for (int i = 1; i <= nSteps; ++i) {
                interpolator.setInterpolatedTime(prev + stepSize * i);
                final double[] msI = interpolator.getInterpolatedDerivatives().clone();
                for (int j = 0; j < n; ++j) {
                    msI[j] *= stepSize;
                }
                multistep[i - 1] = msI;
            }
            nordsieck = initializeHighOrderDerivatives(scaled, multistep);

            // stop the integrator after the first step has been handled
            throw new InitializationCompletedMarkerException();

        }

        /** {@inheritDoc} */
        public boolean requiresDenseOutput() {
            return true;
        }

        /** {@inheritDoc} */
        public void reset() {
            // nothing to do
        }

    }

    /** Marker exception used ONLY to stop the starter integrator after first step. */
    private static class InitializationCompletedMarkerException
        extends DerivativeException {

        /** Serializable version identifier. */
        private static final long serialVersionUID = -4105805787353488365L;

        /** Simple constructor. */
        public InitializationCompletedMarkerException() {
            super((Throwable) null);
        }

    }

    /** Wrapper for differential equations, ensuring start evaluations are counted. */
    private class CountingDifferentialEquations implements FirstOrderDifferentialEquations {

        /** Dimension of the problem. */
        private final int dimension;

        /** Simple constructor.
         * @param dimension dimension of the problem
         */
        public CountingDifferentialEquations(final int dimension) {
            this.dimension = dimension;
        }

        /** {@inheritDoc} */
        public void computeDerivatives(double t, double[] y, double[] dot)
                throws DerivativeException {
            MultistepIntegrator.this.computeDerivatives(t, y, dot);
        }

        /** {@inheritDoc} */
        public int getDimension() {
            return dimension;
        }
    }

}