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