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