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.  *      http://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.apache.commons.math4.legacy.ode.nonstiff;

  20. import org.apache.commons.math4.legacy.core.Field;
  21. import org.apache.commons.math4.legacy.core.RealFieldElement;
  22. import org.apache.commons.math4.legacy.exception.DimensionMismatchException;
  23. import org.apache.commons.math4.legacy.exception.MaxCountExceededException;
  24. import org.apache.commons.math4.legacy.exception.NumberIsTooSmallException;
  25. import org.apache.commons.math4.legacy.exception.util.LocalizedFormats;
  26. import org.apache.commons.math4.legacy.ode.AbstractFieldIntegrator;
  27. import org.apache.commons.math4.legacy.ode.FieldEquationsMapper;
  28. import org.apache.commons.math4.legacy.ode.FieldODEState;
  29. import org.apache.commons.math4.legacy.ode.FieldODEStateAndDerivative;
  30. import org.apache.commons.math4.core.jdkmath.JdkMath;
  31. import org.apache.commons.math4.legacy.core.MathArrays;

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

  69. public abstract class AdaptiveStepsizeFieldIntegrator<T extends RealFieldElement<T>>
  70.     extends AbstractFieldIntegrator<T> {

  72.     /** Allowed absolute scalar error. */
  73.     protected double scalAbsoluteTolerance;

  75.     /** Allowed relative scalar error. */
  76.     protected double scalRelativeTolerance;

  78.     /** Allowed absolute vectorial error. */
  79.     protected double[] vecAbsoluteTolerance;

  81.     /** Allowed relative vectorial error. */
  82.     protected double[] vecRelativeTolerance;

  84.     /** Main set dimension. */
  85.     protected int mainSetDimension;

  87.     /** User supplied initial step. */
  88.     private T initialStep;

  90.     /** Minimal step. */
  91.     private T minStep;

  93.     /** Maximal step. */
  94.     private T maxStep;

  96.     /** Build an integrator with the given stepsize bounds.
  97.      * The default step handler does nothing.
  98.      * @param field field to which the time and state vector elements belong
  99.      * @param name name of the method
  100.      * @param minStep minimal step (sign is irrelevant, regardless of
  101.      * integration direction, forward or backward), the last step can
  102.      * be smaller than this
  103.      * @param maxStep maximal step (sign is irrelevant, regardless of
  104.      * integration direction, forward or backward), the last step can
  105.      * be smaller than this
  106.      * @param scalAbsoluteTolerance allowed absolute error
  107.      * @param scalRelativeTolerance allowed relative error
  108.      */
  109.     public AdaptiveStepsizeFieldIntegrator(final Field<T> field, final String name,
  110.                                            final double minStep, final double maxStep,
  111.                                            final double scalAbsoluteTolerance,
  112.                                            final double scalRelativeTolerance) {

  114.         super(field, name);
  115.         setStepSizeControl(minStep, maxStep, scalAbsoluteTolerance, scalRelativeTolerance);
  116.         resetInternalState();
  117.     }

  119.     /** Build an integrator with the given stepsize bounds.
  120.      * The default step handler does nothing.
  121.      * @param field field to which the time and state vector elements belong
  122.      * @param name name of the method
  123.      * @param minStep minimal step (sign is irrelevant, regardless of
  124.      * integration direction, forward or backward), the last step can
  125.      * be smaller than this
  126.      * @param maxStep maximal step (sign is irrelevant, regardless of
  127.      * integration direction, forward or backward), the last step can
  128.      * be smaller than this
  129.      * @param vecAbsoluteTolerance allowed absolute error
  130.      * @param vecRelativeTolerance allowed relative error
  131.      */
  132.     public AdaptiveStepsizeFieldIntegrator(final Field<T> field, final String name,
  133.                                            final double minStep, final double maxStep,
  134.                                            final double[] vecAbsoluteTolerance,
  135.                                            final double[] vecRelativeTolerance) {

  137.         super(field, name);
  138.         setStepSizeControl(minStep, maxStep, vecAbsoluteTolerance, vecRelativeTolerance);
  139.         resetInternalState();
  140.     }

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

  160.         minStep     = getField().getZero().add(JdkMath.abs(minimalStep));
  161.         maxStep     = getField().getZero().add(JdkMath.abs(maximalStep));
  162.         initialStep = getField().getOne().negate();

  164.         scalAbsoluteTolerance = absoluteTolerance;
  165.         scalRelativeTolerance = relativeTolerance;
  166.         vecAbsoluteTolerance  = null;
  167.         vecRelativeTolerance  = null;
  168.     }

  170.     /** Set the adaptive step size control parameters.
  171.      * <p>
  172.      * A side effect of this method is to also reset the initial
  173.      * step so it will be automatically computed by the integrator
  174.      * if {@link #setInitialStepSize(RealFieldElement) setInitialStepSize}
  175.      * is not called by the user.
  176.      * </p>
  177.      * @param minimalStep minimal step (must be positive even for backward
  178.      * integration), the last step can be smaller than this
  179.      * @param maximalStep maximal step (must be positive even for backward
  180.      * integration)
  181.      * @param absoluteTolerance allowed absolute error
  182.      * @param relativeTolerance allowed relative error
  183.      */
  184.     public void setStepSizeControl(final double minimalStep, final double maximalStep,
  185.                                    final double[] absoluteTolerance,
  186.                                    final double[] relativeTolerance) {

  188.         minStep     = getField().getZero().add(JdkMath.abs(minimalStep));
  189.         maxStep     = getField().getZero().add(JdkMath.abs(maximalStep));
  190.         initialStep = getField().getOne().negate();

  192.         scalAbsoluteTolerance = 0;
  193.         scalRelativeTolerance = 0;
  194.         vecAbsoluteTolerance  = absoluteTolerance.clone();
  195.         vecRelativeTolerance  = relativeTolerance.clone();
  196.     }

  198.     /** Set the initial step size.
  199.      * <p>This method allows the user to specify an initial positive
  200.      * step size instead of letting the integrator guess it by
  201.      * itself. If this method is not called before integration is
  202.      * started, the initial step size will be estimated by the
  203.      * integrator.</p>
  204.      * @param initialStepSize initial step size to use (must be positive even
  205.      * for backward integration ; providing a negative value or a value
  206.      * outside of the min/max step interval will lead the integrator to
  207.      * ignore the value and compute the initial step size by itself)
  208.      */
  209.     public void setInitialStepSize(final T initialStepSize) {
  210.         if (initialStepSize.subtract(minStep).getReal() < 0 ||
  211.             initialStepSize.subtract(maxStep).getReal() > 0) {
  212.             initialStep = getField().getOne().negate();
  213.         } else {
  214.             initialStep = initialStepSize;
  215.         }
  216.     }

  218.     /** {@inheritDoc} */
  219.     @Override
  220.     protected void sanityChecks(final FieldODEState<T> eqn, final T t)
  221.         throws DimensionMismatchException, NumberIsTooSmallException {

  223.         super.sanityChecks(eqn, t);

  225.         mainSetDimension = eqn.getStateDimension();

  227.         if (vecAbsoluteTolerance != null && vecAbsoluteTolerance.length != mainSetDimension) {
  228.             throw new DimensionMismatchException(mainSetDimension, vecAbsoluteTolerance.length);
  229.         }

  231.         if (vecRelativeTolerance != null && vecRelativeTolerance.length != mainSetDimension) {
  232.             throw new DimensionMismatchException(mainSetDimension, vecRelativeTolerance.length);
  233.         }
  234.     }

  236.     /** Initialize the integration step.
  237.      * @param forward forward integration indicator
  238.      * @param order order of the method
  239.      * @param scale scaling vector for the state vector (can be shorter than state vector)
  240.      * @param state0 state at integration start time
  241.      * @param mapper mapper for all the equations
  242.      * @return first integration step
  243.      * @exception MaxCountExceededException if the number of functions evaluations is exceeded
  244.      * @exception DimensionMismatchException if arrays dimensions do not match equations settings
  245.      */
  246.     public T initializeStep(final boolean forward, final int order, final T[] scale,
  247.                             final FieldODEStateAndDerivative<T> state0,
  248.                             final FieldEquationsMapper<T> mapper)
  249.         throws MaxCountExceededException, DimensionMismatchException {

  251.         if (initialStep.getReal() > 0) {
  252.             // use the user provided value
  253.             return forward ? initialStep : initialStep.negate();
  254.         }

  256.         // very rough first guess : h = 0.01 * ||y/scale|| / ||y'/scale||
  257.         // this guess will be used to perform an Euler step
  258.         final T[] y0    = mapper.mapState(state0);
  259.         final T[] yDot0 = mapper.mapDerivative(state0);
  260.         T yOnScale2    = getField().getZero();
  261.         T yDotOnScale2 = getField().getZero();
  262.         for (int j = 0; j < scale.length; ++j) {
  263.             final T ratio    = y0[j].divide(scale[j]);
  264.             yOnScale2        = yOnScale2.add(ratio.multiply(ratio));
  265.             final T ratioDot = yDot0[j].divide(scale[j]);
  266.             yDotOnScale2     = yDotOnScale2.add(ratioDot.multiply(ratioDot));
  267.         }

  269.         T h = (yOnScale2.getReal() < 1.0e-10 || yDotOnScale2.getReal() < 1.0e-10) ?
  270.               getField().getZero().add(1.0e-6) :
  271.               yOnScale2.divide(yDotOnScale2).sqrt().multiply(0.01);
  272.         if (! forward) {
  273.             h = h.negate();
  274.         }

  276.         // perform an Euler step using the preceding rough guess
  277.         final T[] y1 = MathArrays.buildArray(getField(), y0.length);
  278.         for (int j = 0; j < y0.length; ++j) {
  279.             y1[j] = y0[j].add(yDot0[j].multiply(h));
  280.         }
  281.         final T[] yDot1 = computeDerivatives(state0.getTime().add(h), y1);

  283.         // estimate the second derivative of the solution
  284.         T yDDotOnScale = getField().getZero();
  285.         for (int j = 0; j < scale.length; ++j) {
  286.             final T ratioDotDot = yDot1[j].subtract(yDot0[j]).divide(scale[j]);
  287.             yDDotOnScale = yDDotOnScale.add(ratioDotDot.multiply(ratioDotDot));
  288.         }
  289.         yDDotOnScale = yDDotOnScale.sqrt().divide(h);

  291.         // step size is computed such that
  292.         // h^order * max (||y'/tol||, ||y''/tol||) = 0.01
  293.         final T maxInv2 = RealFieldElement.max(yDotOnScale2.sqrt(), yDDotOnScale);
  294.         final T h1 = maxInv2.getReal() < 1.0e-15 ?
  295.                      RealFieldElement.max(getField().getZero().add(1.0e-6), h.abs().multiply(0.001)) :
  296.                      maxInv2.multiply(100).reciprocal().pow(1.0 / order);
  297.         h = RealFieldElement.min(h.abs().multiply(100), h1);
  298.         h = RealFieldElement.max(h, state0.getTime().abs().multiply(1.0e-12));  // avoids cancellation when computing t1 - t0
  299.         h = RealFieldElement.max(minStep, RealFieldElement.min(maxStep, h));
  300.         if (! forward) {
  301.             h = h.negate();
  302.         }

  304.         return h;
  305.     }

  307.     /** Filter the integration step.
  308.      * @param h signed step
  309.      * @param forward forward integration indicator
  310.      * @param acceptSmall if true, steps smaller than the minimal value
  311.      * are silently increased up to this value, if false such small
  312.      * steps generate an exception
  313.      * @return a bounded integration step (h if no bound is reach, or a bounded value)
  314.      * @exception NumberIsTooSmallException if the step is too small and acceptSmall is false
  315.      */
  316.     protected T filterStep(final T h, final boolean forward, final boolean acceptSmall)
  317.         throws NumberIsTooSmallException {

  319.         T filteredH = h;
  320.         if (h.abs().subtract(minStep).getReal() < 0) {
  321.             if (acceptSmall) {
  322.                 filteredH = forward ? minStep : minStep.negate();
  323.             } else {
  324.                 throw new NumberIsTooSmallException(LocalizedFormats.MINIMAL_STEPSIZE_REACHED_DURING_INTEGRATION,
  325.                                                     h.abs().getReal(), minStep.getReal(), true);
  326.             }
  327.         }

  329.         if (filteredH.subtract(maxStep).getReal() > 0) {
  330.             filteredH = maxStep;
  331.         } else if (filteredH.add(maxStep).getReal() < 0) {
  332.             filteredH = maxStep.negate();
  333.         }

  335.         return filteredH;
  336.     }

  338.     /** Reset internal state to dummy values. */
  339.     protected void resetInternalState() {
  340.         setStepStart(null);
  341.         setStepSize(minStep.multiply(maxStep).sqrt());
  342.     }

  344.     /** Get the minimal step.
  345.      * @return minimal step
  346.      */
  347.     public T getMinStep() {
  348.         return minStep;
  349.     }

  351.     /** Get the maximal step.
  352.      * @return maximal step
  353.      */
  354.     public T getMaxStep() {
  355.         return maxStep;
  356.     }
  357. }
Created with JaCoCo 0.8.10.202304240956Code Coverage Report for Apache Commons Math 4.0-SNAPSHOT