AdamsIntegrator.java

  1. /*
  2.  * Licensed to the Hipparchus project under one or more
  3.  * contributor license agreements.  See the NOTICE file distributed with
  4.  * this work for additional information regarding copyright ownership.
  5.  * The Hipparchus project licenses this file to You under the Apache License, Version 2.0
  6.  * (the "License"); you may not use this file except in compliance with
  7.  * the License.  You may obtain a copy of the License at
  8.  *
  9.  *      https://www.apache.org/licenses/LICENSE-2.0
  10.  *
  11.  * Unless required by applicable law or agreed to in writing, software
  12.  * distributed under the License is distributed on an "AS IS" BASIS,
  13.  * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
  14.  * See the License for the specific language governing permissions and
  15.  * limitations under the License.
  16.  */

  18. package org.hipparchus.ode.nonstiff;

  20. import org.hipparchus.exception.MathIllegalArgumentException;
  21. import org.hipparchus.exception.MathIllegalStateException;
  22. import org.hipparchus.linear.Array2DRowRealMatrix;
  23. import org.hipparchus.linear.RealMatrix;
  24. import org.hipparchus.ode.EquationsMapper;
  25. import org.hipparchus.ode.ExpandableODE;
  26. import org.hipparchus.ode.LocalizedODEFormats;
  27. import org.hipparchus.ode.MultistepIntegrator;
  28. import org.hipparchus.ode.ODEState;
  29. import org.hipparchus.ode.ODEStateAndDerivative;
  30. import org.hipparchus.ode.nonstiff.interpolators.AdamsStateInterpolator;


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

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

  41.     /**
  42.      * Build an Adams integrator with the given order and step control parameters.
  43.      * @param name name of the method
  44.      * @param nSteps number of steps of the method excluding the one being computed
  45.      * @param order order of the method
  46.      * @param minStep minimal step (sign is irrelevant, regardless of
  47.      * integration direction, forward or backward), the last step can
  48.      * be smaller than this
  49.      * @param maxStep maximal step (sign is irrelevant, regardless of
  50.      * integration direction, forward or backward), the last step can
  51.      * be smaller than this
  52.      * @param scalAbsoluteTolerance allowed absolute error
  53.      * @param scalRelativeTolerance allowed relative error
  54.      * @exception MathIllegalArgumentException if order is 1 or less
  55.      */
  56.     protected AdamsIntegrator(final String name, final int nSteps, final int order,
  57.                               final double minStep, final double maxStep,
  58.                               final double scalAbsoluteTolerance,
  59.                               final double scalRelativeTolerance)
  60.         throws MathIllegalArgumentException {
  61.         super(name, nSteps, order, minStep, maxStep,
  62.               scalAbsoluteTolerance, scalRelativeTolerance);
  63.         transformer = AdamsNordsieckTransformer.getInstance(nSteps);
  64.     }

  66.     /**
  67.      * Build an Adams integrator with the given order and step control parameters.
  68.      * @param name name of the method
  69.      * @param nSteps number of steps of the method excluding the one being computed
  70.      * @param order order of the method
  71.      * @param minStep minimal step (sign is irrelevant, regardless of
  72.      * integration direction, forward or backward), the last step can
  73.      * be smaller than this
  74.      * @param maxStep maximal step (sign is irrelevant, regardless of
  75.      * integration direction, forward or backward), the last step can
  76.      * be smaller than this
  77.      * @param vecAbsoluteTolerance allowed absolute error
  78.      * @param vecRelativeTolerance allowed relative error
  79.      * @exception IllegalArgumentException if order is 1 or less
  80.      */
  81.     protected AdamsIntegrator(final String name, final int nSteps, final int order,
  82.                               final double minStep, final double maxStep,
  83.                               final double[] vecAbsoluteTolerance, final double[] vecRelativeTolerance)
  84.         throws IllegalArgumentException {
  85.         super(name, nSteps, order, minStep, maxStep,
  86.               vecAbsoluteTolerance, vecRelativeTolerance);
  87.         transformer = AdamsNordsieckTransformer.getInstance(nSteps);
  88.     }

  90.     /** {@inheritDoc} */
  91.     @Override
  92.     public ODEStateAndDerivative integrate(final ExpandableODE equations,
  93.                                            final ODEState initialState,
  94.                                            final double finalTime)
  95.         throws MathIllegalArgumentException, MathIllegalStateException {

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

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

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

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

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

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

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

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

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

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

  151.                 }
  152.             }

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

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

  163.             if (!isLastStep()) {

  165.                 if (resetOccurred()) {

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

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

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

  180.                 } else {

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

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

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

  200.                 }

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

  205.             }

  207.         } while (!isLastStep());

  209.         final ODEStateAndDerivative finalState = getStepStart();
  210.         setStepStart(null);
  211.         setStepSize(Double.NaN);
  212.         return finalState;

  214.     }

  216.     /** {@inheritDoc} */
  217.     @Override
  218.     protected Array2DRowRealMatrix initializeHighOrderDerivatives(final double h, final double[] t,
  219.                                                                   final double[][] y,
  220.                                                                   final double[][] yDot) {
  221.         return transformer.initializeHighOrderDerivatives(h, t, y, yDot);
  222.     }

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

  239.     /** Update the high order scaled derivatives Adams integrators (phase 2).
  240.      * <p>The complete update of high order derivatives has a form similar to:
  241.      * \[
  242.      * r_{n+1} = (s_1(n) - s_1(n+1)) P^{-1} u + P^{-1} A P r_n
  243.      * \]
  244.      * this method computes the (s<sub>1</sub>(n) - s<sub>1</sub>(n+1)) P<sup>-1</sup> u part.</p>
  245.      * <p>Phase 1 of the update must already have been performed.</p>
  246.      * @param start first order scaled derivatives at step start
  247.      * @param end first order scaled derivatives at step end
  248.      * @param highOrder high order scaled derivatives, will be modified
  249.      * (h<sup>2</sup>/2 y'', ... h<sup>k</sup>/k! y(k))
  250.      * @see #updateHighOrderDerivativesPhase1(Array2DRowRealMatrix)
  251.      */
  252.     public void updateHighOrderDerivativesPhase2(final double[] start,
  253.                                                  final double[] end,
  254.                                                  final Array2DRowRealMatrix highOrder) {
  255.         transformer.updateHighOrderDerivativesPhase2(start, end, highOrder);
  256.     }

  258.     /** Estimate error.
  259.      * @param previousState state vector at step start
  260.      * @param predictedTime time at step end
  261.      * @param predictedState predicted state vector at step end
  262.      * @param predictedScaled predicted value of the scaled derivatives at step end
  263.      * @param predictedNordsieck predicted value of the Nordsieck vector at step end
  264.      * @return estimated normalized local discretization error
  265.      * @since 2.0
  266.      */
  267.     protected abstract double errorEstimation(double[] previousState, double predictedTime,
  268.                                               double[] predictedState, double[] predictedScaled,
  269.                                               RealMatrix predictedNordsieck);

  271.     /** Finalize the step.
  272.      * @param stepSize step size used in the scaled and Nordsieck arrays
  273.      * @param predictedState predicted state at end of step
  274.      * @param predictedScaled predicted first scaled derivative
  275.      * @param predictedNordsieck predicted Nordsieck vector
  276.      * @param isForward integration direction indicator
  277.      * @param globalPreviousState start of the global step
  278.      * @param globalCurrentState end of the global step
  279.      * @param equationsMapper mapper for ODE equations primary and secondary components
  280.      * @return step interpolator
  281.      * @since 2.0
  282.      */
  283.     protected abstract AdamsStateInterpolator finalizeStep(double stepSize, double[] predictedState,
  284.                                                            double[] predictedScaled, Array2DRowRealMatrix predictedNordsieck,
  285.                                                            boolean isForward,
  286.                                                            ODEStateAndDerivative globalPreviousState,
  287.                                                            ODEStateAndDerivative globalCurrentState,
  288.                                                            EquationsMapper equationsMapper);

  290. }
Created with JaCoCo 0.8.12.202403310830