View Javadoc

1   /*
2    * Licensed to the Apache Software Foundation (ASF) under one or more
3    * contributor license agreements.  See the NOTICE file distributed with
4    * this work for additional information regarding copyright ownership.
5    * The ASF licenses this file to You under the Apache License, Version 2.0
6    * (the "License"); you may not use this file except in compliance with
7    * the License.  You may obtain a copy of the License at
8    *
9    *      http://www.apache.org/licenses/LICENSE-2.0
10   *
11   * Unless required by applicable law or agreed to in writing, software
12   * distributed under the License is distributed on an "AS IS" BASIS,
13   * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
14   * See the License for the specific language governing permissions and
15   * limitations under the License.
16   */
17  
18  package org.apache.commons.math.ode;
19  
20  /**
21   * This abstract class holds the common part of all adaptive
22   * stepsize integrators for Ordinary Differential Equations.
23   *
24   * <p>These algorithms perform integration with stepsize control, which
25   * means the user does not specify the integration step but rather a
26   * tolerance on error. The error threshold is computed as
27   * <pre>
28   * threshold_i = absTol_i + relTol_i * max (abs (ym), abs (ym+1))
29   * </pre>
30   * where absTol_i is the absolute tolerance for component i of the
31   * state vector and relTol_i is the relative tolerance for the same
32   * component. The user can also use only two scalar values absTol and
33   * relTol which will be used for all components.</p>
34   *
35   * <p>If the estimated error for ym+1 is such that
36   * <pre>
37   * sqrt((sum (errEst_i / threshold_i)^2 ) / n) < 1
38   * </pre>
39   *
40   * (where n is the state vector dimension) then the step is accepted,
41   * otherwise the step is rejected and a new attempt is made with a new
42   * stepsize.</p>
43   *
44   * @version $Revision: 620312 $ $Date: 2008-02-10 12:28:59 -0700 (Sun, 10 Feb 2008) $
45   * @since 1.2
46   *
47   */
48  
49  public abstract class AdaptiveStepsizeIntegrator
50    implements FirstOrderIntegrator {
51  
52    /** Build an integrator with the given stepsize bounds.
53     * The default step handler does nothing.
54     * @param minStep minimal step (must be positive even for backward
55     * integration), the last step can be smaller than this
56     * @param maxStep maximal step (must be positive even for backward
57     * integration)
58     * @param scalAbsoluteTolerance allowed absolute error
59     * @param scalRelativeTolerance allowed relative error
60     */
61    public AdaptiveStepsizeIntegrator(double minStep, double maxStep,
62                                      double scalAbsoluteTolerance,
63                                      double scalRelativeTolerance) {
64  
65      this.minStep     = minStep;
66      this.maxStep     = maxStep;
67      this.initialStep = -1.0;
68  
69      this.scalAbsoluteTolerance = scalAbsoluteTolerance;
70      this.scalRelativeTolerance = scalRelativeTolerance;
71      this.vecAbsoluteTolerance  = null;
72      this.vecRelativeTolerance  = null;
73  
74      // set the default step handler
75      handler = DummyStepHandler.getInstance();
76  
77      switchesHandler = new SwitchingFunctionsHandler();
78  
79      resetInternalState();
80  
81    }
82  
83    /** Build an integrator with the given stepsize bounds.
84     * The default step handler does nothing.
85     * @param minStep minimal step (must be positive even for backward
86     * integration), the last step can be smaller than this
87     * @param maxStep maximal step (must be positive even for backward
88     * integration)
89     * @param vecAbsoluteTolerance allowed absolute error
90     * @param vecRelativeTolerance allowed relative error
91     */
92    public AdaptiveStepsizeIntegrator(double minStep, double maxStep,
93                                      double[] vecAbsoluteTolerance,
94                                      double[] vecRelativeTolerance) {
95  
96      this.minStep     = minStep;
97      this.maxStep     = maxStep;
98      this.initialStep = -1.0;
99  
100     this.scalAbsoluteTolerance = 0;
101     this.scalRelativeTolerance = 0;
102     this.vecAbsoluteTolerance  = vecAbsoluteTolerance;
103     this.vecRelativeTolerance  = vecRelativeTolerance;
104 
105     // set the default step handler
106     handler = DummyStepHandler.getInstance();
107 
108     switchesHandler = new SwitchingFunctionsHandler();
109 
110     resetInternalState();
111 
112   }
113 
114   /** Set the initial step size.
115    * <p>This method allows the user to specify an initial positive
116    * step size instead of letting the integrator guess it by
117    * itself. If this method is not called before integration is
118    * started, the initial step size will be estimated by the
119    * integrator.</p>
120    * @param initialStepSize initial step size to use (must be positive even
121    * for backward integration ; providing a negative value or a value
122    * outside of the min/max step interval will lead the integrator to
123    * ignore the value and compute the initial step size by itself)
124    */
125   public void setInitialStepSize(double initialStepSize) {
126     if ((initialStepSize < minStep) || (initialStepSize > maxStep)) {
127       initialStep = -1.0;
128     } else {
129       initialStep = initialStepSize;
130     }
131   }
132 
133   /** Set the step handler for this integrator.
134    * The handler will be called by the integrator for each accepted
135    * step.
136    * @param handler handler for the accepted steps
137    */
138   public void setStepHandler (StepHandler handler) {
139     this.handler = handler;
140   }
141 
142   /** Get the step handler for this integrator.
143    * @return the step handler for this integrator
144    */
145   public StepHandler getStepHandler() {
146     return handler;
147   }
148 
149   /** Add a switching function to the integrator.
150    * @param function switching function
151    * @param maxCheckInterval maximal time interval between switching
152    * function checks (this interval prevents missing sign changes in
153    * case the integration steps becomes very large)
154    * @param convergence convergence threshold in the event time search
155    * @param maxIterationCount upper limit of the iteration count in
156    * the event time search
157    */
158   public void addSwitchingFunction(SwitchingFunction function,
159                                    double maxCheckInterval,
160                                    double convergence,
161                                    int maxIterationCount) {
162     switchesHandler.add(function, maxCheckInterval, convergence, maxIterationCount);
163   }
164 
165   /** Perform some sanity checks on the integration parameters.
166    * @param equations differential equations set
167    * @param t0 start time
168    * @param y0 state vector at t0
169    * @param t target time for the integration
170    * @param y placeholder where to put the state vector
171    * @exception IntegratorException if some inconsistency is detected
172    */
173   protected void sanityChecks(FirstOrderDifferentialEquations equations,
174                               double t0, double[] y0, double t, double[] y)
175     throws IntegratorException {
176       if (equations.getDimension() != y0.length) {
177           throw new IntegratorException("dimensions mismatch: ODE problem has dimension {0}," +
178                                         " initial state vector has dimension {1}",
179                                         new Object[] {
180                                           new Integer(equations.getDimension()),
181                                           new Integer(y0.length)
182                                         });
183       }
184       if (equations.getDimension() != y.length) {
185           throw new IntegratorException("dimensions mismatch: ODE problem has dimension {0}," +
186                                         " final state vector has dimension {1}",
187                                         new Object[] {
188                                           new Integer(equations.getDimension()),
189                                           new Integer(y.length)
190                                         });
191       }
192       if ((vecAbsoluteTolerance != null) && (vecAbsoluteTolerance.length != y0.length)) {
193           throw new IntegratorException("dimensions mismatch: state vector has dimension {0}," +
194                                         " absolute tolerance vector has dimension {1}",
195                                         new Object[] {
196                                           new Integer(y0.length),
197                                           new Integer(vecAbsoluteTolerance.length)
198                                         });
199       }
200       if ((vecRelativeTolerance != null) && (vecRelativeTolerance.length != y0.length)) {
201           throw new IntegratorException("dimensions mismatch: state vector has dimension {0}," +
202                                         " relative tolerance vector has dimension {1}",
203                                         new Object[] {
204                                           new Integer(y0.length),
205                                           new Integer(vecRelativeTolerance.length)
206                                         });
207       }
208       if (Math.abs(t - t0) <= 1.0e-12 * Math.max(Math.abs(t0), Math.abs(t))) {
209         throw new IntegratorException("too small integration interval: length = {0}",
210                                       new Object[] { new Double(Math.abs(t - t0)) });
211       }
212       
213   }
214 
215   /** Initialize the integration step.
216    * @param equations differential equations set
217    * @param forward forward integration indicator
218    * @param order order of the method
219    * @param scale scaling vector for the state vector
220    * @param t0 start time
221    * @param y0 state vector at t0
222    * @param yDot0 first time derivative of y0
223    * @param y1 work array for a state vector
224    * @param yDot1 work array for the first time derivative of y1
225    * @return first integration step
226    * @exception DerivativeException this exception is propagated to
227    * the caller if the underlying user function triggers one
228    */
229   public double initializeStep(FirstOrderDifferentialEquations equations,
230                                boolean forward, int order, double[] scale,
231                                double t0, double[] y0, double[] yDot0,
232                                double[] y1, double[] yDot1)
233     throws DerivativeException {
234 
235     if (initialStep > 0) {
236       // use the user provided value
237       return forward ? initialStep : -initialStep;
238     }
239 
240     // very rough first guess : h = 0.01 * ||y/scale|| / ||y'/scale||
241     // this guess will be used to perform an Euler step
242     double ratio;
243     double yOnScale2 = 0;
244     double yDotOnScale2 = 0;
245     for (int j = 0; j < y0.length; ++j) {
246       ratio         = y0[j] / scale[j];
247       yOnScale2    += ratio * ratio;
248       ratio         = yDot0[j] / scale[j];
249       yDotOnScale2 += ratio * ratio;
250     }
251 
252     double h = ((yOnScale2 < 1.0e-10) || (yDotOnScale2 < 1.0e-10)) ?
253                1.0e-6 : (0.01 * Math.sqrt(yOnScale2 / yDotOnScale2));
254     if (! forward) {
255       h = -h;
256     }
257 
258     // perform an Euler step using the preceding rough guess
259     for (int j = 0; j < y0.length; ++j) {
260       y1[j] = y0[j] + h * yDot0[j];
261     }
262     equations.computeDerivatives(t0 + h, y1, yDot1);
263 
264     // estimate the second derivative of the solution
265     double yDDotOnScale = 0;
266     for (int j = 0; j < y0.length; ++j) {
267       ratio         = (yDot1[j] - yDot0[j]) / scale[j];
268       yDDotOnScale += ratio * ratio;
269     }
270     yDDotOnScale = Math.sqrt(yDDotOnScale) / h;
271 
272     // step size is computed such that
273     // h^order * max (||y'/tol||, ||y''/tol||) = 0.01
274     double maxInv2 = Math.max(Math.sqrt(yDotOnScale2), yDDotOnScale);
275     double h1 = (maxInv2 < 1.0e-15) ?
276                 Math.max(1.0e-6, 0.001 * Math.abs(h)) :
277                 Math.pow(0.01 / maxInv2, 1.0 / order);
278     h = Math.min(100.0 * Math.abs(h), h1);
279     h = Math.max(h, 1.0e-12 * Math.abs(t0));  // avoids cancellation when computing t1 - t0
280     if (h < getMinStep()) {
281       h = getMinStep();
282     }
283     if (h > getMaxStep()) {
284       h = getMaxStep();
285     }
286     if (! forward) {
287       h = -h;
288     }
289 
290     return h;
291 
292   }
293 
294   /** Filter the integration step.
295    * @param h signed step
296    * @param acceptSmall if true, steps smaller than the minimal value
297    * are silently increased up to this value, if false such small
298    * steps generate an exception
299    * @return a bounded integration step (h if no bound is reach, or a bounded value)
300    * @exception IntegratorException if the step is too small and acceptSmall is false
301    */
302   protected double filterStep(double h, boolean acceptSmall)
303     throws IntegratorException {
304 
305     if (Math.abs(h) < minStep) {
306       if (acceptSmall) {
307         h = (h < 0) ? -minStep : minStep;
308       } else {
309         throw new IntegratorException("minimal step size ({0}) reached," +
310                                       " integration needs {1}",
311                                       new Object[] {
312                                         new Double(minStep),
313                                         new Double(Math.abs(h))
314                                       });
315       }
316     }
317 
318     if (h > maxStep) {
319       h = maxStep;
320     } else if (h < -maxStep) {
321       h = -maxStep;
322     }
323 
324     return h;
325 
326   }
327 
328   /** Integrate the differential equations up to the given time.
329    * <p>This method solves an Initial Value Problem (IVP).</p>
330    * <p>Since this method stores some internal state variables made
331    * available in its public interface during integration ({@link
332    * #getCurrentSignedStepsize()}), it is <em>not</em> thread-safe.</p>
333    * @param equations differential equations to integrate
334    * @param t0 initial time
335    * @param y0 initial value of the state vector at t0
336    * @param t target time for the integration
337    * (can be set to a value smaller than <code>t0</code> for backward integration)
338    * @param y placeholder where to put the state vector at each successful
339    *  step (and hence at the end of integration), can be the same object as y0
340    * @throws IntegratorException if the integrator cannot perform integration
341    * @throws DerivativeException this exception is propagated to the caller if
342    * the underlying user function triggers one
343    */
344   public abstract void integrate (FirstOrderDifferentialEquations equations,
345                                   double t0, double[] y0,
346                                   double t, double[] y)
347     throws DerivativeException, IntegratorException;
348 
349   /** Get the current value of the step start time t<sub>i</sub>.
350    * <p>This method can be called during integration (typically by
351    * the object implementing the {@link FirstOrderDifferentialEquations
352    * differential equations} problem) if the value of the current step that
353    * is attempted is needed.</p>
354    * <p>The result is undefined if the method is called outside of
355    * calls to {@link #integrate}</p>
356    * @return current value of the step start time t<sub>i</sub>
357    */
358   public double getCurrentStepStart() {
359     return stepStart;
360   }
361 
362   /** Get the current signed value of the integration stepsize.
363    * <p>This method can be called during integration (typically by
364    * the object implementing the {@link FirstOrderDifferentialEquations
365    * differential equations} problem) if the signed value of the current stepsize
366    * that is tried is needed.</p>
367    * <p>The result is undefined if the method is called outside of
368    * calls to {@link #integrate}</p>
369    * @return current signed value of the stepsize
370    */
371   public double getCurrentSignedStepsize() {
372     return stepSize;
373   }
374 
375   /** Reset internal state to dummy values. */
376   protected void resetInternalState() {
377     stepStart = Double.NaN;
378     stepSize  = Math.sqrt(minStep * maxStep);
379   }
380 
381   /** Get the minimal step.
382    * @return minimal step
383    */
384   public double getMinStep() {
385     return minStep;
386   }
387 
388   /** Get the maximal step.
389    * @return maximal step
390    */
391   public double getMaxStep() {
392     return maxStep;
393   }
394 
395   /** Minimal step. */
396   private double minStep;
397 
398   /** Maximal step. */
399   private double maxStep;
400 
401   /** User supplied initial step. */
402   private double initialStep;
403 
404   /** Allowed absolute scalar error. */
405   protected double scalAbsoluteTolerance;
406 
407   /** Allowed relative scalar error. */
408   protected double scalRelativeTolerance;
409 
410   /** Allowed absolute vectorial error. */
411   protected double[] vecAbsoluteTolerance;
412 
413   /** Allowed relative vectorial error. */
414   protected double[] vecRelativeTolerance;
415 
416   /** Step handler. */
417   protected StepHandler handler;
418 
419   /** Switching functions handler. */
420   protected SwitchingFunctionsHandler switchesHandler;
421 
422   /** Current step start time. */
423   protected double stepStart;
424 
425   /** Current stepsize. */
426   protected double stepSize;
427 
428 }