← 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/DormandPrince853Integrator.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/DormandPrince853Integrator.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
 
 
21
/**
 
22
 * This class implements the 8(5,3) Dormand-Prince integrator for Ordinary
 
23
 * Differential Equations.
 
24
 *
 
25
 * <p>This integrator is an embedded Runge-Kutta integrator
 
26
 * of order 8(5,3) used in local extrapolation mode (i.e. the solution
 
27
 * is computed using the high order formula) with stepsize control
 
28
 * (and automatic step initialization) and continuous output. This
 
29
 * method uses 12 functions evaluations per step for integration and 4
 
30
 * evaluations for interpolation. However, since the first
 
31
 * interpolation evaluation is the same as the first integration
 
32
 * evaluation of the next step, we have included it in the integrator
 
33
 * rather than in the interpolator and specified the method was an
 
34
 * <i>fsal</i>. Hence, despite we have 13 stages here, the cost is
 
35
 * really 12 evaluations per step even if no interpolation is done,
 
36
 * and the overcost of interpolation is only 3 evaluations.</p>
 
37
 *
 
38
 * <p>This method is based on an 8(6) method by Dormand and Prince
 
39
 * (i.e. order 8 for the integration and order 6 for error estimation)
 
40
 * modified by Hairer and Wanner to use a 5th order error estimator
 
41
 * with 3rd order correction. This modification was introduced because
 
42
 * the original method failed in some cases (wrong steps can be
 
43
 * accepted when step size is too large, for example in the
 
44
 * Brusselator problem) and also had <i>severe difficulties when
 
45
 * applied to problems with discontinuities</i>. This modification is
 
46
 * explained in the second edition of the first volume (Nonstiff
 
47
 * Problems) of the reference book by Hairer, Norsett and Wanner:
 
48
 * <i>Solving Ordinary Differential Equations</i> (Springer-Verlag,
 
49
 * ISBN 3-540-56670-8).</p>
 
50
 *
 
51
 * @version $Revision: 786881 $ $Date: 2009-06-20 14:53:08 -0400 (Sat, 20 Jun 2009) $
 
52
 * @since 1.2
 
53
 */
 
54
 
 
55
public class DormandPrince853Integrator extends EmbeddedRungeKuttaIntegrator {
 
56
 
 
57
  /** Integrator method name. */
 
58
  private static final String METHOD_NAME = "Dormand-Prince 8 (5, 3)";
 
59
 
 
60
  /** Time steps Butcher array. */
 
61
  private static final double[] staticC = {
 
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
    (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
    6.0/7.0, 1.0, 1.0
 
65
  };
 
66
 
 
67
  /** Internal weights Butcher array. */
 
68
  private static final double[][] staticA = {
 
69
 
 
70
    // k2
 
71
    {(12.0 - 2.0 * Math.sqrt(6.0)) / 135.0},
 
72
 
 
73
    // k3
 
74
    {(6.0 - Math.sqrt(6.0)) / 180.0, (6.0 - Math.sqrt(6.0)) / 60.0},
 
75
 
 
76
    // k4
 
77
    {(6.0 - Math.sqrt(6.0)) / 120.0, 0.0, (6.0 - Math.sqrt(6.0)) / 40.0},
 
78
 
 
79
    // k5
 
80
    {(462.0 + 107.0 * Math.sqrt(6.0)) / 3000.0, 0.0,
 
81
     (-402.0 - 197.0 * Math.sqrt(6.0)) / 1000.0, (168.0 + 73.0 * Math.sqrt(6.0)) / 375.0},
 
82
 
 
83
    // k6
 
84
    {1.0 / 27.0, 0.0, 0.0, (16.0 + Math.sqrt(6.0)) / 108.0, (16.0 - Math.sqrt(6.0)) / 108.0},
 
85
 
 
86
    // k7
 
87
    {19.0 / 512.0, 0.0, 0.0, (118.0 + 23.0 * Math.sqrt(6.0)) / 1024.0,
 
88
     (118.0 - 23.0 * Math.sqrt(6.0)) / 1024.0, -9.0 / 512.0},
 
89
 
 
90
    // k8
 
91
    {13772.0 / 371293.0, 0.0, 0.0, (51544.0 + 4784.0 * Math.sqrt(6.0)) / 371293.0,
 
92
     (51544.0 - 4784.0 * Math.sqrt(6.0)) / 371293.0, -5688.0 / 371293.0, 3072.0 / 371293.0},
 
93
 
 
94
    // k9
 
95
    {58656157643.0 / 93983540625.0, 0.0, 0.0,
 
96
     (-1324889724104.0 - 318801444819.0 * Math.sqrt(6.0)) / 626556937500.0,
 
97
     (-1324889724104.0 + 318801444819.0 * Math.sqrt(6.0)) / 626556937500.0,
 
98
     96044563816.0 / 3480871875.0, 5682451879168.0 / 281950621875.0,
 
99
     -165125654.0 / 3796875.0},
 
100
 
 
101
    // k10
 
102
    {8909899.0 / 18653125.0, 0.0, 0.0,
 
103
     (-4521408.0 - 1137963.0 * Math.sqrt(6.0)) / 2937500.0,
 
104
     (-4521408.0 + 1137963.0 * Math.sqrt(6.0)) / 2937500.0,
 
105
     96663078.0 / 4553125.0, 2107245056.0 / 137915625.0,
 
106
     -4913652016.0 / 147609375.0, -78894270.0 / 3880452869.0},
 
107
 
 
108
    // k11
 
109
    {-20401265806.0 / 21769653311.0, 0.0, 0.0,
 
110
     (354216.0 + 94326.0 * Math.sqrt(6.0)) / 112847.0,
 
111
     (354216.0 - 94326.0 * Math.sqrt(6.0)) / 112847.0,
 
112
     -43306765128.0 / 5313852383.0, -20866708358144.0 / 1126708119789.0,
 
113
     14886003438020.0 / 654632330667.0, 35290686222309375.0 / 14152473387134411.0,
 
114
     -1477884375.0 / 485066827.0},
 
115
 
 
116
    // k12
 
117
    {39815761.0 / 17514443.0, 0.0, 0.0,
 
118
     (-3457480.0 - 960905.0 * Math.sqrt(6.0)) / 551636.0,
 
119
     (-3457480.0 + 960905.0 * Math.sqrt(6.0)) / 551636.0,
 
120
     -844554132.0 / 47026969.0, 8444996352.0 / 302158619.0,
 
121
     -2509602342.0 / 877790785.0, -28388795297996250.0 / 3199510091356783.0,
 
122
     226716250.0 / 18341897.0, 1371316744.0 / 2131383595.0},
 
123
 
 
124
    // k13 should be for interpolation only, but since it is the same
 
125
    // stage as the first evaluation of the next step, we perform it
 
126
    // here at no cost by specifying this is an fsal method
 
127
    {104257.0/1920240.0, 0.0, 0.0, 0.0, 0.0, 3399327.0/763840.0,
 
128
     66578432.0/35198415.0, -1674902723.0/288716400.0,
 
129
     54980371265625.0/176692375811392.0, -734375.0/4826304.0,
 
130
     171414593.0/851261400.0, 137909.0/3084480.0}
 
131
 
 
132
  };
 
133
 
 
134
  /** Propagation weights Butcher array. */
 
135
  private static final double[] staticB = {
 
136
      104257.0/1920240.0,
 
137
      0.0,
 
138
      0.0,
 
139
      0.0,
 
140
      0.0,
 
141
      3399327.0/763840.0,
 
142
      66578432.0/35198415.0,
 
143
      -1674902723.0/288716400.0,
 
144
      54980371265625.0/176692375811392.0,
 
145
      -734375.0/4826304.0,
 
146
      171414593.0/851261400.0,
 
147
      137909.0/3084480.0,
 
148
      0.0
 
149
  };
 
150
 
 
151
  /** First error weights array, element 1. */
 
152
  private static final double e1_01 =         116092271.0 / 8848465920.0;
 
153
 
 
154
  // elements 2 to 5 are zero, so they are neither stored nor used
 
155
 
 
156
  /** First error weights array, element 6. */
 
157
  private static final double e1_06 =          -1871647.0 / 1527680.0;
 
158
 
 
159
  /** First error weights array, element 7. */
 
160
  private static final double e1_07 =         -69799717.0 / 140793660.0;
 
161
 
 
162
  /** First error weights array, element 8. */
 
163
  private static final double e1_08 =     1230164450203.0 / 739113984000.0;
 
164
 
 
165
  /** First error weights array, element 9. */
 
166
  private static final double e1_09 = -1980813971228885.0 / 5654156025964544.0;
 
167
 
 
168
  /** First error weights array, element 10. */
 
169
  private static final double e1_10 =         464500805.0 / 1389975552.0;
 
170
 
 
171
  /** First error weights array, element 11. */
 
172
  private static final double e1_11 =     1606764981773.0 / 19613062656000.0;
 
173
 
 
174
  /** First error weights array, element 12. */
 
175
  private static final double e1_12 =           -137909.0 / 6168960.0;
 
176
 
 
177
 
 
178
  /** Second error weights array, element 1. */
 
179
  private static final double e2_01 =           -364463.0 / 1920240.0;
 
180
 
 
181
  // elements 2 to 5 are zero, so they are neither stored nor used
 
182
 
 
183
  /** Second error weights array, element 6. */
 
184
  private static final double e2_06 =           3399327.0 / 763840.0;
 
185
 
 
186
  /** Second error weights array, element 7. */
 
187
  private static final double e2_07 =          66578432.0 / 35198415.0;
 
188
 
 
189
  /** Second error weights array, element 8. */
 
190
  private static final double e2_08 =       -1674902723.0 / 288716400.0;
 
191
 
 
192
  /** Second error weights array, element 9. */
 
193
  private static final double e2_09 =   -74684743568175.0 / 176692375811392.0;
 
194
 
 
195
  /** Second error weights array, element 10. */
 
196
  private static final double e2_10 =           -734375.0 / 4826304.0;
 
197
 
 
198
  /** Second error weights array, element 11. */
 
199
  private static final double e2_11 =         171414593.0 / 851261400.0;
 
200
 
 
201
  /** Second error weights array, element 12. */
 
202
  private static final double e2_12 =             69869.0 / 3084480.0;
 
203
 
 
204
  /** Simple constructor.
 
205
   * Build an eighth order Dormand-Prince integrator with the given step bounds
 
206
   * @param minStep minimal step (must be positive even for backward
 
207
   * integration), the last step can be smaller than this
 
208
   * @param maxStep maximal step (must be positive even for backward
 
209
   * integration)
 
210
   * @param scalAbsoluteTolerance allowed absolute error
 
211
   * @param scalRelativeTolerance allowed relative error
 
212
   */
 
213
  public DormandPrince853Integrator(final double minStep, final double maxStep,
 
214
                                    final double scalAbsoluteTolerance,
 
215
                                    final double scalRelativeTolerance) {
 
216
    super(METHOD_NAME, true, staticC, staticA, staticB,
 
217
          new DormandPrince853StepInterpolator(),
 
218
          minStep, maxStep, scalAbsoluteTolerance, scalRelativeTolerance);
 
219
  }
 
220
 
 
221
  /** Simple constructor.
 
222
   * Build an eighth order Dormand-Prince integrator with the given step bounds
 
223
   * @param minStep minimal step (must be positive even for backward
 
224
   * integration), the last step can be smaller than this
 
225
   * @param maxStep maximal step (must be positive even for backward
 
226
   * integration)
 
227
   * @param vecAbsoluteTolerance allowed absolute error
 
228
   * @param vecRelativeTolerance allowed relative error
 
229
   */
 
230
  public DormandPrince853Integrator(final double minStep, final double maxStep,
 
231
                                    final double[] vecAbsoluteTolerance,
 
232
                                    final double[] vecRelativeTolerance) {
 
233
    super(METHOD_NAME, true, staticC, staticA, staticB,
 
234
          new DormandPrince853StepInterpolator(),
 
235
          minStep, maxStep, vecAbsoluteTolerance, vecRelativeTolerance);
 
236
  }
 
237
 
 
238
  /** {@inheritDoc} */
 
239
  @Override
 
240
  public int getOrder() {
 
241
    return 8;
 
242
  }
 
243
 
 
244
  /** {@inheritDoc} */
 
245
  @Override
 
246
  protected double estimateError(final double[][] yDotK,
 
247
                                 final double[] y0, final double[] y1,
 
248
                                 final double h) {
 
249
    double error1 = 0;
 
250
    double error2 = 0;
 
251
 
 
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];
 
261
 
 
262
      final double yScale = Math.max(Math.abs(y0[j]), Math.abs(y1[j]));
 
263
      final double tol = (vecAbsoluteTolerance == null) ?
 
264
                         (scalAbsoluteTolerance + scalRelativeTolerance * yScale) :
 
265
                         (vecAbsoluteTolerance[j] + vecRelativeTolerance[j] * yScale);
 
266
      final double ratio1  = errSum1 / tol;
 
267
      error1        += ratio1 * ratio1;
 
268
      final double ratio2  = errSum2 / tol;
 
269
      error2        += ratio2 * ratio2;
 
270
    }
 
271
 
 
272
    double den = error1 + 0.01 * error2;
 
273
    if (den <= 0.0) {
 
274
      den = 1.0;
 
275
    }
 
276
 
 
277
    return Math.abs(h) * error1 / Math.sqrt(y0.length * den);
 
278
 
 
279
  }
 
280
 
 
281
}