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
9
* http://www.apache.org/licenses/LICENSE-2.0
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.
18
package org.apache.commons.math.ode;
20
import java.io.ObjectOutput;
21
import java.io.ObjectInput;
22
import java.io.IOException;
25
* This class represents an interpolator over the last step during an
26
* ODE integration for the 8(5,3) Dormand-Prince integrator.
28
* @see DormandPrince853Integrator
30
* @version $Revision: 620312 $ $Date: 2008-02-10 12:28:59 -0700 (Sun, 10 Feb 2008) $
34
class DormandPrince853StepInterpolator
35
extends RungeKuttaStepInterpolator {
37
/** Simple constructor.
38
* This constructor builds an instance that is not usable yet, the
39
* {@link #reinitialize} method should be called before using the
40
* instance in order to initialize the internal arrays. This
41
* constructor is used only in order to delay the initialization in
42
* some cases. The {@link EmbeddedRungeKuttaIntegrator} uses the
43
* prototyping design pattern to create the step interpolators by
44
* cloning an uninitialized model and latter initializing the copy.
46
public DormandPrince853StepInterpolator() {
50
vectorsInitialized = false;
54
* @param interpolator interpolator to copy from. The copy is a deep
55
* copy: its arrays are separated from the original arrays of the
58
public DormandPrince853StepInterpolator(DormandPrince853StepInterpolator interpolator) {
62
if (interpolator.currentState == null) {
66
vectorsInitialized = false;
70
int dimension = interpolator.currentState.length;
72
yDotKLast = new double[3][];
73
for (int k = 0; k < yDotKLast.length; ++k) {
74
yDotKLast[k] = new double[dimension];
75
System.arraycopy(interpolator.yDotKLast[k], 0, yDotKLast[k], 0,
80
for (int k = 0; k < v.length; ++k) {
81
v[k] = new double[dimension];
82
System.arraycopy(interpolator.v[k], 0, v[k], 0, dimension);
85
vectorsInitialized = interpolator.vectorsInitialized;
91
/** Really copy the finalized instance.
92
* @return a copy of the finalized instance
94
protected StepInterpolator doCopy() {
95
return new DormandPrince853StepInterpolator(this);
98
/** Reinitialize the instance
99
* Some embedded Runge-Kutta integrators need fewer functions
100
* evaluations than their counterpart step interpolators. So the
101
* interpolator should perform the last evaluations they need by
102
* themselves. The {@link EmbeddedRungeKuttaIntegrator
103
* EmbeddedRungeKuttaIntegrator} abstract class calls this method in
104
* order to let the step interpolator perform the evaluations it
105
* needs. These evaluations will be performed during the call to
106
* <code>doFinalize</code> if any, i.e. only if the step handler
107
* either calls the {@link AbstractStepInterpolator#finalizeStep
108
* finalizeStep} method or the {@link
109
* AbstractStepInterpolator#getInterpolatedState getInterpolatedState}
110
* method (for an interpolator which needs a finalization) or if it clones
111
* the step interpolator.
112
* @param equations set of differential equations being integrated
113
* @param y reference to the integrator array holding the state at
114
* the end of the step
115
* @param yDotK reference to the integrator array holding all the
116
* intermediate slopes
117
* @param forward integration direction indicator
119
public void reinitialize(FirstOrderDifferentialEquations equations,
120
double[] y, double[][] yDotK, boolean forward) {
122
super.reinitialize(equations, y, yDotK, forward);
124
int dimension = currentState.length;
126
yDotKLast = new double[3][];
127
for (int k = 0; k < yDotKLast.length; ++k) {
128
yDotKLast[k] = new double[dimension];
132
for (int k = 0; k < v.length; ++k) {
133
v[k] = new double[dimension];
136
vectorsInitialized = false;
140
/** Store the current step time.
141
* @param t current time
143
public void storeTime(double t) {
145
vectorsInitialized = false;
148
/** Compute the state at the interpolated time.
149
* This is the main processing method that should be implemented by
150
* the derived classes to perform the interpolation.
151
* @param theta normalized interpolation abscissa within the step
152
* (theta is zero at the previous time step and one at the current time step)
153
* @param oneMinusThetaH time gap between the interpolated time and
155
* @throws DerivativeException this exception is propagated to the caller if the
156
* underlying user function triggers one
158
protected void computeInterpolatedState(double theta,
159
double oneMinusThetaH)
160
throws DerivativeException {
162
if (! vectorsInitialized) {
166
for (int k = 0; k < 7; ++k) {
167
v[k] = new double[interpolatedState.length];
171
// perform the last evaluations if they have not been done yet
174
// compute the interpolation vectors for this time step
175
for (int i = 0; i < interpolatedState.length; ++i) {
176
v[0][i] = h * (b_01 * yDotK[0][i] + b_06 * yDotK[5][i] + b_07 * yDotK[6][i] +
177
b_08 * yDotK[7][i] + b_09 * yDotK[8][i] + b_10 * yDotK[9][i] +
178
b_11 * yDotK[10][i] + b_12 * yDotK[11][i]);
179
v[1][i] = h * yDotK[0][i] - v[0][i];
180
v[2][i] = v[0][i] - v[1][i] - h * yDotK[12][i];
181
for (int k = 0; k < d.length; ++k) {
182
v[k+3][i] = h * (d[k][0] * yDotK[0][i] + d[k][1] * yDotK[5][i] + d[k][2] * yDotK[6][i] +
183
d[k][3] * yDotK[7][i] + d[k][4] * yDotK[8][i] + d[k][5] * yDotK[9][i] +
184
d[k][6] * yDotK[10][i] + d[k][7] * yDotK[11][i] + d[k][8] * yDotK[12][i] +
185
d[k][9] * yDotKLast[0][i] +
186
d[k][10] * yDotKLast[1][i] +
187
d[k][11] * yDotKLast[2][i]);
191
vectorsInitialized = true;
195
double eta = oneMinusThetaH / h;
197
for (int i = 0; i < interpolatedState.length; ++i) {
198
interpolatedState[i] =
199
currentState[i] - eta * (v[0][i] - theta * (v[1][i] +
200
theta * (v[2][i] + eta * (v[3][i] + theta * (v[4][i] +
201
eta * (v[5][i] + theta * (v[6][i])))))));
207
* Really finalize the step.
208
* Perform the last 3 functions evaluations (k14, k15, k16)
209
* @throws DerivativeException this exception is propagated to the caller if the
210
* underlying user function triggers one
212
protected void doFinalize()
213
throws DerivativeException {
215
if (currentState == null) {
216
// we are finalizing an uninitialized instance
221
double[] yTmp = new double[currentState.length];
224
for (int j = 0; j < currentState.length; ++j) {
225
s = k14_01 * yDotK[0][j] + k14_06 * yDotK[5][j] + k14_07 * yDotK[6][j] +
226
k14_08 * yDotK[7][j] + k14_09 * yDotK[8][j] + k14_10 * yDotK[9][j] +
227
k14_11 * yDotK[10][j] + k14_12 * yDotK[11][j] + k14_13 * yDotK[12][j];
228
yTmp[j] = currentState[j] + h * s;
230
equations.computeDerivatives(previousTime + c14 * h, yTmp, yDotKLast[0]);
233
for (int j = 0; j < currentState.length; ++j) {
234
s = k15_01 * yDotK[0][j] + k15_06 * yDotK[5][j] + k15_07 * yDotK[6][j] +
235
k15_08 * yDotK[7][j] + k15_09 * yDotK[8][j] + k15_10 * yDotK[9][j] +
236
k15_11 * yDotK[10][j] + k15_12 * yDotK[11][j] + k15_13 * yDotK[12][j] +
237
k15_14 * yDotKLast[0][j];
238
yTmp[j] = currentState[j] + h * s;
240
equations.computeDerivatives(previousTime + c15 * h, yTmp, yDotKLast[1]);
243
for (int j = 0; j < currentState.length; ++j) {
244
s = k16_01 * yDotK[0][j] + k16_06 * yDotK[5][j] + k16_07 * yDotK[6][j] +
245
k16_08 * yDotK[7][j] + k16_09 * yDotK[8][j] + k16_10 * yDotK[9][j] +
246
k16_11 * yDotK[10][j] + k16_12 * yDotK[11][j] + k16_13 * yDotK[12][j] +
247
k16_14 * yDotKLast[0][j] + k16_15 * yDotKLast[1][j];
248
yTmp[j] = currentState[j] + h * s;
250
equations.computeDerivatives(previousTime + c16 * h, yTmp, yDotKLast[2]);
254
/** Save the state of the instance.
255
* @param out stream where to save the state
256
* @exception IOException in case of write error
258
public void writeExternal(ObjectOutput out)
262
// save the local attributes
264
} catch (DerivativeException e) {
265
throw new IOException(e.getMessage());
267
out.writeInt(currentState.length);
268
for (int i = 0; i < currentState.length; ++i) {
269
out.writeDouble(yDotKLast[0][i]);
270
out.writeDouble(yDotKLast[1][i]);
271
out.writeDouble(yDotKLast[2][i]);
274
// save the state of the base class
275
super.writeExternal(out);
279
/** Read the state of the instance.
280
* @param in stream where to read the state from
281
* @exception IOException in case of read error
283
public void readExternal(ObjectInput in)
286
// read the local attributes
287
yDotKLast = new double[3][];
288
int dimension = in.readInt();
289
yDotKLast[0] = new double[dimension];
290
yDotKLast[1] = new double[dimension];
291
yDotKLast[2] = new double[dimension];
293
for (int i = 0; i < dimension; ++i) {
294
yDotKLast[0][i] = in.readDouble();
295
yDotKLast[1][i] = in.readDouble();
296
yDotKLast[2][i] = in.readDouble();
299
// read the base state
300
super.readExternal(in);
304
/** Last evaluations. */
305
private double[][] yDotKLast;
307
/** Vectors for interpolation. */
308
private double[][] v;
310
/** Initialization indicator for the interpolation vectors. */
311
private boolean vectorsInitialized;
313
/** Propagation weights, element 1. */
314
private static final double b_01 = 104257.0 / 1920240.0;
316
// elements 2 to 5 are zero, so they are neither stored nor used
318
/** Propagation weights, element 6. */
319
private static final double b_06 = 3399327.0 / 763840.0;
321
/** Propagation weights, element 7. */
322
private static final double b_07 = 66578432.0 / 35198415.0;
324
/** Propagation weights, element 8. */
325
private static final double b_08 = -1674902723.0 / 288716400.0;
327
/** Propagation weights, element 9. */
328
private static final double b_09 = 54980371265625.0 / 176692375811392.0;
330
/** Propagation weights, element 10. */
331
private static final double b_10 = -734375.0 / 4826304.0;
333
/** Propagation weights, element 11. */
334
private static final double b_11 = 171414593.0 / 851261400.0;
336
/** Propagation weights, element 12. */
337
private static final double b_12 = 137909.0 / 3084480.0;
339
/** Time step for stage 14 (interpolation only). */
340
private static final double c14 = 1.0 / 10.0;
342
/** Internal weights for stage 14, element 1. */
343
private static final double k14_01 = 13481885573.0 / 240030000000.0 - b_01;
345
// elements 2 to 5 are zero, so they are neither stored nor used
347
/** Internal weights for stage 14, element 6. */
348
private static final double k14_06 = 0.0 - b_06;
350
/** Internal weights for stage 14, element 7. */
351
private static final double k14_07 = 139418837528.0 / 549975234375.0 - b_07;
353
/** Internal weights for stage 14, element 8. */
354
private static final double k14_08 = -11108320068443.0 / 45111937500000.0 - b_08;
356
/** Internal weights for stage 14, element 9. */
357
private static final double k14_09 = -1769651421925959.0 / 14249385146080000.0 - b_09;
359
/** Internal weights for stage 14, element 10. */
360
private static final double k14_10 = 57799439.0 / 377055000.0 - b_10;
362
/** Internal weights for stage 14, element 11. */
363
private static final double k14_11 = 793322643029.0 / 96734250000000.0 - b_11;
365
/** Internal weights for stage 14, element 12. */
366
private static final double k14_12 = 1458939311.0 / 192780000000.0 - b_12;
368
/** Internal weights for stage 14, element 13. */
369
private static final double k14_13 = -4149.0 / 500000.0;
371
/** Time step for stage 15 (interpolation only). */
372
private static final double c15 = 1.0 / 5.0;
375
/** Internal weights for stage 15, element 1. */
376
private static final double k15_01 = 1595561272731.0 / 50120273500000.0 - b_01;
378
// elements 2 to 5 are zero, so they are neither stored nor used
380
/** Internal weights for stage 15, element 6. */
381
private static final double k15_06 = 975183916491.0 / 34457688031250.0 - b_06;
383
/** Internal weights for stage 15, element 7. */
384
private static final double k15_07 = 38492013932672.0 / 718912673015625.0 - b_07;
386
/** Internal weights for stage 15, element 8. */
387
private static final double k15_08 = -1114881286517557.0 / 20298710767500000.0 - b_08;
389
/** Internal weights for stage 15, element 9. */
390
private static final double k15_09 = 0.0 - b_09;
392
/** Internal weights for stage 15, element 10. */
393
private static final double k15_10 = 0.0 - b_10;
395
/** Internal weights for stage 15, element 11. */
396
private static final double k15_11 = -2538710946863.0 / 23431227861250000.0 - b_11;
398
/** Internal weights for stage 15, element 12. */
399
private static final double k15_12 = 8824659001.0 / 23066716781250.0 - b_12;
401
/** Internal weights for stage 15, element 13. */
402
private static final double k15_13 = -11518334563.0 / 33831184612500.0;
404
/** Internal weights for stage 15, element 14. */
405
private static final double k15_14 = 1912306948.0 / 13532473845.0;
407
/** Time step for stage 16 (interpolation only). */
408
private static final double c16 = 7.0 / 9.0;
411
/** Internal weights for stage 16, element 1. */
412
private static final double k16_01 = -13613986967.0 / 31741908048.0 - b_01;
414
// elements 2 to 5 are zero, so they are neither stored nor used
416
/** Internal weights for stage 16, element 6. */
417
private static final double k16_06 = -4755612631.0 / 1012344804.0 - b_06;
419
/** Internal weights for stage 16, element 7. */
420
private static final double k16_07 = 42939257944576.0 / 5588559685701.0 - b_07;
422
/** Internal weights for stage 16, element 8. */
423
private static final double k16_08 = 77881972900277.0 / 19140370552944.0 - b_08;
425
/** Internal weights for stage 16, element 9. */
426
private static final double k16_09 = 22719829234375.0 / 63689648654052.0 - b_09;
428
/** Internal weights for stage 16, element 10. */
429
private static final double k16_10 = 0.0 - b_10;
431
/** Internal weights for stage 16, element 11. */
432
private static final double k16_11 = 0.0 - b_11;
434
/** Internal weights for stage 16, element 12. */
435
private static final double k16_12 = 0.0 - b_12;
437
/** Internal weights for stage 16, element 13. */
438
private static final double k16_13 = -1199007803.0 / 857031517296.0;
440
/** Internal weights for stage 16, element 14. */
441
private static final double k16_14 = 157882067000.0 / 53564469831.0;
443
/** Internal weights for stage 16, element 15. */
444
private static final double k16_15 = -290468882375.0 / 31741908048.0;
446
/** Interpolation weights.
447
* (beware that only the non-null values are in the table)
449
private static final double[][] d = {
451
{ -17751989329.0 / 2106076560.0, 4272954039.0 / 7539864640.0,
452
-118476319744.0 / 38604839385.0, 755123450731.0 / 316657731600.0,
453
3692384461234828125.0 / 1744130441634250432.0, -4612609375.0 / 5293382976.0,
454
2091772278379.0 / 933644586600.0, 2136624137.0 / 3382989120.0,
455
-126493.0 / 1421424.0, 98350000.0 / 5419179.0,
456
-18878125.0 / 2053168.0, -1944542619.0 / 438351368.0},
458
{ 32941697297.0 / 3159114840.0, 456696183123.0 / 1884966160.0,
459
19132610714624.0 / 115814518155.0, -177904688592943.0 / 474986597400.0,
460
-4821139941836765625.0 / 218016305204281304.0, 30702015625.0 / 3970037232.0,
461
-85916079474274.0 / 2800933759800.0, -5919468007.0 / 634310460.0,
462
2479159.0 / 157936.0, -18750000.0 / 602131.0,
463
-19203125.0 / 2053168.0, 15700361463.0 / 438351368.0},
465
{ 12627015655.0 / 631822968.0, -72955222965.0 / 188496616.0,
466
-13145744952320.0 / 69488710893.0, 30084216194513.0 / 56998391688.0,
467
-296858761006640625.0 / 25648977082856624.0, 569140625.0 / 82709109.0,
468
-18684190637.0 / 18672891732.0, 69644045.0 / 89549712.0,
469
-11847025.0 / 4264272.0, -978650000.0 / 16257537.0,
470
519371875.0 / 6159504.0, 5256837225.0 / 438351368.0},
472
{ -450944925.0 / 17550638.0, -14532122925.0 / 94248308.0,
473
-595876966400.0 / 2573655959.0, 188748653015.0 / 527762886.0,
474
2545485458115234375.0 / 27252038150535163.0, -1376953125.0 / 36759604.0,
475
53995596795.0 / 518691437.0, 210311225.0 / 7047894.0,
476
-1718875.0 / 39484.0, 58000000.0 / 602131.0,
477
-1546875.0 / 39484.0, -1262172375.0 / 8429834.0}
481
/** Serializable version identifier */
482
private static final long serialVersionUID = 7152276390558450974L;