← Back to branch summary

~ubuntu-branches/ubuntu/oneiric/commons-math/oneiric

« back to all changes in this revision

Viewing changes to src/java/org/apache/commons/math/ode/AdaptiveStepsizeIntegrator.java

  • Committer: Bazaar Package Importer
  • Author(s): Damien Raude-Morvan
  • Date: 2009-08-22 01:13:25 UTC
  • mfrom: (1.1.1 upstream)
  • Revision ID: james.westby@ubuntu.com-20090822011325-hi4peq1ua5weguwn
Tags: 2.0-1
* New upstream release.
* Set Maintainer field to Debian Java Team
* Add myself as Uploaders
* Switch to Quilt patch system:
  - Refresh all patchs
  - Remove B-D on dpatch, Add B-D on quilt
  - Include patchsys-quilt.mk in debian/rules
* Bump Standards-Version to 3.8.3:
  - Add a README.source to describe patch system
* Maven POMs:
  - Add a Build-Depends-Indep dependency on maven-repo-helper
  - Use mh_installpom and mh_installjar to install the POM and the jar to the
    Maven repository
* Use default-jdk/jre:
  - Depends on java5-runtime-headless
  - Build-Depends on default-jdk
  - Use /usr/lib/jvm/default-java as JAVA_HOME
* Move api documentation to /usr/share/doc/libcommons-math-java/api
* Build-Depends on junit4 instead of junit

Show diffs side-by-side

added added

removed removed

Lines of Context:
src/java/org/apache/commons/math/ode/AdaptiveStepsizeIntegrator.java
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
 
}