← Back to branch summary

~ubuntu-branches/ubuntu/maverick/commons-math/maverick

« back to all changes in this revision

Viewing changes to src/main/java/org/apache/commons/math/ode/nonstiff/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/main/java/org/apache/commons/math/ode/nonstiff/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.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
}