AdaptiveStepsizeIntegrator.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.  */

  17. package org.apache.commons.math4.legacy.ode.nonstiff;

  18. import org.apache.commons.math4.legacy.exception.DimensionMismatchException;
  19. import org.apache.commons.math4.legacy.exception.MaxCountExceededException;
  20. import org.apache.commons.math4.legacy.exception.NoBracketingException;
  21. import org.apache.commons.math4.legacy.exception.NumberIsTooSmallException;
  22. import org.apache.commons.math4.legacy.exception.util.LocalizedFormats;
  23. import org.apache.commons.math4.legacy.ode.AbstractIntegrator;
  24. import org.apache.commons.math4.legacy.ode.ExpandableStatefulODE;
  25. import org.apache.commons.math4.core.jdkmath.JdkMath;

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

  61. public abstract class AdaptiveStepsizeIntegrator
  62.   extends AbstractIntegrator {

  63.     /** Allowed absolute scalar error. */
  64.     protected double scalAbsoluteTolerance;

  65.     /** Allowed relative scalar error. */
  66.     protected double scalRelativeTolerance;

  67.     /** Allowed absolute vectorial error. */
  68.     protected double[] vecAbsoluteTolerance;

  69.     /** Allowed relative vectorial error. */
  70.     protected double[] vecRelativeTolerance;

  71.     /** Main set dimension. */
  72.     protected int mainSetDimension;

  73.     /** User supplied initial step. */
  74.     private double initialStep;

  75.     /** Minimal step. */
  76.     private double minStep;

  77.     /** Maximal step. */
  78.     private double maxStep;

  79.   /** Build an integrator with the given stepsize bounds.
  80.    * The default step handler does nothing.
  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.   public AdaptiveStepsizeIntegrator(final String name,
  92.                                     final double minStep, final double maxStep,
  93.                                     final double scalAbsoluteTolerance,
  94.                                     final double scalRelativeTolerance) {

  95.     super(name);
  96.     setStepSizeControl(minStep, maxStep, scalAbsoluteTolerance, scalRelativeTolerance);
  97.     resetInternalState();
  98.   }

  99.   /** Build an integrator with the given stepsize bounds.
  100.    * The default step handler does nothing.
  101.    * @param name name of the method
  102.    * @param minStep minimal step (sign is irrelevant, regardless of
  103.    * integration direction, forward or backward), the last step can
  104.    * be smaller than this
  105.    * @param maxStep maximal step (sign is irrelevant, regardless of
  106.    * integration direction, forward or backward), the last step can
  107.    * be smaller than this
  108.    * @param vecAbsoluteTolerance allowed absolute error
  109.    * @param vecRelativeTolerance allowed relative error
  110.    */
  111.   public AdaptiveStepsizeIntegrator(final String name,
  112.                                     final double minStep, final double maxStep,
  113.                                     final double[] vecAbsoluteTolerance,
  114.                                     final double[] vecRelativeTolerance) {

  115.     super(name);
  116.     setStepSizeControl(minStep, maxStep, vecAbsoluteTolerance, vecRelativeTolerance);
  117.     resetInternalState();
  118.   }

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

  136.       minStep     = JdkMath.abs(minimalStep);
  137.       maxStep     = JdkMath.abs(maximalStep);
  138.       initialStep = -1;

  139.       scalAbsoluteTolerance = absoluteTolerance;
  140.       scalRelativeTolerance = relativeTolerance;
  141.       vecAbsoluteTolerance  = null;
  142.       vecRelativeTolerance  = null;
  143.   }

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

  161.       minStep     = JdkMath.abs(minimalStep);
  162.       maxStep     = JdkMath.abs(maximalStep);
  163.       initialStep = -1;

  164.       scalAbsoluteTolerance = 0;
  165.       scalRelativeTolerance = 0;
  166.       vecAbsoluteTolerance  = absoluteTolerance.clone();
  167.       vecRelativeTolerance  = relativeTolerance.clone();
  168.   }

  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.     if (initialStepSize < minStep || initialStepSize > maxStep) {
  182.       initialStep = -1.0;
  183.     } else {
  184.       initialStep = initialStepSize;
  185.     }
  186.   }

  187.   /** {@inheritDoc} */
  188.   @Override
  189.   protected void sanityChecks(final ExpandableStatefulODE equations, final double t)
  190.       throws DimensionMismatchException, NumberIsTooSmallException {

  191.       super.sanityChecks(equations, t);

  192.       mainSetDimension = equations.getPrimaryMapper().getDimension();

  193.       if (vecAbsoluteTolerance != null && vecAbsoluteTolerance.length != mainSetDimension) {
  194.           throw new DimensionMismatchException(mainSetDimension, vecAbsoluteTolerance.length);
  195.       }

  196.       if (vecRelativeTolerance != null && vecRelativeTolerance.length != mainSetDimension) {
  197.           throw new DimensionMismatchException(mainSetDimension, vecRelativeTolerance.length);
  198.       }
  199.   }

  200.   /** Initialize the integration step.
  201.    * @param forward forward integration indicator
  202.    * @param order order of the method
  203.    * @param scale scaling vector for the state vector (can be shorter than state vector)
  204.    * @param t0 start time
  205.    * @param y0 state vector at t0
  206.    * @param yDot0 first time derivative of y0
  207.    * @param y1 work array for a state vector
  208.    * @param yDot1 work array for the first time derivative of y1
  209.    * @return first integration step
  210.    * @exception MaxCountExceededException if the number of functions evaluations is exceeded
  211.    * @exception DimensionMismatchException if arrays dimensions do not match equations settings
  212.    */
  213.   public double initializeStep(final boolean forward, final int order, final double[] scale,
  214.                                final double t0, final double[] y0, final double[] yDot0,
  215.                                final double[] y1, final double[] yDot1)
  216.       throws MaxCountExceededException, DimensionMismatchException {

  217.     if (initialStep > 0) {
  218.       // use the user provided value
  219.       return forward ? initialStep : -initialStep;
  220.     }

  221.     // very rough first guess : h = 0.01 * ||y/scale|| / ||y'/scale||
  222.     // this guess will be used to perform an Euler step
  223.     double ratio;
  224.     double yOnScale2 = 0;
  225.     double yDotOnScale2 = 0;
  226.     for (int j = 0; j < scale.length; ++j) {
  227.       ratio         = y0[j] / scale[j];
  228.       yOnScale2    += ratio * ratio;
  229.       ratio         = yDot0[j] / scale[j];
  230.       yDotOnScale2 += ratio * ratio;
  231.     }

  232.     double h = (yOnScale2 < 1.0e-10 || yDotOnScale2 < 1.0e-10) ?
  233.                1.0e-6 : (0.01 * JdkMath.sqrt(yOnScale2 / yDotOnScale2));
  234.     if (! forward) {
  235.       h = -h;
  236.     }

  237.     // perform an Euler step using the preceding rough guess
  238.     for (int j = 0; j < y0.length; ++j) {
  239.       y1[j] = y0[j] + h * yDot0[j];
  240.     }
  241.     computeDerivatives(t0 + h, y1, yDot1);

  242.     // estimate the second derivative of the solution
  243.     double yDDotOnScale = 0;
  244.     for (int j = 0; j < scale.length; ++j) {
  245.       ratio         = (yDot1[j] - yDot0[j]) / scale[j];
  246.       yDDotOnScale += ratio * ratio;
  247.     }
  248.     yDDotOnScale = JdkMath.sqrt(yDDotOnScale) / h;

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

  266.     return h;
  267.   }

  268.   /** Filter the integration step.
  269.    * @param h signed step
  270.    * @param forward forward integration indicator
  271.    * @param acceptSmall if true, steps smaller than the minimal value
  272.    * are silently increased up to this value, if false such small
  273.    * steps generate an exception
  274.    * @return a bounded integration step (h if no bound is reach, or a bounded value)
  275.    * @exception NumberIsTooSmallException if the step is too small and acceptSmall is false
  276.    */
  277.   protected double filterStep(final double h, final boolean forward, final boolean acceptSmall)
  278.     throws NumberIsTooSmallException {

  279.       double filteredH = h;
  280.       if (JdkMath.abs(h) < minStep) {
  281.           if (acceptSmall) {
  282.               filteredH = forward ? minStep : -minStep;
  283.           } else {
  284.               throw new NumberIsTooSmallException(LocalizedFormats.MINIMAL_STEPSIZE_REACHED_DURING_INTEGRATION,
  285.                                                   JdkMath.abs(h), minStep, true);
  286.           }
  287.       }

  288.       if (filteredH > maxStep) {
  289.           filteredH = maxStep;
  290.       } else if (filteredH < -maxStep) {
  291.           filteredH = -maxStep;
  292.       }

  293.       return filteredH;
  294.   }

  295.   /** {@inheritDoc} */
  296.   @Override
  297.   public abstract void integrate (ExpandableStatefulODE equations, double t)
  298.       throws NumberIsTooSmallException, DimensionMismatchException,
  299.              MaxCountExceededException, NoBracketingException;

  300.   /** {@inheritDoc} */
  301.   @Override
  302.   public double getCurrentStepStart() {
  303.     return stepStart;
  304.   }

  305.   /** Reset internal state to dummy values. */
  306.   protected void resetInternalState() {
  307.     stepStart = Double.NaN;
  308.     stepSize  = JdkMath.sqrt(minStep * maxStep);
  309.   }

  310.   /** Get the minimal step.
  311.    * @return minimal step
  312.    */
  313.   public double getMinStep() {
  314.     return minStep;
  315.   }

  316.   /** Get the maximal step.
  317.    * @return maximal step
  318.    */
  319.   public double getMaxStep() {
  320.     return maxStep;
  321.   }
  322. }