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    }