AdamsIntegrator.java

/*
 * Licensed to the Hipparchus project under one or more
 * contributor license agreements.  See the NOTICE file distributed with
 * this work for additional information regarding copyright ownership.
 * The Hipparchus project 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
 *
 *      https://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.hipparchus.ode.nonstiff;

import org.hipparchus.exception.MathIllegalArgumentException;
import org.hipparchus.exception.MathIllegalStateException;
import org.hipparchus.linear.Array2DRowRealMatrix;
import org.hipparchus.linear.RealMatrix;
import org.hipparchus.ode.EquationsMapper;
import org.hipparchus.ode.ExpandableODE;
import org.hipparchus.ode.LocalizedODEFormats;
import org.hipparchus.ode.MultistepIntegrator;
import org.hipparchus.ode.ODEState;
import org.hipparchus.ode.ODEStateAndDerivative;


/** Base class for {@link AdamsBashforthIntegrator Adams-Bashforth} and
 * {@link AdamsMoultonIntegrator Adams-Moulton} integrators.
 */
public abstract class AdamsIntegrator extends MultistepIntegrator {

    /** Transformer. */
    private final AdamsNordsieckTransformer transformer;

    /**
     * Build an Adams integrator with the given order and step control parameters.
     * @param name name of the method
     * @param nSteps number of steps of the method excluding the one being computed
     * @param order order of the method
     * @param minStep minimal step (sign is irrelevant, regardless of
     * integration direction, forward or backward), the last step can
     * be smaller than this
     * @param maxStep maximal step (sign is irrelevant, regardless of
     * integration direction, forward or backward), the last step can
     * be smaller than this
     * @param scalAbsoluteTolerance allowed absolute error
     * @param scalRelativeTolerance allowed relative error
     * @exception MathIllegalArgumentException if order is 1 or less
     */
    public AdamsIntegrator(final String name, final int nSteps, final int order,
                           final double minStep, final double maxStep,
                           final double scalAbsoluteTolerance,
                           final double scalRelativeTolerance)
        throws MathIllegalArgumentException {
        super(name, nSteps, order, minStep, maxStep,
              scalAbsoluteTolerance, scalRelativeTolerance);
        transformer = AdamsNordsieckTransformer.getInstance(nSteps);
    }

    /**
     * Build an Adams integrator with the given order and step control parameters.
     * @param name name of the method
     * @param nSteps number of steps of the method excluding the one being computed
     * @param order order of the method
     * @param minStep minimal step (sign is irrelevant, regardless of
     * integration direction, forward or backward), the last step can
     * be smaller than this
     * @param maxStep maximal step (sign is irrelevant, regardless of
     * integration direction, forward or backward), the last step can
     * be smaller than this
     * @param vecAbsoluteTolerance allowed absolute error
     * @param vecRelativeTolerance allowed relative error
     * @exception IllegalArgumentException if order is 1 or less
     */
    public AdamsIntegrator(final String name, final int nSteps, final int order,
                           final double minStep, final double maxStep,
                           final double[] vecAbsoluteTolerance,
                           final double[] vecRelativeTolerance)
        throws IllegalArgumentException {
        super(name, nSteps, order, minStep, maxStep,
              vecAbsoluteTolerance, vecRelativeTolerance);
        transformer = AdamsNordsieckTransformer.getInstance(nSteps);
    }

    /** {@inheritDoc} */
    @Override
    public ODEStateAndDerivative integrate(final ExpandableODE equations,
                                           final ODEState initialState,
                                           final double finalTime)
        throws MathIllegalArgumentException, MathIllegalStateException {

        sanityChecks(initialState, finalTime);
        setStepStart(initIntegration(equations, initialState, finalTime));
        final boolean forward = finalTime > initialState.getTime();

        // compute the initial Nordsieck vector using the configured starter integrator
        start(equations, getStepStart(), finalTime);

        // reuse the step that was chosen by the starter integrator
        ODEStateAndDerivative stepEnd   =
                        AdamsStateInterpolator.taylor(equations.getMapper(), getStepStart(),
                                                      getStepStart().getTime() + getStepSize(),
                                                      getStepSize(), scaled, nordsieck);

        // main integration loop
        setIsLastStep(false);
        final double[] y  = getStepStart().getCompleteState();
        do {

            double[] predictedY  = null;
            final double[] predictedScaled = new double[y.length];
            Array2DRowRealMatrix predictedNordsieck = null;
            double error = 10;
            while (error >= 1.0) {

                // predict a first estimate of the state at step end
                predictedY = stepEnd.getCompleteState();

                // evaluate the derivative
                final double[] yDot = computeDerivatives(stepEnd.getTime(), predictedY);

                // predict Nordsieck vector at step end
                for (int j = 0; j < predictedScaled.length; ++j) {
                    predictedScaled[j] = getStepSize() * yDot[j];
                }
                predictedNordsieck = updateHighOrderDerivativesPhase1(nordsieck);
                updateHighOrderDerivativesPhase2(scaled, predictedScaled, predictedNordsieck);

                // evaluate error
                error = errorEstimation(y, stepEnd.getTime(), predictedY, predictedScaled, predictedNordsieck);
                if (Double.isNaN(error)) {
                    throw new MathIllegalStateException(LocalizedODEFormats.NAN_APPEARING_DURING_INTEGRATION,
                                                        stepEnd.getTime());
                }

                if (error >= 1.0) {
                    // reject the step and attempt to reduce error by stepsize control
                    final double factor = computeStepGrowShrinkFactor(error);
                    rescale(getStepSizeHelper().filterStep(getStepSize() * factor, forward, false));
                    stepEnd = AdamsStateInterpolator.taylor(equations.getMapper(), getStepStart(),
                                                            getStepStart().getTime() + getStepSize(),
                                                            getStepSize(),
                                                            scaled,
                                                            nordsieck);

                }
            }

            final AdamsStateInterpolator interpolator =
                            finalizeStep(getStepSize(), predictedY, predictedScaled, predictedNordsieck,
                                         forward, getStepStart(), stepEnd, equations.getMapper());

            // discrete events handling
            setStepStart(acceptStep(interpolator, finalTime));
            scaled    = interpolator.getScaled();
            nordsieck = interpolator.getNordsieck();

            if (!isLastStep()) {

                if (resetOccurred()) {

                    // some events handler has triggered changes that
                    // invalidate the derivatives, we need to restart from scratch
                    start(equations, getStepStart(), finalTime);

                    final double  nextT      = getStepStart().getTime() + getStepSize();
                    final boolean nextIsLast = forward ?
                                               (nextT >= finalTime) :
                                               (nextT <= finalTime);
                    final double hNew = nextIsLast ? finalTime - getStepStart().getTime() : getStepSize();

                    rescale(hNew);
                    System.arraycopy(getStepStart().getCompleteState(), 0, y, 0, y.length);

                } else {

                    // stepsize control for next step
                    final double  factor     = computeStepGrowShrinkFactor(error);
                    final double  scaledH    = getStepSize() * factor;
                    final double  nextT      = getStepStart().getTime() + scaledH;
                    final boolean nextIsLast = forward ?
                                               (nextT >= finalTime) :
                                               (nextT <= finalTime);
                    double hNew = getStepSizeHelper().filterStep(scaledH, forward, nextIsLast);

                    final double  filteredNextT      = getStepStart().getTime() + hNew;
                    final boolean filteredNextIsLast = forward ? (filteredNextT >= finalTime) : (filteredNextT <= finalTime);
                    if (filteredNextIsLast) {
                        hNew = finalTime - getStepStart().getTime();
                    }

                    rescale(hNew);
                    System.arraycopy(predictedY, 0, y, 0, y.length);

                }

                stepEnd = AdamsStateInterpolator.taylor(equations.getMapper(), getStepStart(), getStepStart().getTime() + getStepSize(),
                                                        getStepSize(), scaled, nordsieck);

            }

        } while (!isLastStep());

        final ODEStateAndDerivative finalState = getStepStart();
        setStepStart(null);
        setStepSize(Double.NaN);
        return finalState;

    }

    /** {@inheritDoc} */
    @Override
    protected Array2DRowRealMatrix initializeHighOrderDerivatives(final double h, final double[] t,
                                                                  final double[][] y,
                                                                  final double[][] yDot) {
        return transformer.initializeHighOrderDerivatives(h, t, y, yDot);
    }

    /** Update the high order scaled derivatives for Adams integrators (phase 1).
     * <p>The complete update of high order derivatives has a form similar to:
     * \[
     * r_{n+1} = (s_1(n) - s_1(n+1)) P^{-1} u + P^{-1} A P r_n
     * \]
     * this method computes the P<sup>-1</sup> A P r<sub>n</sub> part.</p>
     * @param highOrder high order scaled derivatives
     * (h<sup>2</sup>/2 y'', ... h<sup>k</sup>/k! y(k))
     * @return updated high order derivatives
     * @see #updateHighOrderDerivativesPhase2(double[], double[], Array2DRowRealMatrix)
     */
    public Array2DRowRealMatrix updateHighOrderDerivativesPhase1(final Array2DRowRealMatrix highOrder) {
        return transformer.updateHighOrderDerivativesPhase1(highOrder);
    }

    /** Update the high order scaled derivatives Adams integrators (phase 2).
     * <p>The complete update of high order derivatives has a form similar to:
     * \[
     * r_{n+1} = (s_1(n) - s_1(n+1)) P^{-1} u + P^{-1} A P r_n
     * \]
     * this method computes the (s<sub>1</sub>(n) - s<sub>1</sub>(n+1)) P<sup>-1</sup> u part.</p>
     * <p>Phase 1 of the update must already have been performed.</p>
     * @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</sup>/2 y'', ... h<sup>k</sup>/k! y(k))
     * @see #updateHighOrderDerivativesPhase1(Array2DRowRealMatrix)
     */
    public void updateHighOrderDerivativesPhase2(final double[] start,
                                                 final double[] end,
                                                 final Array2DRowRealMatrix highOrder) {
        transformer.updateHighOrderDerivativesPhase2(start, end, highOrder);
    }

    /** Estimate error.
     * @param previousState state vector at step start
     * @param predictedTime time at step end
     * @param predictedState predicted state vector at step end
     * @param predictedScaled predicted value of the scaled derivatives at step end
     * @param predictedNordsieck predicted value of the Nordsieck vector at step end
     * @return estimated normalized local discretization error
     * @since 2.0
     */
    protected abstract double errorEstimation(double[] previousState, double predictedTime,
                                              double[] predictedState, double[] predictedScaled,
                                              RealMatrix predictedNordsieck);

    /** Finalize the step.
     * @param stepSize step size used in the scaled and Nordsieck arrays
     * @param predictedState predicted state at end of step
     * @param predictedScaled predicted first scaled derivative
     * @param predictedNordsieck predicted Nordsieck vector
     * @param isForward integration direction indicator
     * @param globalPreviousState start of the global step
     * @param globalCurrentState end of the global step
     * @param equationsMapper mapper for ODE equations primary and secondary components
     * @return step interpolator
     * @since 2.0
     */
    protected abstract AdamsStateInterpolator finalizeStep(double stepSize, double[] predictedState,
                                                           double[] predictedScaled, Array2DRowRealMatrix predictedNordsieck,
                                                           boolean isForward,
                                                           ODEStateAndDerivative globalPreviousState,
                                                           ODEStateAndDerivative globalCurrentState,
                                                           EquationsMapper equationsMapper);

}