← Back to branch summary

~ubuntu-branches/ubuntu/quantal/commons-math/quantal

« back to all changes in this revision

Viewing changes to src/main/java/org/apache/commons/math/ode/nonstiff/DormandPrince853Integrator.java

  • Committer: Bazaar Package Importer
  • Author(s): Damien Raude-Morvan
  • Date: 2010-04-05 23:33:02 UTC
  • mfrom: (1.1.2 upstream)
  • Revision ID: james.westby@ubuntu.com-20100405233302-gpqlceked76nw28a
Tags: 2.1-1
* New upstream release.
* Bump Standards-Version to 3.8.4: no changes needed
* Bump debhelper to >= 7
* Switch to 3.0 (quilt) source format:
  - Remove B-D on quilt
  - Add d/source/format
  - Remove d/README.source

Show diffs side-by-side

added added

removed removed

Lines of Context:
src/main/java/org/apache/commons/math/ode/nonstiff/DormandPrince853Integrator.java
48
48
 * <i>Solving Ordinary Differential Equations</i> (Springer-Verlag,
49
49
 * ISBN 3-540-56670-8).</p>
50
50
 *
51
 
 * @version $Revision: 786881 $ $Date: 2009-06-20 14:53:08 -0400 (Sat, 20 Jun 2009) $
 
51
 * @version $Revision: 810196 $ $Date: 2009-09-01 15:47:46 -0400 (Tue, 01 Sep 2009) $
52
52
 * @since 1.2
53
53
 */
54
54
 
58
58
  private static final String METHOD_NAME = "Dormand-Prince 8 (5, 3)";
59
59
 
60
60
  /** Time steps Butcher array. */
61
 
  private static final double[] staticC = {
 
61
  private static final double[] STATIC_C = {
62
62
    (12.0 - 2.0 * Math.sqrt(6.0)) / 135.0, (6.0 - Math.sqrt(6.0)) / 45.0, (6.0 - Math.sqrt(6.0)) / 30.0,
63
63
    (6.0 + Math.sqrt(6.0)) / 30.0, 1.0/3.0, 1.0/4.0, 4.0/13.0, 127.0/195.0, 3.0/5.0,
64
64
    6.0/7.0, 1.0, 1.0
65
65
  };
66
66
 
67
67
  /** Internal weights Butcher array. */
68
 
  private static final double[][] staticA = {
 
68
  private static final double[][] STATIC_A = {
69
69
 
70
70
    // k2
71
71
    {(12.0 - 2.0 * Math.sqrt(6.0)) / 135.0},
132
132
  };
133
133
 
134
134
  /** Propagation weights Butcher array. */
135
 
  private static final double[] staticB = {
 
135
  private static final double[] STATIC_B = {
136
136
      104257.0/1920240.0,
137
137
      0.0,
138
138
      0.0,
149
149
  };
150
150
 
151
151
  /** First error weights array, element 1. */
152
 
  private static final double e1_01 =         116092271.0 / 8848465920.0;
 
152
  private static final double E1_01 =         116092271.0 / 8848465920.0;
153
153
 
154
154
  // elements 2 to 5 are zero, so they are neither stored nor used
155
155
 
156
156
  /** First error weights array, element 6. */
157
 
  private static final double e1_06 =          -1871647.0 / 1527680.0;
 
157
  private static final double E1_06 =          -1871647.0 / 1527680.0;
158
158
 
159
159
  /** First error weights array, element 7. */
160
 
  private static final double e1_07 =         -69799717.0 / 140793660.0;
 
160
  private static final double E1_07 =         -69799717.0 / 140793660.0;
161
161
 
162
162
  /** First error weights array, element 8. */
163
 
  private static final double e1_08 =     1230164450203.0 / 739113984000.0;
 
163
  private static final double E1_08 =     1230164450203.0 / 739113984000.0;
164
164
 
165
165
  /** First error weights array, element 9. */
166
 
  private static final double e1_09 = -1980813971228885.0 / 5654156025964544.0;
 
166
  private static final double E1_09 = -1980813971228885.0 / 5654156025964544.0;
167
167
 
168
168
  /** First error weights array, element 10. */
169
 
  private static final double e1_10 =         464500805.0 / 1389975552.0;
 
169
  private static final double E1_10 =         464500805.0 / 1389975552.0;
170
170
 
171
171
  /** First error weights array, element 11. */
172
 
  private static final double e1_11 =     1606764981773.0 / 19613062656000.0;
 
172
  private static final double E1_11 =     1606764981773.0 / 19613062656000.0;
173
173
 
174
174
  /** First error weights array, element 12. */
175
 
  private static final double e1_12 =           -137909.0 / 6168960.0;
 
175
  private static final double E1_12 =           -137909.0 / 6168960.0;
176
176
 
177
177
 
178
178
  /** Second error weights array, element 1. */
179
 
  private static final double e2_01 =           -364463.0 / 1920240.0;
 
179
  private static final double E2_01 =           -364463.0 / 1920240.0;
180
180
 
181
181
  // elements 2 to 5 are zero, so they are neither stored nor used
182
182
 
183
183
  /** Second error weights array, element 6. */
184
 
  private static final double e2_06 =           3399327.0 / 763840.0;
 
184
  private static final double E2_06 =           3399327.0 / 763840.0;
185
185
 
186
186
  /** Second error weights array, element 7. */
187
 
  private static final double e2_07 =          66578432.0 / 35198415.0;
 
187
  private static final double E2_07 =          66578432.0 / 35198415.0;
188
188
 
189
189
  /** Second error weights array, element 8. */
190
 
  private static final double e2_08 =       -1674902723.0 / 288716400.0;
 
190
  private static final double E2_08 =       -1674902723.0 / 288716400.0;
191
191
 
192
192
  /** Second error weights array, element 9. */
193
 
  private static final double e2_09 =   -74684743568175.0 / 176692375811392.0;
 
193
  private static final double E2_09 =   -74684743568175.0 / 176692375811392.0;
194
194
 
195
195
  /** Second error weights array, element 10. */
196
 
  private static final double e2_10 =           -734375.0 / 4826304.0;
 
196
  private static final double E2_10 =           -734375.0 / 4826304.0;
197
197
 
198
198
  /** Second error weights array, element 11. */
199
 
  private static final double e2_11 =         171414593.0 / 851261400.0;
 
199
  private static final double E2_11 =         171414593.0 / 851261400.0;
200
200
 
201
201
  /** Second error weights array, element 12. */
202
 
  private static final double e2_12 =             69869.0 / 3084480.0;
 
202
  private static final double E2_12 =             69869.0 / 3084480.0;
203
203
 
204
204
  /** Simple constructor.
205
205
   * Build an eighth order Dormand-Prince integrator with the given step bounds
213
213
  public DormandPrince853Integrator(final double minStep, final double maxStep,
214
214
                                    final double scalAbsoluteTolerance,
215
215
                                    final double scalRelativeTolerance) {
216
 
    super(METHOD_NAME, true, staticC, staticA, staticB,
 
216
    super(METHOD_NAME, true, STATIC_C, STATIC_A, STATIC_B,
217
217
          new DormandPrince853StepInterpolator(),
218
218
          minStep, maxStep, scalAbsoluteTolerance, scalRelativeTolerance);
219
219
  }
230
230
  public DormandPrince853Integrator(final double minStep, final double maxStep,
231
231
                                    final double[] vecAbsoluteTolerance,
232
232
                                    final double[] vecRelativeTolerance) {
233
 
    super(METHOD_NAME, true, staticC, staticA, staticB,
 
233
    super(METHOD_NAME, true, STATIC_C, STATIC_A, STATIC_B,
234
234
          new DormandPrince853StepInterpolator(),
235
235
          minStep, maxStep, vecAbsoluteTolerance, vecRelativeTolerance);
236
236
  }
250
250
    double error2 = 0;
251
251
 
252
252
    for (int j = 0; j < y0.length; ++j) {
253
 
      final double errSum1 = e1_01 * yDotK[0][j]  + e1_06 * yDotK[5][j] +
254
 
                             e1_07 * yDotK[6][j]  + e1_08 * yDotK[7][j] +
255
 
                             e1_09 * yDotK[8][j]  + e1_10 * yDotK[9][j] +
256
 
                             e1_11 * yDotK[10][j] + e1_12 * yDotK[11][j];
257
 
      final double errSum2 = e2_01 * yDotK[0][j]  + e2_06 * yDotK[5][j] +
258
 
                             e2_07 * yDotK[6][j]  + e2_08 * yDotK[7][j] +
259
 
                             e2_09 * yDotK[8][j]  + e2_10 * yDotK[9][j] +
260
 
                             e2_11 * yDotK[10][j] + e2_12 * yDotK[11][j];
 
253
      final double errSum1 = E1_01 * yDotK[0][j]  + E1_06 * yDotK[5][j] +
 
254
                             E1_07 * yDotK[6][j]  + E1_08 * yDotK[7][j] +
 
255
                             E1_09 * yDotK[8][j]  + E1_10 * yDotK[9][j] +
 
256
                             E1_11 * yDotK[10][j] + E1_12 * yDotK[11][j];
 
257
      final double errSum2 = E2_01 * yDotK[0][j]  + E2_06 * yDotK[5][j] +
 
258
                             E2_07 * yDotK[6][j]  + E2_08 * yDotK[7][j] +
 
259
                             E2_09 * yDotK[8][j]  + E2_10 * yDotK[9][j] +
 
260
                             E2_11 * yDotK[10][j] + E2_12 * yDotK[11][j];
261
261
 
262
262
      final double yScale = Math.max(Math.abs(y0[j]), Math.abs(y1[j]));
263
263
      final double tol = (vecAbsoluteTolerance == null) ?