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.nonstiff;
19  
20  import org.apache.commons.math.ode.AbstractIntegrator;
21  import org.apache.commons.math.ode.DerivativeException;
22  import org.apache.commons.math.ode.FirstOrderDifferentialEquations;
23  import org.apache.commons.math.ode.IntegratorException;
24  
25  /**
26   * This abstract class holds the common part of all adaptive
27   * stepsize integrators for Ordinary Differential Equations.
28   *
29   * <p>These algorithms perform integration with stepsize control, which
30   * means the user does not specify the integration step but rather a
31   * tolerance on error. The error threshold is computed as
32   * <pre>
33   * threshold_i = absTol_i + relTol_i * max (abs (ym), abs (ym+1))
34   * </pre>
35   * where absTol_i is the absolute tolerance for component i of the
36   * state vector and relTol_i is the relative tolerance for the same
37   * component. The user can also use only two scalar values absTol and
38   * relTol which will be used for all components.</p>
39   *
40   * <p>If the estimated error for ym+1 is such that
41   * <pre>
42   * sqrt((sum (errEst_i / threshold_i)^2 ) / n) < 1
43   * </pre>
44   *
45   * (where n is the state vector dimension) then the step is accepted,
46   * otherwise the step is rejected and a new attempt is made with a new
47   * stepsize.</p>
48   *
49   * @version $Revision: 795591 $ $Date: 2009-07-19 14:36:46 -0400 (Sun, 19 Jul 2009) $
50   * @since 1.2
51   *
52   */
53  
54  public abstract class AdaptiveStepsizeIntegrator
55    extends AbstractIntegrator {
56  
57    
58    /** Build an integrator with the given stepsize bounds.
59     * The default step handler does nothing.
60     * @param name name of the method
61     * @param minStep minimal step (must be positive even for backward
62     * integration), the last step can be smaller than this
63     * @param maxStep maximal step (must be positive even for backward
64     * integration)
65     * @param scalAbsoluteTolerance allowed absolute error
66     * @param scalRelativeTolerance allowed relative error
67     */
68    public AdaptiveStepsizeIntegrator(final String name,
69                                      final double minStep, final double maxStep,
70                                      final double scalAbsoluteTolerance,
71                                      final double scalRelativeTolerance) {
72  
73      super(name);
74  
75      this.minStep     = Math.abs(minStep);
76      this.maxStep     = Math.abs(maxStep);
77      this.initialStep = -1.0;
78  
79      this.scalAbsoluteTolerance = scalAbsoluteTolerance;
80      this.scalRelativeTolerance = scalRelativeTolerance;
81      this.vecAbsoluteTolerance  = null;
82      this.vecRelativeTolerance  = null;
83  
84      resetInternalState();
85  
86    }
87  
88    /** Build an integrator with the given stepsize bounds.
89     * The default step handler does nothing.
90     * @param name name of the method
91     * @param minStep minimal step (must be positive even for backward
92     * integration), the last step can be smaller than this
93     * @param maxStep maximal step (must be positive even for backward
94     * integration)
95     * @param vecAbsoluteTolerance allowed absolute error
96     * @param vecRelativeTolerance allowed relative error
97     */
98    public AdaptiveStepsizeIntegrator(final String name,
99                                      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 }