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    }