StepsizeHelper.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.CalculusFieldElement;
  21. import org.hipparchus.exception.LocalizedCoreFormats;
  22. import org.hipparchus.exception.MathIllegalArgumentException;
  23. import org.hipparchus.ode.LocalizedODEFormats;
  24. import org.hipparchus.util.FastMath;

  26. /** Helper for adaptive stepsize control.
  27.  * @since 2.0
  28.  */

  30. public class StepsizeHelper {

  32.     /** Allowed absolute scalar error. */
  33.     private double scalAbsoluteTolerance;

  35.     /** Allowed relative scalar error. */
  36.     private double scalRelativeTolerance;

  38.     /** Allowed absolute vectorial error. */
  39.     private double[] vecAbsoluteTolerance;

  41.     /** Allowed relative vectorial error. */
  42.     private double[] vecRelativeTolerance;

  44.     /** Main set dimension. */
  45.     private int mainSetDimension;

  47.     /** User supplied initial step. */
  48.     private double initialStep;

  50.     /** Minimal step. */
  51.     private double minStep;

  53.     /** Maximal step. */
  54.     private double maxStep;

  56.     /** Simple constructor.
  57.      * @param minStep minimal step (sign is irrelevant, regardless of
  58.      * integration direction, forward or backward), the last step can
  59.      * be smaller than this
  60.      * @param maxStep maximal step (sign is irrelevant, regardless of
  61.      * integration direction, forward or backward), the last step can
  62.      * be smaller than this
  63.      * @param scalAbsoluteTolerance allowed absolute error
  64.      * @param scalRelativeTolerance allowed relative error
  65.      */
  66.     public StepsizeHelper(final double minStep, final double maxStep,
  67.                           final double scalAbsoluteTolerance,
  68.                           final double scalRelativeTolerance) {
  69.         this.minStep     = FastMath.abs(minStep);
  70.         this.maxStep     = FastMath.abs(maxStep);
  71.         this.initialStep = -1;

  73.         this.scalAbsoluteTolerance = scalAbsoluteTolerance;
  74.         this.scalRelativeTolerance = scalRelativeTolerance;
  75.         this.vecAbsoluteTolerance  = null;
  76.         this.vecRelativeTolerance  = null;
  77.     }

  79.     /** Simple constructor..
  80.      * @param minStep minimal step (sign is irrelevant, regardless of
  81.      * integration direction, forward or backward), the last step can
  82.      * be smaller than this
  83.      * @param maxStep maximal step (sign is irrelevant, regardless of
  84.      * integration direction, forward or backward), the last step can
  85.      * be smaller than this
  86.      * @param vecAbsoluteTolerance allowed absolute error
  87.      * @param vecRelativeTolerance allowed relative error
  88.      */
  89.     public StepsizeHelper(final double minStep, final double maxStep,
  90.                           final double[] vecAbsoluteTolerance,
  91.                           final double[] vecRelativeTolerance) {

  93.         this.minStep     = FastMath.abs(minStep);
  94.         this.maxStep     = FastMath.abs(maxStep);
  95.         this.initialStep = -1;

  97.        this.scalAbsoluteTolerance = 0;
  98.        this.scalRelativeTolerance = 0;
  99.        this.vecAbsoluteTolerance  = vecAbsoluteTolerance.clone();
  100.        this.vecRelativeTolerance  = vecRelativeTolerance.clone();

  102.     }

  104.     /** Set main set dimension.
  105.      * @param mainSetDimension dimension of the main set
  106.      * @exception MathIllegalArgumentException if adaptive step size integrators
  107.      * tolerance arrays dimensions are not compatible with equations settings
  108.      */
  109.     protected void setMainSetDimension(final int mainSetDimension) throws MathIllegalArgumentException {
  110.         this.mainSetDimension = mainSetDimension;

  112.         if (vecAbsoluteTolerance != null && vecAbsoluteTolerance.length != mainSetDimension) {
  113.             throw new MathIllegalArgumentException(LocalizedCoreFormats.DIMENSIONS_MISMATCH,
  114.                                                    mainSetDimension, vecAbsoluteTolerance.length);
  115.         }

  117.         if (vecRelativeTolerance != null && vecRelativeTolerance.length != mainSetDimension) {
  118.             throw new MathIllegalArgumentException(LocalizedCoreFormats.DIMENSIONS_MISMATCH,
  119.                                                    mainSetDimension, vecRelativeTolerance.length);
  120.         }
  121.     }

  123.     /** Get the main set dimension.
  124.      * @return main set dimension
  125.      */
  126.     public int getMainSetDimension() {
  127.         return mainSetDimension;
  128.     }

  130.     /** Get the relative tolerance for one component.
  131.      * @param i component to select
  132.      * @return relative tolerance for selected component
  133.      */
  134.     public double getRelativeTolerance(final int i) {
  135.         return vecAbsoluteTolerance == null ? scalRelativeTolerance : vecRelativeTolerance[i];
  136.     }

  138.     /** Get the tolerance for one component.
  139.      * @param i component to select
  140.      * @param scale scale factor for relative tolerance (i.e. y[i])
  141.      * @return tolerance for selected component
  142.      */
  143.     public double getTolerance(final int i, final double scale) {
  144.         return vecAbsoluteTolerance == null ?
  145.                scalAbsoluteTolerance   + scalRelativeTolerance   * scale :
  146.                vecAbsoluteTolerance[i] + vecRelativeTolerance[i] * scale;
  147.     }

  149.     /** Get the tolerance for one component.
  150.      * @param i component to select
  151.      * @param scale scale factor for relative tolerance (i.e. y[i])
  152.      * @param <T> type of the field elements
  153.      * @return tolerance for selected component
  154.      */
  155.     public <T extends CalculusFieldElement<T>> T getTolerance(final int i, final T scale) {
  156.         return vecAbsoluteTolerance == null ?
  157.                scale.multiply(scalRelativeTolerance).add(scalAbsoluteTolerance) :
  158.                scale.multiply(vecRelativeTolerance[i]).add(vecAbsoluteTolerance[i]);
  159.     }

  161.     /** Filter the integration step.
  162.      * @param h signed step
  163.      * @param forward forward integration indicator
  164.      * @param acceptSmall if true, steps smaller than the minimal value
  165.      * are silently increased up to this value, if false such small
  166.      * steps generate an exception
  167.      * @return a bounded integration step (h if no bound is reach, or a bounded value)
  168.      * @exception MathIllegalArgumentException if the step is too small and acceptSmall is false
  169.      */
  170.     public double filterStep(final double h, final boolean forward, final boolean acceptSmall)
  171.         throws MathIllegalArgumentException {

  173.         double filteredH = h;
  174.         if (FastMath.abs(h) < minStep) {
  175.             if (acceptSmall) {
  176.                 filteredH = forward ? minStep : -minStep;
  177.             } else {
  178.                 throw new MathIllegalArgumentException(LocalizedODEFormats.MINIMAL_STEPSIZE_REACHED_DURING_INTEGRATION,
  179.                                                        FastMath.abs(h), minStep, true);
  180.             }
  181.         }

  183.         if (filteredH > maxStep) {
  184.             filteredH = maxStep;
  185.         } else if (filteredH < -maxStep) {
  186.             filteredH = -maxStep;
  187.         }

  189.         return filteredH;

  191.     }

  193.     /** Filter the integration step.
  194.      * @param h signed step
  195.      * @param forward forward integration indicator
  196.      * @param acceptSmall if true, steps smaller than the minimal value
  197.      * are silently increased up to this value, if false such small
  198.      * steps generate an exception
  199.      * @param <T> type of the field elements
  200.      * @return a bounded integration step (h if no bound is reach, or a bounded value)
  201.      * @exception MathIllegalArgumentException if the step is too small and acceptSmall is false
  202.      */
  203.     public <T extends CalculusFieldElement<T>> T filterStep(final T h, final boolean forward, final boolean acceptSmall)
  204.         throws MathIllegalArgumentException {

  206.         T filteredH = h;
  207.         if (h.abs().subtract(minStep).getReal() < 0) {
  208.             if (acceptSmall) {
  209.                 filteredH = h.getField().getZero().add(forward ? minStep : -minStep);
  210.             } else {
  211.                 throw new MathIllegalArgumentException(LocalizedODEFormats.MINIMAL_STEPSIZE_REACHED_DURING_INTEGRATION,
  212.                                                        FastMath.abs(h.getReal()), minStep, true);
  213.             }
  214.         }

  216.         if (filteredH.subtract(maxStep).getReal() > 0) {
  217.             filteredH = h.getField().getZero().newInstance(maxStep);
  218.         } else if (filteredH.add(maxStep).getReal() < 0) {
  219.             filteredH = h.getField().getZero().newInstance(-maxStep);
  220.         }

  222.         return filteredH;

  224.     }

  226.     /** Set the initial step size.
  227.      * <p>This method allows the user to specify an initial positive
  228.      * step size instead of letting the integrator guess it by
  229.      * itself. If this method is not called before integration is
  230.      * started, the initial step size will be estimated by the
  231.      * integrator.</p>
  232.      * @param initialStepSize initial step size to use (must be positive even
  233.      * for backward integration ; providing a negative value or a value
  234.      * outside of the min/max step interval will lead the integrator to
  235.      * ignore the value and compute the initial step size by itself)
  236.      */
  237.     public void setInitialStepSize(final double initialStepSize) {
  238.         if ((initialStepSize < minStep) || (initialStepSize > maxStep)) {
  239.             initialStep = -1.0;
  240.         } else {
  241.             initialStep = initialStepSize;
  242.         }
  243.     }

  245.     /** Get the initial step.
  246.      * @return initial step
  247.      */
  248.     public double getInitialStep() {
  249.         return initialStep;
  250.     }

  252.     /** Get the minimal step.
  253.      * @return minimal step
  254.      */
  255.     public double getMinStep() {
  256.         return minStep;
  257.     }

  259.     /** Get the maximal step.
  260.      * @return maximal step
  261.      */
  262.     public double getMaxStep() {
  263.         return maxStep;
  264.     }

  266.     /** Get a dummy step size.
  267.      * @return geometric mean of {@link #getMinStep()} and {@link #getMaxStep()}
  268.      */
  269.     public double getDummyStepsize() {
  270.         return FastMath.sqrt(minStep * maxStep);
  271.     }

  273. }
Created with JaCoCo 0.8.12.202403310830