AdaptiveStepsizeFieldIntegrator.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.CalculusFieldElement;
  26. import org.hipparchus.Field;
  27. import org.hipparchus.exception.MathIllegalArgumentException;
  28. import org.hipparchus.exception.MathIllegalStateException;
  29. import org.hipparchus.ode.AbstractFieldIntegrator;
  30. import org.hipparchus.ode.FieldODEState;
  31. import org.hipparchus.ode.FieldODEStateAndDerivative;
  32. import org.hipparchus.util.FastMath;
  33. import org.hipparchus.util.MathArrays;

  35. /**
  36.  * This abstract class holds the common part of all adaptive
  37.  * stepsize integrators for Ordinary Differential Equations.
  38.  *
  39.  * <p>These algorithms perform integration with stepsize control, which
  40.  * means the user does not specify the integration step but rather a
  41.  * tolerance on error. The error threshold is computed as
  42.  * </p>
  43.  * <pre>
  44.  * threshold_i = absTol_i + relTol_i * max (abs (ym), abs (ym+1))
  45.  * </pre>
  46.  * <p>
  47.  * where absTol_i is the absolute tolerance for component i of the
  48.  * state vector and relTol_i is the relative tolerance for the same
  49.  * component. The user can also use only two scalar values absTol and
  50.  * relTol which will be used for all components.
  51.  * </p>
  52.  * <p>
  53.  * Note that <em>only</em> the {@link FieldODEState#getPrimaryState() main part}
  54.  * of the state vector is used for stepsize control. The {@link
  55.  * FieldODEState#getSecondaryState(int) secondary parts} of the state
  56.  * vector are explicitly ignored for stepsize control.
  57.  * </p>
  58.  *
  59.  * <p>If the estimated error for ym+1 is such that</p>
  60.  * <pre>
  61.  * sqrt((sum (errEst_i / threshold_i)^2 ) / n) &lt; 1
  62.  * </pre>
  63.  *
  64.  * <p>(where n is the main set dimension) then the step is accepted,
  65.  * otherwise the step is rejected and a new attempt is made with a new
  66.  * stepsize.</p>
  67.  *
  68.  * @param <T> the type of the field elements
  69.  *
  70.  */

  72. public abstract class AdaptiveStepsizeFieldIntegrator<T extends CalculusFieldElement<T>>
  73.     extends AbstractFieldIntegrator<T> {

  75.     /** Helper for step size control. */
  76.     private StepsizeHelper stepsizeHelper;

  78.     /** Build an integrator with the given stepsize bounds.
  79.      * The default step handler does nothing.
  80.      * @param field field to which the time and state vector elements belong
  81.      * @param name name of the method
  82.      * @param minStep minimal step (sign is irrelevant, regardless of
  83.      * integration direction, forward or backward), the last step can
  84.      * be smaller than this
  85.      * @param maxStep maximal step (sign is irrelevant, regardless of
  86.      * integration direction, forward or backward), the last step can
  87.      * be smaller than this
  88.      * @param scalAbsoluteTolerance allowed absolute error
  89.      * @param scalRelativeTolerance allowed relative error
  90.      */
  91.     protected AdaptiveStepsizeFieldIntegrator(final Field<T> field, final String name,
  92.                                               final double minStep, final double maxStep,
  93.                                               final double scalAbsoluteTolerance,
  94.                                               final double scalRelativeTolerance) {
  95.         super(field, name);
  96.         stepsizeHelper = new StepsizeHelper(minStep, maxStep, scalAbsoluteTolerance, scalRelativeTolerance);
  97.         resetInternalState();
  98.     }

  100.     /** Build an integrator with the given stepsize bounds.
  101.      * The default step handler does nothing.
  102.      * @param field field to which the time and state vector elements belong
  103.      * @param name name of the method
  104.      * @param minStep minimal step (sign is irrelevant, regardless of
  105.      * integration direction, forward or backward), the last step can
  106.      * be smaller than this
  107.      * @param maxStep maximal step (sign is irrelevant, regardless of
  108.      * integration direction, forward or backward), the last step can
  109.      * be smaller than this
  110.      * @param vecAbsoluteTolerance allowed absolute error
  111.      * @param vecRelativeTolerance allowed relative error
  112.      */
  113.     protected AdaptiveStepsizeFieldIntegrator(final Field<T> field, final String name,
  114.                                               final double minStep, final double maxStep,
  115.                                               final double[] vecAbsoluteTolerance,
  116.                                               final double[] vecRelativeTolerance) {
  117.         super(field, name);
  118.         stepsizeHelper = new StepsizeHelper(minStep, maxStep, vecAbsoluteTolerance, vecRelativeTolerance);
  119.         resetInternalState();
  120.     }

  122.     /** Set the adaptive step size control parameters.
  123.      * <p>
  124.      * A side effect of this method is to also reset the initial
  125.      * step so it will be automatically computed by the integrator
  126.      * if {@link #setInitialStepSize(double) setInitialStepSize}
  127.      * is not called by the user.
  128.      * </p>
  129.      * @param minimalStep minimal step (must be positive even for backward
  130.      * integration), the last step can be smaller than this
  131.      * @param maximalStep maximal step (must be positive even for backward
  132.      * integration)
  133.      * @param absoluteTolerance allowed absolute error
  134.      * @param relativeTolerance allowed relative error
  135.      */
  136.     public void setStepSizeControl(final double minimalStep, final double maximalStep,
  137.                                    final double absoluteTolerance, final double relativeTolerance) {
  138.         stepsizeHelper = new StepsizeHelper(minimalStep, maximalStep, absoluteTolerance, relativeTolerance);
  139.     }

  141.     /** Set the adaptive step size control parameters.
  142.      * <p>
  143.      * A side effect of this method is to also reset the initial
  144.      * step so it will be automatically computed by the integrator
  145.      * if {@link #setInitialStepSize(double) setInitialStepSize}
  146.      * is not called by the user.
  147.      * </p>
  148.      * @param minimalStep minimal step (must be positive even for backward
  149.      * integration), the last step can be smaller than this
  150.      * @param maximalStep maximal step (must be positive even for backward
  151.      * integration)
  152.      * @param absoluteTolerance allowed absolute error
  153.      * @param relativeTolerance allowed relative error
  154.      */
  155.     public void setStepSizeControl(final double minimalStep, final double maximalStep,
  156.                                    final double[] absoluteTolerance,
  157.                                    final double[] relativeTolerance) {
  158.         stepsizeHelper = new StepsizeHelper(minimalStep, maximalStep, absoluteTolerance, relativeTolerance);
  159.     }

  161.     /** Get the stepsize helper.
  162.      * @return stepsize helper
  163.      * @since 2.0
  164.      */
  165.     protected StepsizeHelper getStepSizeHelper() {
  166.         return stepsizeHelper;
  167.     }

  169.     /** Set the initial step size.
  170.      * <p>This method allows the user to specify an initial positive
  171.      * step size instead of letting the integrator guess it by
  172.      * itself. If this method is not called before integration is
  173.      * started, the initial step size will be estimated by the
  174.      * integrator.</p>
  175.      * @param initialStepSize initial step size to use (must be positive even
  176.      * for backward integration ; providing a negative value or a value
  177.      * outside of the min/max step interval will lead the integrator to
  178.      * ignore the value and compute the initial step size by itself)
  179.      */
  180.     public void setInitialStepSize(final double initialStepSize) {
  181.         stepsizeHelper.setInitialStepSize(initialStepSize);
  182.     }

  184.     /** {@inheritDoc} */
  185.     @Override
  186.     protected void sanityChecks(final FieldODEState<T> initialState, final T t)
  187.         throws MathIllegalArgumentException {
  188.         super.sanityChecks(initialState, t);
  189.         stepsizeHelper.setMainSetDimension(initialState.getPrimaryStateDimension());
  190.     }

  192.     /** Initialize the integration step.
  193.      * @param forward forward integration indicator
  194.      * @param order order of the method
  195.      * @param scale scaling vector for the state vector (can be shorter than state vector)
  196.      * @param state0 state at integration start time
  197.      * @return first integration step
  198.      * @exception MathIllegalStateException if the number of functions evaluations is exceeded
  199.      * @exception MathIllegalArgumentException if arrays dimensions do not match equations settings
  200.      */
  201.     public double initializeStep(final boolean forward, final int order, final T[] scale,
  202.                                  final FieldODEStateAndDerivative<T> state0)
  203.         throws MathIllegalArgumentException, MathIllegalStateException {

  205.         if (stepsizeHelper.getInitialStep() > 0) {
  206.             // use the user provided value
  207.             return forward ? stepsizeHelper.getInitialStep() : -stepsizeHelper.getInitialStep();
  208.         }

  210.         // very rough first guess : h = 0.01 * ||y/scale|| / ||y'/scale||
  211.         // this guess will be used to perform an Euler step
  212.         final T[] y0    = state0.getCompleteState();
  213.         final T[] yDot0 = state0.getCompleteDerivative();
  214.         double yOnScale2     = 0;
  215.         double yDotOnScale2  = 0;
  216.         for (int j = 0; j < scale.length; ++j) {
  217.             final double ratio    = y0[j].getReal() / scale[j].getReal();
  218.             yOnScale2            += ratio * ratio;
  219.             final double ratioDot = yDot0[j].getReal() / scale[j].getReal();
  220.             yDotOnScale2         += ratioDot * ratioDot;
  221.         }

  223.         double h = ((yOnScale2 < 1.0e-10) || (yDotOnScale2 < 1.0e-10)) ?
  224.                    1.0e-6 : (0.01 * FastMath.sqrt(yOnScale2 / yDotOnScale2));
  225.         if (h > getMaxStep()) {
  226.             h = getMaxStep();
  227.         }
  228.         if (! forward) {
  229.             h = -h;
  230.         }

  232.         // perform an Euler step using the preceding rough guess
  233.         final T[] y1 = MathArrays.buildArray(getField(), y0.length);
  234.         for (int j = 0; j < y0.length; ++j) {
  235.             y1[j] = y0[j].add(yDot0[j].multiply(h));
  236.         }
  237.         final T[] yDot1 = computeDerivatives(state0.getTime().add(h), y1);

  239.         // estimate the second derivative of the solution
  240.         double yDDotOnScale = 0;
  241.         for (int j = 0; j < scale.length; ++j) {
  242.             final double ratioDotDot = (yDot1[j].getReal() - yDot0[j].getReal()) / scale[j].getReal();
  243.             yDDotOnScale += ratioDotDot * ratioDotDot;
  244.         }
  245.         yDDotOnScale = FastMath.sqrt(yDDotOnScale) / h;

  247.         // step size is computed such that
  248.         // h^order * max (||y'/tol||, ||y''/tol||) = 0.01
  249.         final double maxInv2 = FastMath.max(FastMath.sqrt(yDotOnScale2), yDDotOnScale);
  250.         final double h1 = (maxInv2 < 1.0e-15) ?
  251.                           FastMath.max(1.0e-6, 0.001 * FastMath.abs(h)) :
  252.                           FastMath.pow(0.01 / maxInv2, 1.0 / order);
  253.         h = FastMath.min(100.0 * FastMath.abs(h), h1);
  254.         h = FastMath.max(h, 1.0e-12 * FastMath.abs(state0.getTime().getReal()));  // avoids cancellation when computing t1 - t0
  255.         if (h < getMinStep()) {
  256.             h = getMinStep();
  257.         }
  258.         if (h > getMaxStep()) {
  259.             h = getMaxStep();
  260.         }

  262.         if (! forward) {
  263.             h = -h;
  264.         }

  266.         return h;

  268.     }

  270.     /** Reset internal state to dummy values. */
  271.     protected void resetInternalState() {
  272.         setStepStart(null);
  273.         setStepSize(getField().getZero().add(stepsizeHelper.getDummyStepsize()));
  274.     }

  276.     /** Get the minimal step.
  277.      * @return minimal step
  278.      */
  279.     public double getMinStep() {
  280.         return stepsizeHelper.getMinStep();
  281.     }

  283.     /** Get the maximal step.
  284.      * @return maximal step
  285.      */
  286.     public double getMaxStep() {
  287.         return stepsizeHelper.getMaxStep();
  288.     }

  290. }
Created with JaCoCo 0.8.12.202403310830