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