AdamsFieldIntegrator.java

  1. /*
  2.  * Licensed to the Apache Software Foundation (ASF) 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 ASF 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. /*
  19.  * This is not the original file distributed by the Apache Software Foundation
  20.  * It has been modified by the Hipparchus project
  21.  */

  23. package org.hipparchus.ode.nonstiff;

  25. import org.hipparchus.Field;
  26. import org.hipparchus.CalculusFieldElement;
  27. import org.hipparchus.exception.MathIllegalArgumentException;
  28. import org.hipparchus.exception.MathIllegalStateException;
  29. import org.hipparchus.linear.Array2DRowFieldMatrix;
  30. import org.hipparchus.linear.FieldMatrix;
  31. import org.hipparchus.ode.FieldEquationsMapper;
  32. import org.hipparchus.ode.FieldExpandableODE;
  33. import org.hipparchus.ode.FieldODEState;
  34. import org.hipparchus.ode.FieldODEStateAndDerivative;
  35. import org.hipparchus.ode.LocalizedODEFormats;
  36. import org.hipparchus.ode.MultistepFieldIntegrator;
  37. import org.hipparchus.ode.nonstiff.interpolators.AdamsFieldStateInterpolator;
  38. import org.hipparchus.util.MathArrays;


  41. /** Base class for {@link AdamsBashforthFieldIntegrator Adams-Bashforth} and
  42.  * {@link AdamsMoultonFieldIntegrator Adams-Moulton} integrators.
  43.  * @param <T> the type of the field elements
  44.  */
  45. public abstract class AdamsFieldIntegrator<T extends CalculusFieldElement<T>> extends MultistepFieldIntegrator<T> {

  47.     /** Transformer. */
  48.     private final AdamsNordsieckFieldTransformer<T> transformer;

  50.     /**
  51.      * Build an Adams integrator with the given order and step control parameters.
  52.      * @param field field to which the time and state vector elements belong
  53.      * @param name name of the method
  54.      * @param nSteps number of steps of the method excluding the one being computed
  55.      * @param order order of the method
  56.      * @param minStep minimal step (sign is irrelevant, regardless of
  57.      * integration direction, forward or backward), the last step can
  58.      * be smaller than this
  59.      * @param maxStep maximal step (sign is irrelevant, regardless of
  60.      * integration direction, forward or backward), the last step can
  61.      * be smaller than this
  62.      * @param scalAbsoluteTolerance allowed absolute error
  63.      * @param scalRelativeTolerance allowed relative error
  64.      * @exception MathIllegalArgumentException if order is 1 or less
  65.      */
  66.     protected AdamsFieldIntegrator(final Field<T> field, final String name, final int nSteps, final int order,
  67.                                    final double minStep, final double maxStep, final double scalAbsoluteTolerance,
  68.                                    final double scalRelativeTolerance)
  69.         throws MathIllegalArgumentException {
  70.         super(field, name, nSteps, order, minStep, maxStep,
  71.               scalAbsoluteTolerance, scalRelativeTolerance);
  72.         transformer = AdamsNordsieckFieldTransformer.getInstance(field, nSteps);
  73.     }

  75.     /**
  76.      * Build an Adams integrator with the given order and step control parameters.
  77.      * @param field field to which the time and state vector elements belong
  78.      * @param name name of the method
  79.      * @param nSteps number of steps of the method excluding the one being computed
  80.      * @param order order of the method
  81.      * @param minStep minimal step (sign is irrelevant, regardless of
  82.      * integration direction, forward or backward), the last step can
  83.      * be smaller than this
  84.      * @param maxStep maximal step (sign is irrelevant, regardless of
  85.      * integration direction, forward or backward), the last step can
  86.      * be smaller than this
  87.      * @param vecAbsoluteTolerance allowed absolute error
  88.      * @param vecRelativeTolerance allowed relative error
  89.      * @exception IllegalArgumentException if order is 1 or less
  90.      */
  91.     protected AdamsFieldIntegrator(final Field<T> field, final String name, final int nSteps, final int order,
  92.                                    final double minStep, final double maxStep, final double[] vecAbsoluteTolerance,
  93.                                    final double[] vecRelativeTolerance)
  94.         throws IllegalArgumentException {
  95.         super(field, name, nSteps, order, minStep, maxStep,
  96.               vecAbsoluteTolerance, vecRelativeTolerance);
  97.         transformer = AdamsNordsieckFieldTransformer.getInstance(field, nSteps);
  98.     }

  100.     /** {@inheritDoc} */
  101.     @Override
  102.     public FieldODEStateAndDerivative<T> integrate(final FieldExpandableODE<T> equations,
  103.                                                    final FieldODEState<T> initialState,
  104.                                                    final T finalTime)
  105.         throws MathIllegalArgumentException, MathIllegalStateException {

  107.         sanityChecks(initialState, finalTime);
  108.         setStepStart(initIntegration(equations, initialState, finalTime));
  109.         final boolean forward = finalTime.subtract(initialState.getTime()).getReal() > 0;

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

  114.         // reuse the step that was chosen by the starter integrator
  115.         FieldODEStateAndDerivative<T> stepStart = getStepStart();
  116.         FieldODEStateAndDerivative<T> stepEnd   =
  117.                         AdamsFieldStateInterpolator.taylor(equations.getMapper(), stepStart,
  118.                                                            stepStart.getTime().add(getStepSize()),
  119.                                                            getStepSize(), scaled, nordsieck);

  121.         // main integration loop
  122.         setIsLastStep(false);
  123.         final T[] y = stepStart.getCompleteState();
  124.         do {

  126.             T[] predictedY = null;
  127.             final T[] predictedScaled = MathArrays.buildArray(getField(), y.length);
  128.             Array2DRowFieldMatrix<T> predictedNordsieck = null;
  129.             double error = 10;
  130.             while (error >= 1.0) {

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

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

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

  145.                 // evaluate error
  146.                 error = errorEstimation(y, stepEnd.getTime(), predictedY, predictedScaled, predictedNordsieck);
  147.                 if (Double.isNaN(error)) {
  148.                     throw new MathIllegalStateException(LocalizedODEFormats.NAN_APPEARING_DURING_INTEGRATION,
  149.                                                         stepEnd.getTime().getReal());
  150.                 }

  152.                 if (error >= 1.0) {
  153.                     // reject the step and attempt to reduce error by stepsize control
  154.                     final double factor = computeStepGrowShrinkFactor(error);
  155.                     rescale(getStepSizeHelper().filterStep(getStepSize().multiply(factor), forward, false));
  156.                     stepEnd = AdamsFieldStateInterpolator.taylor(equations.getMapper(), getStepStart(),
  157.                                                                  getStepStart().getTime().add(getStepSize()),
  158.                                                                  getStepSize(),
  159.                                                                  scaled,
  160.                                                                  nordsieck);

  162.                 }
  163.             }

  165.             final AdamsFieldStateInterpolator<T> interpolator =
  166.                             finalizeStep(getStepSize(), predictedY, predictedScaled, predictedNordsieck,
  167.                                          forward, getStepStart(), stepEnd, equations.getMapper());

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

  174.             if (!isLastStep()) {

  176.                 if (resetOccurred()) {

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

  182.                     final T  nextT      = getStepStart().getTime().add(getStepSize());
  183.                     final boolean nextIsLast = forward ?
  184.                                                nextT.subtract(finalTime).getReal() >= 0 :
  185.                                                nextT.subtract(finalTime).getReal() <= 0;
  186.                     final T hNew = nextIsLast ? finalTime.subtract(getStepStart().getTime()) : getStepSize();

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

  191.                 } else {

  193.                     // stepsize control for next step
  194.                     final double  factor     = computeStepGrowShrinkFactor(error);
  195.                     final T       scaledH    = getStepSize().multiply(factor);
  196.                     final T       nextT      = getStepStart().getTime().add(scaledH);
  197.                     final boolean nextIsLast = forward ?
  198.                                                nextT.subtract(finalTime).getReal() >= 0 :
  199.                                                nextT.subtract(finalTime).getReal() <= 0;
  200.                     T hNew = getStepSizeHelper().filterStep(scaledH, forward, nextIsLast);

  202.                     final T       filteredNextT      = getStepStart().getTime().add(hNew);
  203.                     final boolean filteredNextIsLast = forward ?
  204.                                                        filteredNextT.subtract(finalTime).getReal() >= 0 :
  205.                                                        filteredNextT.subtract(finalTime).getReal() <= 0;
  206.                     if (filteredNextIsLast) {
  207.                         hNew = finalTime.subtract(getStepStart().getTime());
  208.                     }

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

  213.                 }

  215.                 stepEnd = AdamsFieldStateInterpolator.taylor(equations.getMapper(), getStepStart(),
  216.                                                              getStepStart().getTime().add(getStepSize()),
  217.                                                              getStepSize(), scaled, nordsieck);

  219.             }

  221.         } while (!isLastStep());

  223.         final FieldODEStateAndDerivative<T> finalState = getStepStart();
  224.         setStepStart(null);
  225.         setStepSize(null);
  226.         return finalState;

  228.     }

  230.     /** {@inheritDoc} */
  231.     @Override
  232.     protected Array2DRowFieldMatrix<T> initializeHighOrderDerivatives(final T h, final T[] t,
  233.                                                                       final T[][] y,
  234.                                                                       final T[][] yDot) {
  235.         return transformer.initializeHighOrderDerivatives(h, t, y, yDot);
  236.     }

  238.     /** Update the high order scaled derivatives for Adams integrators (phase 1).
  239.      * <p>The complete update of high order derivatives has a form similar to:
  240.      * \[
  241.      * r_{n+1} = (s_1(n) - s_1(n+1)) P^{-1} u + P^{-1} A P r_n
  242.      * \]
  243.      * this method computes the P<sup>-1</sup> A P r<sub>n</sub> part.</p>
  244.      * @param highOrder high order scaled derivatives
  245.      * (h<sup>2</sup>/2 y'', ... h<sup>k</sup>/k! y(k))
  246.      * @return updated high order derivatives
  247.      * @see #updateHighOrderDerivativesPhase2(CalculusFieldElement[], CalculusFieldElement[], Array2DRowFieldMatrix)
  248.      */
  249.     public Array2DRowFieldMatrix<T> updateHighOrderDerivativesPhase1(final Array2DRowFieldMatrix<T> highOrder) {
  250.         return transformer.updateHighOrderDerivativesPhase1(highOrder);
  251.     }

  253.     /** Update the high order scaled derivatives Adams integrators (phase 2).
  254.      * <p>The complete update of high order derivatives has a form similar to:
  255.      * \[
  256.      * r_{n+1} = (s_1(n) - s_1(n+1)) P^{-1} u + P^{-1} A P r_n
  257.      * \]
  258.      * this method computes the (s<sub>1</sub>(n) - s<sub>1</sub>(n+1)) P<sup>-1</sup> u part.</p>
  259.      * <p>Phase 1 of the update must already have been performed.</p>
  260.      * @param start first order scaled derivatives at step start
  261.      * @param end first order scaled derivatives at step end
  262.      * @param highOrder high order scaled derivatives, will be modified
  263.      * (h<sup>2</sup>/2 y'', ... h<sup>k</sup>/k! y(k))
  264.      * @see #updateHighOrderDerivativesPhase1(Array2DRowFieldMatrix)
  265.      */
  266.     public void updateHighOrderDerivativesPhase2(final T[] start, final T[] end,
  267.                                                  final Array2DRowFieldMatrix<T> highOrder) {
  268.         transformer.updateHighOrderDerivativesPhase2(start, end, highOrder);
  269.     }

  271.     /** Estimate error.
  272.      * @param previousState state vector at step start
  273.      * @param predictedTime time at step end
  274.      * @param predictedState predicted state vector at step end
  275.      * @param predictedScaled predicted value of the scaled derivatives at step end
  276.      * @param predictedNordsieck predicted value of the Nordsieck vector at step end
  277.      * @return estimated normalized local discretization error
  278.      * @since 2.0
  279.      */
  280.     protected abstract double errorEstimation(T[] previousState, T predictedTime,
  281.                                               T[] predictedState, T[] predictedScaled,
  282.                                               FieldMatrix<T> predictedNordsieck);

  284.     /** Finalize the step.
  285.      * @param stepSize step size used in the scaled and Nordsieck arrays
  286.      * @param predictedState predicted state at end of step
  287.      * @param predictedScaled predicted first scaled derivative
  288.      * @param predictedNordsieck predicted Nordsieck vector
  289.      * @param isForward integration direction indicator
  290.      * @param globalPreviousState start of the global step
  291.      * @param globalCurrentState end of the global step
  292.      * @param equationsMapper mapper for ODE equations primary and secondary components
  293.      * @return step interpolator
  294.      * @since 2.0
  295.      */
  296.     protected abstract AdamsFieldStateInterpolator<T> finalizeStep(T stepSize, T[] predictedState,
  297.                                                                    T[] predictedScaled, Array2DRowFieldMatrix<T> predictedNordsieck,
  298.                                                                    boolean isForward,
  299.                                                                    FieldODEStateAndDerivative<T> globalPreviousState,
  300.                                                                    FieldODEStateAndDerivative<T> globalCurrentState,
  301.                                                                    FieldEquationsMapper<T> equationsMapper);

  303. }
Created with JaCoCo 0.8.12.202403310830