    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.math.ode.nonstiff;
    019
    020    import org.apache.commons.math.exception.util.LocalizedFormats;
    021    import org.apache.commons.math.ode.AbstractIntegrator;
    022    import org.apache.commons.math.ode.DerivativeException;
    023    import org.apache.commons.math.ode.ExtendedFirstOrderDifferentialEquations;
    024    import org.apache.commons.math.ode.FirstOrderDifferentialEquations;
    025    import org.apache.commons.math.ode.IntegratorException;
    026    import org.apache.commons.math.util.FastMath;
    027
    028    /**
    029     * This abstract class holds the common part of all adaptive
    030     * stepsize integrators for Ordinary Differential Equations.
    031     *
    032     * <p>These algorithms perform integration with stepsize control, which
    033     * means the user does not specify the integration step but rather a
    034     * tolerance on error. The error threshold is computed as
    035     * <pre>
    036     * threshold_i = absTol_i + relTol_i * max (abs (ym), abs (ym+1))
    037     * </pre>
    038     * where absTol_i is the absolute tolerance for component i of the
    039     * state vector and relTol_i is the relative tolerance for the same
    040     * component. The user can also use only two scalar values absTol and
    041     * relTol which will be used for all components.
    042     * </p>
    043     *
    044     * <p>If the Ordinary Differential Equations is an {@link ExtendedFirstOrderDifferentialEquations
    045     * extended ODE} rather than a {@link FirstOrderDifferentialEquations basic ODE},
    046     * then <em>only</em> the {@link ExtendedFirstOrderDifferentialEquations#getMainSetDimension()
    047     * main set} part of the state vector is used for stepsize control, not the complete
    048     * state vector.
    049     * </p>
    050     *
    051     * <p>If the estimated error for ym+1 is such that
    052     * <pre>
    053     * sqrt((sum (errEst_i / threshold_i)^2 ) / n) < 1
    054     * </pre>
    055     *
    056     * (where n is the main set dimension) then the step is accepted,
    057     * otherwise the step is rejected and a new attempt is made with a new
    058     * stepsize.</p>
    059     *
    060     * @version $Revision: 1073158 $ $Date: 2011-02-21 22:46:52 +0100 (lun. 21 f??vr. 2011) $
    061     * @since 1.2
    062     *
    063     */
    064
    065    public abstract class AdaptiveStepsizeIntegrator
    066      extends AbstractIntegrator {
    067
    068        /** Allowed absolute scalar error. */
    069        protected final double scalAbsoluteTolerance;
    070
    071        /** Allowed relative scalar error. */
    072        protected final double scalRelativeTolerance;
    073
    074        /** Allowed absolute vectorial error. */
    075        protected final double[] vecAbsoluteTolerance;
    076
    077        /** Allowed relative vectorial error. */
    078        protected final double[] vecRelativeTolerance;
    079
    080        /** Main set dimension. */
    081        protected int mainSetDimension;
    082
    083        /** User supplied initial step. */
    084        private double initialStep;
    085
    086        /** Minimal step. */
    087        private final double minStep;
    088
    089        /** Maximal step. */
    090        private final double maxStep;
    091
    092      /** Build an integrator with the given stepsize bounds.
    093       * The default step handler does nothing.
    094       * @param name name of the method
    095       * @param minStep minimal step (must be positive even for backward
    096       * integration), the last step can be smaller than this
    097       * @param maxStep maximal step (must be positive even for backward
    098       * integration)
    099       * @param scalAbsoluteTolerance allowed absolute error
    100       * @param scalRelativeTolerance allowed relative error
    101       */
    102      public AdaptiveStepsizeIntegrator(final String name,
    103                                        final double minStep, final double maxStep,
    104                                        final double scalAbsoluteTolerance,
    105                                        final double scalRelativeTolerance) {
    106
    107        super(name);
    108
    109        this.minStep     = FastMath.abs(minStep);
    110        this.maxStep     = FastMath.abs(maxStep);
    111        this.initialStep = -1.0;
    112
    113        this.scalAbsoluteTolerance = scalAbsoluteTolerance;
    114        this.scalRelativeTolerance = scalRelativeTolerance;
    115        this.vecAbsoluteTolerance  = null;
    116        this.vecRelativeTolerance  = null;
    117
    118        resetInternalState();
    119
    120      }
    121
    122      /** Build an integrator with the given stepsize bounds.
    123       * The default step handler does nothing.
    124       * @param name name of the method
    125       * @param minStep minimal step (must be positive even for backward
    126       * integration), the last step can be smaller than this
    127       * @param maxStep maximal step (must be positive even for backward
    128       * integration)
    129       * @param vecAbsoluteTolerance allowed absolute error
    130       * @param vecRelativeTolerance allowed relative error
    131       */
    132      public AdaptiveStepsizeIntegrator(final String name,
    133                                        final double minStep, final double maxStep,
    134                                        final double[] vecAbsoluteTolerance,
    135                                        final double[] vecRelativeTolerance) {
    136
    137        super(name);
    138
    139        this.minStep     = minStep;
    140        this.maxStep     = maxStep;
    141        this.initialStep = -1.0;
    142
    143        this.scalAbsoluteTolerance = 0;
    144        this.scalRelativeTolerance = 0;
    145        this.vecAbsoluteTolerance  = vecAbsoluteTolerance.clone();
    146        this.vecRelativeTolerance  = vecRelativeTolerance.clone();
    147
    148        resetInternalState();
    149
    150      }
    151
    152      /** Set the initial step size.
    153       * <p>This method allows the user to specify an initial positive
    154       * step size instead of letting the integrator guess it by
    155       * itself. If this method is not called before integration is
    156       * started, the initial step size will be estimated by the
    157       * integrator.</p>
    158       * @param initialStepSize initial step size to use (must be positive even
    159       * for backward integration ; providing a negative value or a value
    160       * outside of the min/max step interval will lead the integrator to
    161       * ignore the value and compute the initial step size by itself)
    162       */
    163      public void setInitialStepSize(final double initialStepSize) {
    164        if ((initialStepSize < minStep) || (initialStepSize > maxStep)) {
    165          initialStep = -1.0;
    166        } else {
    167          initialStep = initialStepSize;
    168        }
    169      }
    170
    171      /** Perform some sanity checks on the integration parameters.
    172       * @param equations differential equations set
    173       * @param t0 start time
    174       * @param y0 state vector at t0
    175       * @param t target time for the integration
    176       * @param y placeholder where to put the state vector
    177       * @exception IntegratorException if some inconsistency is detected
    178       */
    179      @Override
    180      protected void sanityChecks(final FirstOrderDifferentialEquations equations,
    181                                  final double t0, final double[] y0,
    182                                  final double t, final double[] y)
    183          throws IntegratorException {
    184
    185          super.sanityChecks(equations, t0, y0, t, y);
    186
    187          if (equations instanceof ExtendedFirstOrderDifferentialEquations) {
    188              mainSetDimension = ((ExtendedFirstOrderDifferentialEquations) equations).getMainSetDimension();
    189          } else {
    190              mainSetDimension = equations.getDimension();
    191          }
    192
    193          if ((vecAbsoluteTolerance != null) && (vecAbsoluteTolerance.length != mainSetDimension)) {
    194              throw new IntegratorException(
    195                      LocalizedFormats.DIMENSIONS_MISMATCH_SIMPLE, mainSetDimension, vecAbsoluteTolerance.length);
    196          }
    197
    198          if ((vecRelativeTolerance != null) && (vecRelativeTolerance.length != mainSetDimension)) {
    199              throw new IntegratorException(
    200                      LocalizedFormats.DIMENSIONS_MISMATCH_SIMPLE, mainSetDimension, vecRelativeTolerance.length);
    201          }
    202
    203      }
    204
    205      /** Initialize the integration step.
    206       * @param equations differential equations set
    207       * @param forward forward integration indicator
    208       * @param order order of the method
    209       * @param scale scaling vector for the state vector (can be shorter than state vector)
    210       * @param t0 start time
    211       * @param y0 state vector at t0
    212       * @param yDot0 first time derivative of y0
    213       * @param y1 work array for a state vector
    214       * @param yDot1 work array for the first time derivative of y1
    215       * @return first integration step
    216       * @exception DerivativeException this exception is propagated to
    217       * the caller if the underlying user function triggers one
    218       */
    219      public double initializeStep(final FirstOrderDifferentialEquations equations,
    220                                   final boolean forward, final int order, final double[] scale,
    221                                   final double t0, final double[] y0, final double[] yDot0,
    222                                   final double[] y1, final double[] yDot1)
    223          throws DerivativeException {
    224
    225        if (initialStep > 0) {
    226          // use the user provided value
    227          return forward ? initialStep : -initialStep;
    228        }
    229
    230        // very rough first guess : h = 0.01 * ||y/scale|| / ||y'/scale||
    231        // this guess will be used to perform an Euler step
    232        double ratio;
    233        double yOnScale2 = 0;
    234        double yDotOnScale2 = 0;
    235        for (int j = 0; j < scale.length; ++j) {
    236          ratio         = y0[j] / scale[j];
    237          yOnScale2    += ratio * ratio;
    238          ratio         = yDot0[j] / scale[j];
    239          yDotOnScale2 += ratio * ratio;
    240        }
    241
    242        double h = ((yOnScale2 < 1.0e-10) || (yDotOnScale2 < 1.0e-10)) ?
    243                   1.0e-6 : (0.01 * FastMath.sqrt(yOnScale2 / yDotOnScale2));
    244        if (! forward) {
    245          h = -h;
    246        }
    247
    248        // perform an Euler step using the preceding rough guess
    249        for (int j = 0; j < y0.length; ++j) {
    250          y1[j] = y0[j] + h * yDot0[j];
    251        }
    252        computeDerivatives(t0 + h, y1, yDot1);
    253
    254        // estimate the second derivative of the solution
    255        double yDDotOnScale = 0;
    256        for (int j = 0; j < scale.length; ++j) {
    257          ratio         = (yDot1[j] - yDot0[j]) / scale[j];
    258          yDDotOnScale += ratio * ratio;
    259        }
    260        yDDotOnScale = FastMath.sqrt(yDDotOnScale) / h;
    261
    262        // step size is computed such that
    263        // h^order * max (||y'/tol||, ||y''/tol||) = 0.01
    264        final double maxInv2 = FastMath.max(FastMath.sqrt(yDotOnScale2), yDDotOnScale);
    265        final double h1 = (maxInv2 < 1.0e-15) ?
    266                          FastMath.max(1.0e-6, 0.001 * FastMath.abs(h)) :
    267                          FastMath.pow(0.01 / maxInv2, 1.0 / order);
    268        h = FastMath.min(100.0 * FastMath.abs(h), h1);
    269        h = FastMath.max(h, 1.0e-12 * FastMath.abs(t0));  // avoids cancellation when computing t1 - t0
    270        if (h < getMinStep()) {
    271          h = getMinStep();
    272        }
    273        if (h > getMaxStep()) {
    274          h = getMaxStep();
    275        }
    276        if (! forward) {
    277          h = -h;
    278        }
    279
    280        return h;
    281
    282      }
    283
    284      /** Filter the integration step.
    285       * @param h signed step
    286       * @param forward forward integration indicator
    287       * @param acceptSmall if true, steps smaller than the minimal value
    288       * are silently increased up to this value, if false such small
    289       * steps generate an exception
    290       * @return a bounded integration step (h if no bound is reach, or a bounded value)
    291       * @exception IntegratorException if the step is too small and acceptSmall is false
    292       */
    293      protected double filterStep(final double h, final boolean forward, final boolean acceptSmall)
    294        throws IntegratorException {
    295
    296          double filteredH = h;
    297          if (FastMath.abs(h) < minStep) {
    298              if (acceptSmall) {
    299                  filteredH = forward ? minStep : -minStep;
    300              } else {
    301                  throw new IntegratorException(
    302                          LocalizedFormats.MINIMAL_STEPSIZE_REACHED_DURING_INTEGRATION,
    303                          minStep, FastMath.abs(h));
    304              }
    305          }
    306
    307          if (filteredH > maxStep) {
    308              filteredH = maxStep;
    309          } else if (filteredH < -maxStep) {
    310              filteredH = -maxStep;
    311          }
    312
    313          return filteredH;
    314
    315      }
    316
    317      /** {@inheritDoc} */
    318      public abstract double integrate (FirstOrderDifferentialEquations equations,
    319                                        double t0, double[] y0,
    320                                        double t, double[] y)
    321        throws DerivativeException, IntegratorException;
    322
    323      /** {@inheritDoc} */
    324      @Override
    325      public double getCurrentStepStart() {
    326        return stepStart;
    327      }
    328
    329      /** Reset internal state to dummy values. */
    330      protected void resetInternalState() {
    331        stepStart = Double.NaN;
    332        stepSize  = FastMath.sqrt(minStep * maxStep);
    333      }
    334
    335      /** Get the minimal step.
    336       * @return minimal step
    337       */
    338      public double getMinStep() {
    339        return minStep;
    340      }
    341
    342      /** Get the maximal step.
    343       * @return maximal step
    344       */
    345      public double getMaxStep() {
    346        return maxStep;
    347      }
    348
    349    }