AdaptiveStepsizeIntegrator.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.ode.AbstractIntegrator;
  23. import org.hipparchus.ode.ODEState;
  24. import org.hipparchus.ode.ODEStateAndDerivative;
  25. import org.hipparchus.util.FastMath;

  27. /**
  28.  * This abstract class holds the common part of all adaptive
  29.  * stepsize integrators for Ordinary Differential Equations.
  30.  *
  31.  * <p>These algorithms perform integration with stepsize control, which
  32.  * means the user does not specify the integration step but rather a
  33.  * tolerance on error. The error threshold is computed as
  34.  * </p>
  35.  * <pre>
  36.  * threshold_i = absTol_i + relTol_i * max (abs (ym), abs (ym+1))
  37.  * </pre>
  38.  * <p>
  39.  * where absTol_i is the absolute tolerance for component i of the
  40.  * state vector and relTol_i is the relative tolerance for the same
  41.  * component. The user can also use only two scalar values absTol and
  42.  * relTol which will be used for all components.
  43.  * </p>
  44.  * <p>
  45.  * If the Ordinary Differential Equations is an {@link org.hipparchus.ode.ExpandableODE
  46.  * extended ODE} rather than a {@link
  47.  * org.hipparchus.ode.OrdinaryDifferentialEquation basic ODE}, then
  48.  * <em>only</em> the {@link org.hipparchus.ode.ExpandableODE#getPrimary() primary part}
  49.  * of the state vector is used for stepsize control, not the complete state vector.
  50.  * </p>
  51.  *
  52.  * <p>If the estimated error for ym+1 is such that</p>
  53.  * <pre>
  54.  * sqrt((sum (errEst_i / threshold_i)^2 ) / n) &lt; 1
  55.  * </pre>
  56.  *
  57.  * <p>(where n is the main set dimension) then the step is accepted,
  58.  * otherwise the step is rejected and a new attempt is made with a new
  59.  * stepsize.</p>
  60.  *
  61.  *
  62.  */

  64. public abstract class AdaptiveStepsizeIntegrator
  65.     extends AbstractIntegrator {

  67.     /** Helper for step size control. */
  68.     private StepsizeHelper stepsizeHelper;

  70.     /** Build an integrator with the given stepsize bounds.
  71.      * The default step handler does nothing.
  72.      * @param name name of the method
  73.      * @param minStep minimal step (sign is irrelevant, regardless of
  74.      * integration direction, forward or backward), the last step can
  75.      * be smaller than this
  76.      * @param maxStep maximal step (sign is irrelevant, regardless of
  77.      * integration direction, forward or backward), the last step can
  78.      * be smaller than this
  79.      * @param scalAbsoluteTolerance allowed absolute error
  80.      * @param scalRelativeTolerance allowed relative error
  81.      */
  82.     protected AdaptiveStepsizeIntegrator(final String name,
  83.                                       final double minStep, final double maxStep,
  84.                                       final double scalAbsoluteTolerance,
  85.                                       final double scalRelativeTolerance) {
  86.         super(name);
  87.         stepsizeHelper = new StepsizeHelper(minStep, maxStep, scalAbsoluteTolerance, scalRelativeTolerance);
  88.         resetInternalState();
  89.     }

  91.     /** Build an integrator with the given stepsize bounds.
  92.      * The default step handler does nothing.
  93.      * @param name name of the method
  94.      * @param minStep minimal step (sign is irrelevant, regardless of
  95.      * integration direction, forward or backward), the last step can
  96.      * be smaller than this
  97.      * @param maxStep maximal step (sign is irrelevant, regardless of
  98.      * integration direction, forward or backward), the last step can
  99.      * be smaller than this
  100.      * @param vecAbsoluteTolerance allowed absolute error
  101.      * @param vecRelativeTolerance allowed relative error
  102.      */
  103.     protected AdaptiveStepsizeIntegrator(final String name,
  104.                                       final double minStep, final double maxStep,
  105.                                       final double[] vecAbsoluteTolerance,
  106.                                       final double[] vecRelativeTolerance) {
  107.         super(name);
  108.         stepsizeHelper = new StepsizeHelper(minStep, maxStep, vecAbsoluteTolerance, vecRelativeTolerance);
  109.         resetInternalState();
  110.     }

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

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

  152.     /** Get the stepsize helper.
  153.      * @return stepsize helper
  154.      * @since 2.0
  155.      */
  156.     protected StepsizeHelper getStepSizeHelper() {
  157.         return stepsizeHelper;
  158.     }

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

  175.     /** {@inheritDoc} */
  176.     @Override
  177.     protected void sanityChecks(final ODEState initialState, final double t)
  178.                     throws MathIllegalArgumentException {
  179.         super.sanityChecks(initialState, t);
  180.         stepsizeHelper.setMainSetDimension(initialState.getPrimaryStateDimension());
  181.     }

  183.     /** Initialize the integration step.
  184.      * @param forward forward integration indicator
  185.      * @param order order of the method
  186.      * @param scale scaling vector for the state vector (can be shorter than state vector)
  187.      * @param state0 state at integration start time
  188.      * @return first integration step
  189.      * @exception MathIllegalStateException if the number of functions evaluations is exceeded
  190.      * @exception MathIllegalArgumentException if arrays dimensions do not match equations settings
  191.      */
  192.     public double initializeStep(final boolean forward, final int order, final double[] scale,
  193.                                  final ODEStateAndDerivative state0)
  194.         throws MathIllegalArgumentException, MathIllegalStateException {

  196.         if (stepsizeHelper.getInitialStep() > 0) {
  197.             // use the user provided value
  198.             return forward ? stepsizeHelper.getInitialStep() : -stepsizeHelper.getInitialStep();
  199.         }

  201.         // very rough first guess : h = 0.01 * ||y/scale|| / ||y'/scale||
  202.         // this guess will be used to perform an Euler step
  203.         final double[] y0    = state0.getCompleteState();
  204.         final double[] yDot0 = state0.getCompleteDerivative();
  205.         double yOnScale2 = 0;
  206.         double yDotOnScale2 = 0;
  207.         for (int j = 0; j < scale.length; ++j) {
  208.             final double ratio    = y0[j] / scale[j];
  209.             yOnScale2            += ratio * ratio;
  210.             final double ratioDot = yDot0[j] / scale[j];
  211.             yDotOnScale2         += ratioDot * ratioDot;
  212.         }

  214.         double h = ((yOnScale2 < 1.0e-10) || (yDotOnScale2 < 1.0e-10)) ?
  215.                    1.0e-6 : (0.01 * FastMath.sqrt(yOnScale2 / yDotOnScale2));
  216.         if (h > getMaxStep()) {
  217.             h = getMaxStep();
  218.         }
  219.         if (! forward) {
  220.             h = -h;
  221.         }

  223.         // perform an Euler step using the preceding rough guess
  224.         final double[] y1 = new double[y0.length];
  225.         for (int j = 0; j < y0.length; ++j) {
  226.             y1[j] = y0[j] + h * yDot0[j];
  227.         }
  228.         final double[] yDot1 = computeDerivatives(state0.getTime() + h, y1);

  230.         // estimate the second derivative of the solution
  231.         double yDDotOnScale = 0;
  232.         for (int j = 0; j < scale.length; ++j) {
  233.             final double ratioDotDot = (yDot1[j] - yDot0[j]) / scale[j];
  234.             yDDotOnScale += ratioDotDot * ratioDotDot;
  235.         }
  236.         yDDotOnScale = FastMath.sqrt(yDDotOnScale) / h;

  238.         // step size is computed such that
  239.         // h^order * max (||y'/tol||, ||y''/tol||) = 0.01
  240.         final double maxInv2 = FastMath.max(FastMath.sqrt(yDotOnScale2), yDDotOnScale);
  241.         final double h1 = (maxInv2 < 1.0e-15) ?
  242.                            FastMath.max(1.0e-6, 0.001 * FastMath.abs(h)) :
  243.                            FastMath.pow(0.01 / maxInv2, 1.0 / order);
  244.         h = FastMath.min(100.0 * FastMath.abs(h), h1);
  245.         h = FastMath.max(h, 1.0e-12 * FastMath.abs(state0.getTime()));  // avoids cancellation when computing t1 - t0
  246.         if (h < getMinStep()) {
  247.             h = getMinStep();
  248.         }
  249.         if (h > getMaxStep()) {
  250.             h = getMaxStep();
  251.         }
  252.         if (! forward) {
  253.             h = -h;
  254.         }

  256.         return h;

  258.     }

  260.     /** Reset internal state to dummy values. */
  261.     protected void resetInternalState() {
  262.         setStepStart(null);
  263.         setStepSize(stepsizeHelper.getDummyStepsize());
  264.     }

  266.     /** Get the minimal step.
  267.      * @return minimal step
  268.      */
  269.     public double getMinStep() {
  270.         return stepsizeHelper.getMinStep();
  271.     }

  273.     /** Get the maximal step.
  274.      * @return maximal step
  275.      */
  276.     public double getMaxStep() {
  277.         return stepsizeHelper.getMaxStep();
  278.     }

  280. }
Created with JaCoCo 0.8.12.202403310830