View Javadoc

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 java.io.IOException;
21  import java.io.ObjectInput;
22  import java.io.ObjectOutput;
23  
24  import org.apache.commons.math.MathRuntimeException;
25  import org.apache.commons.math.ode.AbstractIntegrator;
26  import org.apache.commons.math.ode.DerivativeException;
27  import org.apache.commons.math.ode.sampling.StepInterpolator;
28  
29  /**
30   * This class represents an interpolator over the last step during an
31   * ODE integration for the 8(5,3) Dormand-Prince integrator.
32   *
33   * @see DormandPrince853Integrator
34   *
35   * @version $Revision: 783103 $ $Date: 2009-06-09 15:33:19 -0400 (Tue, 09 Jun 2009) $
36   * @since 1.2
37   */
38  
39  class DormandPrince853StepInterpolator
40    extends RungeKuttaStepInterpolator {
41  
42    /** Simple constructor.
43     * This constructor builds an instance that is not usable yet, the
44     * {@link #reinitialize} method should be called before using the
45     * instance in order to initialize the internal arrays. This
46     * constructor is used only in order to delay the initialization in
47     * some cases. The {@link EmbeddedRungeKuttaIntegrator} uses the
48     * prototyping design pattern to create the step interpolators by
49     * cloning an uninitialized model and latter initializing the copy.
50     */
51    public DormandPrince853StepInterpolator() {
52      super();
53      yDotKLast = null;
54      v         = null;
55      vectorsInitialized = false;
56    }
57  
58    /** Copy constructor.
59     * @param interpolator interpolator to copy from. The copy is a deep
60     * copy: its arrays are separated from the original arrays of the
61     * instance
62     */
63    public DormandPrince853StepInterpolator(final DormandPrince853StepInterpolator interpolator) {
64  
65      super(interpolator);
66  
67      if (interpolator.currentState == null) {
68  
69        yDotKLast = null;
70        v         = null;
71        vectorsInitialized = false;
72  
73      } else {
74  
75        final int dimension = interpolator.currentState.length;
76  
77        yDotKLast    = new double[3][];
78        for (int k = 0; k < yDotKLast.length; ++k) {
79          yDotKLast[k] = new double[dimension];
80          System.arraycopy(interpolator.yDotKLast[k], 0, yDotKLast[k], 0,
81                           dimension);
82        }
83  
84        v = new double[7][];
85        for (int k = 0; k < v.length; ++k) {
86          v[k] = new double[dimension];
87          System.arraycopy(interpolator.v[k], 0, v[k], 0, dimension);
88        }
89  
90        vectorsInitialized = interpolator.vectorsInitialized;
91  
92      }
93  
94    }
95  
96    /** {@inheritDoc} */
97    @Override
98    protected StepInterpolator doCopy() {
99      return new DormandPrince853StepInterpolator(this);
100   }
101 
102   /** {@inheritDoc} */
103   @Override
104   public void reinitialize(final AbstractIntegrator integrator,
105                            final double[] y, final double[][] yDotK, final boolean forward) {
106 
107     super.reinitialize(integrator, y, yDotK, forward);
108 
109     final int dimension = currentState.length;
110 
111     yDotKLast = new double[3][];
112     for (int k = 0; k < yDotKLast.length; ++k) {
113       yDotKLast[k] = new double[dimension];
114     }
115 
116     v = new double[7][];
117     for (int k = 0; k < v.length; ++k) {
118       v[k]  = new double[dimension];
119     }
120 
121     vectorsInitialized = false;
122 
123   }
124 
125   /** {@inheritDoc} */
126   @Override
127   public void storeTime(final double t) {
128     super.storeTime(t);
129     vectorsInitialized = false;
130   }
131 
132   /** {@inheritDoc} */
133   @Override
134   protected void computeInterpolatedStateAndDerivatives(final double theta,
135                                           final double oneMinusThetaH)
136     throws DerivativeException {
137 
138     if (! vectorsInitialized) {
139 
140       if (v == null) {
141         v = new double[7][];
142         for (int k = 0; k < 7; ++k) {
143           v[k] = new double[interpolatedState.length];
144         }
145       }
146 
147       // perform the last evaluations if they have not been done yet
148       finalizeStep();
149 
150       // compute the interpolation vectors for this time step
151       for (int i = 0; i < interpolatedState.length; ++i) {
152           final double yDot1  = yDotK[0][i];
153           final double yDot6  = yDotK[5][i];
154           final double yDot7  = yDotK[6][i];
155           final double yDot8  = yDotK[7][i];
156           final double yDot9  = yDotK[8][i];
157           final double yDot10 = yDotK[9][i];
158           final double yDot11 = yDotK[10][i];
159           final double yDot12 = yDotK[11][i];
160           final double yDot13 = yDotK[12][i];
161           final double yDot14 = yDotKLast[0][i];
162           final double yDot15 = yDotKLast[1][i];
163           final double yDot16 = yDotKLast[2][i];
164           v[0][i] = b_01 * yDot1  + b_06 * yDot6 + b_07 * yDot7 +
165                     b_08 * yDot8  + b_09 * yDot9 + b_10 * yDot10 +
166                     b_11 * yDot11 + b_12 * yDot12;
167           v[1][i] = yDot1 - v[0][i];
168           v[2][i] = v[0][i] - v[1][i] - yDotK[12][i];
169           for (int k = 0; k < d.length; ++k) {
170               v[k+3][i] = d[k][0] * yDot1  + d[k][1]  * yDot6  + d[k][2]  * yDot7  +
171                           d[k][3] * yDot8  + d[k][4]  * yDot9  + d[k][5]  * yDot10 +
172                           d[k][6] * yDot11 + d[k][7]  * yDot12 + d[k][8]  * yDot13 +
173                           d[k][9] * yDot14 + d[k][10] * yDot15 + d[k][11] * yDot16;
174           }
175       }
176 
177       vectorsInitialized = true;
178 
179     }
180 
181     final double eta      = 1 - theta;
182     final double twoTheta = 2 * theta;
183     final double theta2   = theta * theta;
184     final double dot1 = 1 - twoTheta;
185     final double dot2 = theta * (2 - 3 * theta);
186     final double dot3 = twoTheta * (1 + theta * (twoTheta -3));
187     final double dot4 = theta2 * (3 + theta * (5 * theta - 8));
188     final double dot5 = theta2 * (3 + theta * (-12 + theta * (15 - 6 * theta)));
189     final double dot6 = theta2 * theta * (4 + theta * (-15 + theta * (18 - 7 * theta)));
190 
191     for (int i = 0; i < interpolatedState.length; ++i) {
192       interpolatedState[i] = currentState[i] -
193                              oneMinusThetaH * (v[0][i] -
194                                                theta * (v[1][i] +
195                                                         theta * (v[2][i] +
196                                                                  eta * (v[3][i] +
197                                                                         theta * (v[4][i] +
198                                                                                  eta * (v[5][i] +
199                                                                                         theta * (v[6][i])))))));
200       interpolatedDerivatives[i] =  v[0][i] + dot1 * v[1][i] + dot2 * v[2][i] +
201                                     dot3 * v[3][i] + dot4 * v[4][i] +
202                                     dot5 * v[5][i] + dot6 * v[6][i];
203     }
204 
205   }
206  
207   /** {@inheritDoc} */
208   @Override
209   protected void doFinalize()
210     throws DerivativeException {
211 
212     if (currentState == null) {
213       // we are finalizing an uninitialized instance
214       return;
215     }
216 
217     double s;
218     final double[] yTmp = new double[currentState.length];
219 
220     // k14
221     for (int j = 0; j < currentState.length; ++j) {
222       s = k14_01 * yDotK[0][j]  + k14_06 * yDotK[5][j]  + k14_07 * yDotK[6][j] +
223           k14_08 * yDotK[7][j]  + k14_09 * yDotK[8][j]  + k14_10 * yDotK[9][j] +
224           k14_11 * yDotK[10][j] + k14_12 * yDotK[11][j] + k14_13 * yDotK[12][j];
225       yTmp[j] = currentState[j] + h * s;
226     }
227     integrator.computeDerivatives(previousTime + c14 * h, yTmp, yDotKLast[0]);
228 
229     // k15
230     for (int j = 0; j < currentState.length; ++j) {
231      s = k15_01 * yDotK[0][j]  + k15_06 * yDotK[5][j]  + k15_07 * yDotK[6][j] +
232          k15_08 * yDotK[7][j]  + k15_09 * yDotK[8][j]  + k15_10 * yDotK[9][j] +
233          k15_11 * yDotK[10][j] + k15_12 * yDotK[11][j] + k15_13 * yDotK[12][j] +
234          k15_14 * yDotKLast[0][j];
235      yTmp[j] = currentState[j] + h * s;
236     }
237     integrator.computeDerivatives(previousTime + c15 * h, yTmp, yDotKLast[1]);
238 
239     // k16
240     for (int j = 0; j < currentState.length; ++j) {
241       s = k16_01 * yDotK[0][j]  + k16_06 * yDotK[5][j]  + k16_07 * yDotK[6][j] +
242           k16_08 * yDotK[7][j]  + k16_09 * yDotK[8][j]  + k16_10 * yDotK[9][j] +
243           k16_11 * yDotK[10][j] + k16_12 * yDotK[11][j] + k16_13 * yDotK[12][j] +
244           k16_14 * yDotKLast[0][j] +  k16_15 * yDotKLast[1][j];
245       yTmp[j] = currentState[j] + h * s;
246     }
247     integrator.computeDerivatives(previousTime + c16 * h, yTmp, yDotKLast[2]);
248 
249   }
250 
251   /** {@inheritDoc} */
252   @Override
253   public void writeExternal(final ObjectOutput out)
254     throws IOException {
255 
256     try {
257       // save the local attributes
258       finalizeStep();
259     } catch (DerivativeException e) {
260       throw MathRuntimeException.createIOException(e);
261     }
262     final int dimension = (currentState == null) ? -1 : currentState.length;
263     out.writeInt(dimension);
264     for (int i = 0; i < dimension; ++i) {
265       out.writeDouble(yDotKLast[0][i]);
266       out.writeDouble(yDotKLast[1][i]);
267       out.writeDouble(yDotKLast[2][i]);
268     }
269 
270     // save the state of the base class
271     super.writeExternal(out);
272 
273   }
274 
275   /** {@inheritDoc} */
276   @Override
277   public void readExternal(final ObjectInput in)
278     throws IOException {
279 
280     // read the local attributes
281     yDotKLast = new double[3][];
282     final int dimension = in.readInt();
283     yDotKLast[0] = (dimension < 0) ? null : new double[dimension];
284     yDotKLast[1] = (dimension < 0) ? null : new double[dimension];
285     yDotKLast[2] = (dimension < 0) ? null : new double[dimension];
286 
287     for (int i = 0; i < dimension; ++i) {
288       yDotKLast[0][i] = in.readDouble();
289       yDotKLast[1][i] = in.readDouble();
290       yDotKLast[2][i] = in.readDouble();
291     }
292 
293     // read the base state
294     super.readExternal(in);
295 
296   }
297 
298   /** Last evaluations. */
299   private double[][] yDotKLast;
300 
301   /** Vectors for interpolation. */
302   private double[][] v;
303 
304   /** Initialization indicator for the interpolation vectors. */
305   private boolean vectorsInitialized;
306 
307   /** Propagation weights, element 1. */
308   private static final double b_01 =         104257.0 / 1920240.0;
309 
310   // elements 2 to 5 are zero, so they are neither stored nor used
311 
312   /** Propagation weights, element 6. */
313   private static final double b_06 =        3399327.0 / 763840.0;
314 
315   /** Propagation weights, element 7. */
316   private static final double b_07 =       66578432.0 / 35198415.0;
317 
318   /** Propagation weights, element 8. */
319   private static final double b_08 =    -1674902723.0 / 288716400.0;
320 
321   /** Propagation weights, element 9. */
322   private static final double b_09 = 54980371265625.0 / 176692375811392.0;
323 
324   /** Propagation weights, element 10. */
325   private static final double b_10 =        -734375.0 / 4826304.0;
326 
327   /** Propagation weights, element 11. */
328   private static final double b_11 =      171414593.0 / 851261400.0;
329 
330   /** Propagation weights, element 12. */
331   private static final double b_12 =         137909.0 / 3084480.0;
332 
333   /** Time step for stage 14 (interpolation only). */
334   private static final double c14    = 1.0 / 10.0;
335 
336   /** Internal weights for stage 14, element 1. */
337   private static final double k14_01 =       13481885573.0 / 240030000000.0      - b_01;
338 
339   // elements 2 to 5 are zero, so they are neither stored nor used
340 
341   /** Internal weights for stage 14, element 6. */
342   private static final double k14_06 =                 0.0                       - b_06;
343 
344   /** Internal weights for stage 14, element 7. */
345   private static final double k14_07 =      139418837528.0 / 549975234375.0      - b_07;
346 
347   /** Internal weights for stage 14, element 8. */
348   private static final double k14_08 =   -11108320068443.0 / 45111937500000.0    - b_08;
349 
350   /** Internal weights for stage 14, element 9. */
351   private static final double k14_09 = -1769651421925959.0 / 14249385146080000.0 - b_09;
352 
353   /** Internal weights for stage 14, element 10. */
354   private static final double k14_10 =          57799439.0 / 377055000.0         - b_10;
355 
356   /** Internal weights for stage 14, element 11. */
357   private static final double k14_11 =      793322643029.0 / 96734250000000.0    - b_11;
358 
359   /** Internal weights for stage 14, element 12. */
360   private static final double k14_12 =        1458939311.0 / 192780000000.0      - b_12;
361 
362   /** Internal weights for stage 14, element 13. */
363   private static final double k14_13 =             -4149.0 / 500000.0;
364 
365   /** Time step for stage 15 (interpolation only). */
366   private static final double c15    = 1.0 / 5.0;
367 
368 
369   /** Internal weights for stage 15, element 1. */
370   private static final double k15_01 =     1595561272731.0 / 50120273500000.0    - b_01;
371 
372   // elements 2 to 5 are zero, so they are neither stored nor used
373 
374   /** Internal weights for stage 15, element 6. */
375   private static final double k15_06 =      975183916491.0 / 34457688031250.0    - b_06;
376 
377   /** Internal weights for stage 15, element 7. */
378   private static final double k15_07 =    38492013932672.0 / 718912673015625.0   - b_07;
379 
380   /** Internal weights for stage 15, element 8. */
381   private static final double k15_08 = -1114881286517557.0 / 20298710767500000.0 - b_08;
382 
383   /** Internal weights for stage 15, element 9. */
384   private static final double k15_09 =                 0.0                       - b_09;
385 
386   /** Internal weights for stage 15, element 10. */
387   private static final double k15_10 =                 0.0                       - b_10;
388 
389   /** Internal weights for stage 15, element 11. */
390   private static final double k15_11 =    -2538710946863.0 / 23431227861250000.0 - b_11;
391 
392   /** Internal weights for stage 15, element 12. */
393   private static final double k15_12 =        8824659001.0 / 23066716781250.0    - b_12;
394 
395   /** Internal weights for stage 15, element 13. */
396   private static final double k15_13 =      -11518334563.0 / 33831184612500.0;
397 
398   /** Internal weights for stage 15, element 14. */
399   private static final double k15_14 =        1912306948.0 / 13532473845.0;
400 
401   /** Time step for stage 16 (interpolation only). */
402   private static final double c16    = 7.0 / 9.0;
403 
404 
405   /** Internal weights for stage 16, element 1. */
406   private static final double k16_01 =      -13613986967.0 / 31741908048.0       - b_01;
407 
408   // elements 2 to 5 are zero, so they are neither stored nor used
409 
410   /** Internal weights for stage 16, element 6. */
411   private static final double k16_06 =       -4755612631.0 / 1012344804.0        - b_06;
412 
413   /** Internal weights for stage 16, element 7. */
414   private static final double k16_07 =    42939257944576.0 / 5588559685701.0     - b_07;
415 
416   /** Internal weights for stage 16, element 8. */
417   private static final double k16_08 =    77881972900277.0 / 19140370552944.0    - b_08;
418 
419   /** Internal weights for stage 16, element 9. */
420   private static final double k16_09 =    22719829234375.0 / 63689648654052.0    - b_09;
421 
422   /** Internal weights for stage 16, element 10. */
423   private static final double k16_10 =                 0.0                       - b_10;
424 
425   /** Internal weights for stage 16, element 11. */
426   private static final double k16_11 =                 0.0                       - b_11;
427 
428   /** Internal weights for stage 16, element 12. */
429   private static final double k16_12 =                 0.0                       - b_12;
430 
431   /** Internal weights for stage 16, element 13. */
432   private static final double k16_13 =       -1199007803.0 / 857031517296.0;
433 
434   /** Internal weights for stage 16, element 14. */
435   private static final double k16_14 =      157882067000.0 / 53564469831.0;
436 
437   /** Internal weights for stage 16, element 15. */
438   private static final double k16_15 =     -290468882375.0 / 31741908048.0;
439 
440   /** Interpolation weights.
441    * (beware that only the non-null values are in the table)
442    */
443   private static final double[][] d = {
444 
445     {        -17751989329.0 / 2106076560.0,               4272954039.0 / 7539864640.0,
446             -118476319744.0 / 38604839385.0,            755123450731.0 / 316657731600.0,
447       3692384461234828125.0 / 1744130441634250432.0,     -4612609375.0 / 5293382976.0,
448             2091772278379.0 / 933644586600.0,             2136624137.0 / 3382989120.0,
449                   -126493.0 / 1421424.0,                    98350000.0 / 5419179.0,
450                 -18878125.0 / 2053168.0,                 -1944542619.0 / 438351368.0},
451 
452     {         32941697297.0 / 3159114840.0,             456696183123.0 / 1884966160.0,
453            19132610714624.0 / 115814518155.0,       -177904688592943.0 / 474986597400.0,
454      -4821139941836765625.0 / 218016305204281304.0,      30702015625.0 / 3970037232.0,
455           -85916079474274.0 / 2800933759800.0,           -5919468007.0 / 634310460.0,
456                   2479159.0 / 157936.0,                    -18750000.0 / 602131.0,
457                 -19203125.0 / 2053168.0,                 15700361463.0 / 438351368.0},
458 
459     {         12627015655.0 / 631822968.0,              -72955222965.0 / 188496616.0,
460           -13145744952320.0 / 69488710893.0,          30084216194513.0 / 56998391688.0,
461       -296858761006640625.0 / 25648977082856624.0,         569140625.0 / 82709109.0,
462              -18684190637.0 / 18672891732.0,                69644045.0 / 89549712.0,
463                 -11847025.0 / 4264272.0,                  -978650000.0 / 16257537.0,
464                 519371875.0 / 6159504.0,                  5256837225.0 / 438351368.0},
465 
466     {          -450944925.0 / 17550638.0,               -14532122925.0 / 94248308.0,
467             -595876966400.0 / 2573655959.0,             188748653015.0 / 527762886.0,
468       2545485458115234375.0 / 27252038150535163.0,       -1376953125.0 / 36759604.0,
469               53995596795.0 / 518691437.0,                 210311225.0 / 7047894.0,
470                  -1718875.0 / 39484.0,                      58000000.0 / 602131.0,
471                  -1546875.0 / 39484.0,                   -1262172375.0 / 8429834.0}
472 
473   };
474 
475   /** Serializable version identifier */
476   private static final long serialVersionUID = 7152276390558450974L;
477 
478 }