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