    001    /*
    002     * Licensed to the Apache Software Foundation (ASF) under one or more
    003     * contributor license agreements.  See the NOTICE file distributed with
    004     * this work for additional information regarding copyright ownership.
    005     * The ASF licenses this file to You under the Apache License, Version 2.0
    006     * (the "License"); you may not use this file except in compliance with
    007     * the License.  You may obtain a copy of the License at
    008     *
    009     *      http://www.apache.org/licenses/LICENSE-2.0
    010     *
    011     * Unless required by applicable law or agreed to in writing, software
    012     * distributed under the License is distributed on an "AS IS" BASIS,
    013     * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
    014     * See the License for the specific language governing permissions and
    015     * limitations under the License.
    016     */
    017
    018    package org.apache.commons.math3.ode.nonstiff;
    019
    020    import org.apache.commons.math3.exception.DimensionMismatchException;
    021    import org.apache.commons.math3.exception.MaxCountExceededException;
    022    import org.apache.commons.math3.exception.NoBracketingException;
    023    import org.apache.commons.math3.exception.NumberIsTooSmallException;
    024    import org.apache.commons.math3.exception.util.LocalizedFormats;
    025    import org.apache.commons.math3.ode.AbstractIntegrator;
    026    import org.apache.commons.math3.ode.ExpandableStatefulODE;
    027    import org.apache.commons.math3.util.FastMath;
    028
    029    /**
    030     * This abstract class holds the common part of all adaptive
    031     * stepsize integrators for Ordinary Differential Equations.
    032     *
    033     * <p>These algorithms perform integration with stepsize control, which
    034     * means the user does not specify the integration step but rather a
    035     * tolerance on error. The error threshold is computed as
    036     * <pre>
    037     * threshold_i = absTol_i + relTol_i * max (abs (ym), abs (ym+1))
    038     * </pre>
    039     * where absTol_i is the absolute tolerance for component i of the
    040     * state vector and relTol_i is the relative tolerance for the same
    041     * component. The user can also use only two scalar values absTol and
    042     * relTol which will be used for all components.
    043     * </p>
    044     * <p>
    045     * If the Ordinary Differential Equations is an {@link ExpandableStatefulODE
    046     * extended ODE} rather than a {@link
    047     * org.apache.commons.math3.ode.FirstOrderDifferentialEquations basic ODE}, then
    048     * <em>only</em> the {@link ExpandableStatefulODE#getPrimaryState() primary part}
    049     * of the state vector is used for stepsize control, not the complete state vector.
    050     * </p>
    051     *
    052     * <p>If the estimated error for ym+1 is such that
    053     * <pre>
    054     * sqrt((sum (errEst_i / threshold_i)^2 ) / n) < 1
    055     * </pre>
    056     *
    057     * (where n is the main set dimension) then the step is accepted,
    058     * otherwise the step is rejected and a new attempt is made with a new
    059     * stepsize.</p>
    060     *
    061     * @version $Id: AdaptiveStepsizeIntegrator.java 1416643 2012-12-03 19:37:14Z tn $
    062     * @since 1.2
    063     *
    064     */
    065
    066    public abstract class AdaptiveStepsizeIntegrator
    067      extends AbstractIntegrator {
    068
    069        /** Allowed absolute scalar error. */
    070        protected double scalAbsoluteTolerance;
    071
    072        /** Allowed relative scalar error. */
    073        protected double scalRelativeTolerance;
    074
    075        /** Allowed absolute vectorial error. */
    076        protected double[] vecAbsoluteTolerance;
    077
    078        /** Allowed relative vectorial error. */
    079        protected double[] vecRelativeTolerance;
    080
    081        /** Main set dimension. */
    082        protected int mainSetDimension;
    083
    084        /** User supplied initial step. */
    085        private double initialStep;
    086
    087        /** Minimal step. */
    088        private double minStep;
    089
    090        /** Maximal step. */
    091        private double maxStep;
    092
    093      /** Build an integrator with the given stepsize bounds.
    094       * The default step handler does nothing.
    095       * @param name name of the method
    096       * @param minStep minimal step (sign is irrelevant, regardless of
    097       * integration direction, forward or backward), the last step can
    098       * be smaller than this
    099       * @param maxStep maximal step (sign is irrelevant, regardless of
    100       * integration direction, forward or backward), the last step can
    101       * be smaller than this
    102       * @param scalAbsoluteTolerance allowed absolute error
    103       * @param scalRelativeTolerance allowed relative error
    104       */
    105      public AdaptiveStepsizeIntegrator(final String name,
    106                                        final double minStep, final double maxStep,
    107                                        final double scalAbsoluteTolerance,
    108                                        final double scalRelativeTolerance) {
    109
    110        super(name);
    111        setStepSizeControl(minStep, maxStep, scalAbsoluteTolerance, scalRelativeTolerance);
    112        resetInternalState();
    113
    114      }
    115
    116      /** Build an integrator with the given stepsize bounds.
    117       * The default step handler does nothing.
    118       * @param name name of the method
    119       * @param minStep minimal step (sign is irrelevant, regardless of
    120       * integration direction, forward or backward), the last step can
    121       * be smaller than this
    122       * @param maxStep maximal step (sign is irrelevant, regardless of
    123       * integration direction, forward or backward), the last step can
    124       * be smaller than this
    125       * @param vecAbsoluteTolerance allowed absolute error
    126       * @param vecRelativeTolerance allowed relative error
    127       */
    128      public AdaptiveStepsizeIntegrator(final String name,
    129                                        final double minStep, final double maxStep,
    130                                        final double[] vecAbsoluteTolerance,
    131                                        final double[] vecRelativeTolerance) {
    132
    133        super(name);
    134        setStepSizeControl(minStep, maxStep, vecAbsoluteTolerance, vecRelativeTolerance);
    135        resetInternalState();
    136
    137      }
    138
    139      /** Set the adaptive step size control parameters.
    140       * <p>
    141       * A side effect of this method is to also reset the initial
    142       * step so it will be automatically computed by the integrator
    143       * if {@link #setInitialStepSize(double) setInitialStepSize}
    144       * is not called by the user.
    145       * </p>
    146       * @param minimalStep minimal step (must be positive even for backward
    147       * integration), the last step can be smaller than this
    148       * @param maximalStep maximal step (must be positive even for backward
    149       * integration)
    150       * @param absoluteTolerance allowed absolute error
    151       * @param relativeTolerance allowed relative error
    152       */
    153      public void setStepSizeControl(final double minimalStep, final double maximalStep,
    154                                     final double absoluteTolerance,
    155                                     final double relativeTolerance) {
    156
    157          minStep     = FastMath.abs(minimalStep);
    158          maxStep     = FastMath.abs(maximalStep);
    159          initialStep = -1;
    160
    161          scalAbsoluteTolerance = absoluteTolerance;
    162          scalRelativeTolerance = relativeTolerance;
    163          vecAbsoluteTolerance  = null;
    164          vecRelativeTolerance  = null;
    165
    166      }
    167
    168      /** Set the adaptive step size control parameters.
    169       * <p>
    170       * A side effect of this method is to also reset the initial
    171       * step so it will be automatically computed by the integrator
    172       * if {@link #setInitialStepSize(double) setInitialStepSize}
    173       * is not called by the user.
    174       * </p>
    175       * @param minimalStep minimal step (must be positive even for backward
    176       * integration), the last step can be smaller than this
    177       * @param maximalStep maximal step (must be positive even for backward
    178       * integration)
    179       * @param absoluteTolerance allowed absolute error
    180       * @param relativeTolerance allowed relative error
    181       */
    182      public void setStepSizeControl(final double minimalStep, final double maximalStep,
    183                                     final double[] absoluteTolerance,
    184                                     final double[] relativeTolerance) {
    185
    186          minStep     = FastMath.abs(minimalStep);
    187          maxStep     = FastMath.abs(maximalStep);
    188          initialStep = -1;
    189
    190          scalAbsoluteTolerance = 0;
    191          scalRelativeTolerance = 0;
    192          vecAbsoluteTolerance  = absoluteTolerance.clone();
    193          vecRelativeTolerance  = relativeTolerance.clone();
    194
    195      }
    196
    197      /** Set the initial step size.
    198       * <p>This method allows the user to specify an initial positive
    199       * step size instead of letting the integrator guess it by
    200       * itself. If this method is not called before integration is
    201       * started, the initial step size will be estimated by the
    202       * integrator.</p>
    203       * @param initialStepSize initial step size to use (must be positive even
    204       * for backward integration ; providing a negative value or a value
    205       * outside of the min/max step interval will lead the integrator to
    206       * ignore the value and compute the initial step size by itself)
    207       */
    208      public void setInitialStepSize(final double initialStepSize) {
    209        if ((initialStepSize < minStep) || (initialStepSize > maxStep)) {
    210          initialStep = -1.0;
    211        } else {
    212          initialStep = initialStepSize;
    213        }
    214      }
    215
    216      /** {@inheritDoc} */
    217      @Override
    218      protected void sanityChecks(final ExpandableStatefulODE equations, final double t)
    219          throws DimensionMismatchException, NumberIsTooSmallException {
    220
    221          super.sanityChecks(equations, t);
    222
    223          mainSetDimension = equations.getPrimaryMapper().getDimension();
    224
    225          if ((vecAbsoluteTolerance != null) && (vecAbsoluteTolerance.length != mainSetDimension)) {
    226              throw new DimensionMismatchException(mainSetDimension, vecAbsoluteTolerance.length);
    227          }
    228
    229          if ((vecRelativeTolerance != null) && (vecRelativeTolerance.length != mainSetDimension)) {
    230              throw new DimensionMismatchException(mainSetDimension, vecRelativeTolerance.length);
    231          }
    232
    233      }
    234
    235      /** Initialize the integration step.
    236       * @param forward forward integration indicator
    237       * @param order order of the method
    238       * @param scale scaling vector for the state vector (can be shorter than state vector)
    239       * @param t0 start time
    240       * @param y0 state vector at t0
    241       * @param yDot0 first time derivative of y0
    242       * @param y1 work array for a state vector
    243       * @param yDot1 work array for the first time derivative of y1
    244       * @return first integration step
    245       * @exception MaxCountExceededException if the number of functions evaluations is exceeded
    246       * @exception DimensionMismatchException if arrays dimensions do not match equations settings
    247       */
    248      public double initializeStep(final boolean forward, final int order, final double[] scale,
    249                                   final double t0, final double[] y0, final double[] yDot0,
    250                                   final double[] y1, final double[] yDot1)
    251          throws MaxCountExceededException, DimensionMismatchException {
    252
    253        if (initialStep > 0) {
    254          // use the user provided value
    255          return forward ? initialStep : -initialStep;
    256        }
    257
    258        // very rough first guess : h = 0.01 * ||y/scale|| / ||y'/scale||
    259        // this guess will be used to perform an Euler step
    260        double ratio;
    261        double yOnScale2 = 0;
    262        double yDotOnScale2 = 0;
    263        for (int j = 0; j < scale.length; ++j) {
    264          ratio         = y0[j] / scale[j];
    265          yOnScale2    += ratio * ratio;
    266          ratio         = yDot0[j] / scale[j];
    267          yDotOnScale2 += ratio * ratio;
    268        }
    269
    270        double h = ((yOnScale2 < 1.0e-10) || (yDotOnScale2 < 1.0e-10)) ?
    271                   1.0e-6 : (0.01 * FastMath.sqrt(yOnScale2 / yDotOnScale2));
    272        if (! forward) {
    273          h = -h;
    274        }
    275
    276        // perform an Euler step using the preceding rough guess
    277        for (int j = 0; j < y0.length; ++j) {
    278          y1[j] = y0[j] + h * yDot0[j];
    279        }
    280        computeDerivatives(t0 + h, y1, yDot1);
    281
    282        // estimate the second derivative of the solution
    283        double yDDotOnScale = 0;
    284        for (int j = 0; j < scale.length; ++j) {
    285          ratio         = (yDot1[j] - yDot0[j]) / scale[j];
    286          yDDotOnScale += ratio * ratio;
    287        }
    288        yDDotOnScale = FastMath.sqrt(yDDotOnScale) / h;
    289
    290        // step size is computed such that
    291        // h^order * max (||y'/tol||, ||y''/tol||) = 0.01
    292        final double maxInv2 = FastMath.max(FastMath.sqrt(yDotOnScale2), yDDotOnScale);
    293        final double h1 = (maxInv2 < 1.0e-15) ?
    294                          FastMath.max(1.0e-6, 0.001 * FastMath.abs(h)) :
    295                          FastMath.pow(0.01 / maxInv2, 1.0 / order);
    296        h = FastMath.min(100.0 * FastMath.abs(h), h1);
    297        h = FastMath.max(h, 1.0e-12 * FastMath.abs(t0));  // avoids cancellation when computing t1 - t0
    298        if (h < getMinStep()) {
    299          h = getMinStep();
    300        }
    301        if (h > getMaxStep()) {
    302          h = getMaxStep();
    303        }
    304        if (! forward) {
    305          h = -h;
    306        }
    307
    308        return h;
    309
    310      }
    311
    312      /** Filter the integration step.
    313       * @param h signed step
    314       * @param forward forward integration indicator
    315       * @param acceptSmall if true, steps smaller than the minimal value
    316       * are silently increased up to this value, if false such small
    317       * steps generate an exception
    318       * @return a bounded integration step (h if no bound is reach, or a bounded value)
    319       * @exception NumberIsTooSmallException if the step is too small and acceptSmall is false
    320       */
    321      protected double filterStep(final double h, final boolean forward, final boolean acceptSmall)
    322        throws NumberIsTooSmallException {
    323
    324          double filteredH = h;
    325          if (FastMath.abs(h) < minStep) {
    326              if (acceptSmall) {
    327                  filteredH = forward ? minStep : -minStep;
    328              } else {
    329                  throw new NumberIsTooSmallException(LocalizedFormats.MINIMAL_STEPSIZE_REACHED_DURING_INTEGRATION,
    330                                                      FastMath.abs(h), minStep, true);
    331              }
    332          }
    333
    334          if (filteredH > maxStep) {
    335              filteredH = maxStep;
    336          } else if (filteredH < -maxStep) {
    337              filteredH = -maxStep;
    338          }
    339
    340          return filteredH;
    341
    342      }
    343
    344      /** {@inheritDoc} */
    345      @Override
    346      public abstract void integrate (ExpandableStatefulODE equations, double t)
    347          throws NumberIsTooSmallException, DimensionMismatchException,
    348                 MaxCountExceededException, NoBracketingException;
    349
    350      /** {@inheritDoc} */
    351      @Override
    352      public double getCurrentStepStart() {
    353        return stepStart;
    354      }
    355
    356      /** Reset internal state to dummy values. */
    357      protected void resetInternalState() {
    358        stepStart = Double.NaN;
    359        stepSize  = FastMath.sqrt(minStep * maxStep);
    360      }
    361
    362      /** Get the minimal step.
    363       * @return minimal step
    364       */
    365      public double getMinStep() {
    366        return minStep;
    367      }
    368
    369      /** Get the maximal step.
    370       * @return maximal step
    371       */
    372      public double getMaxStep() {
    373        return maxStep;
    374      }
    375
    376    }