1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18 package org.apache.commons.math.ode.nonstiff;
19
20 import org.apache.commons.math.ode.AbstractIntegrator;
21 import org.apache.commons.math.ode.DerivativeException;
22 import org.apache.commons.math.ode.sampling.StepInterpolator;
23
24
25
26
27
28
29
30
31
32
33
34 class DormandPrince54StepInterpolator
35 extends RungeKuttaStepInterpolator {
36
37
38
39
40
41
42
43
44
45
46 public DormandPrince54StepInterpolator() {
47 super();
48 v1 = null;
49 v2 = null;
50 v3 = null;
51 v4 = null;
52 vectorsInitialized = false;
53 }
54
55
56
57
58
59
60 public DormandPrince54StepInterpolator(final DormandPrince54StepInterpolator interpolator) {
61
62 super(interpolator);
63
64 if (interpolator.v1 == null) {
65
66 v1 = null;
67 v2 = null;
68 v3 = null;
69 v4 = null;
70 vectorsInitialized = false;
71
72 } else {
73
74 v1 = interpolator.v1.clone();
75 v2 = interpolator.v2.clone();
76 v3 = interpolator.v3.clone();
77 v4 = interpolator.v4.clone();
78 vectorsInitialized = interpolator.vectorsInitialized;
79
80 }
81
82 }
83
84
85 @Override
86 protected StepInterpolator doCopy() {
87 return new DormandPrince54StepInterpolator(this);
88 }
89
90
91
92 @Override
93 public void reinitialize(final AbstractIntegrator integrator,
94 final double[] y, final double[][] yDotK, final boolean forward) {
95 super.reinitialize(integrator, y, yDotK, forward);
96 v1 = null;
97 v2 = null;
98 v3 = null;
99 v4 = null;
100 vectorsInitialized = false;
101 }
102
103
104 @Override
105 public void storeTime(final double t) {
106 super.storeTime(t);
107 vectorsInitialized = false;
108 }
109
110
111 @Override
112 protected void computeInterpolatedStateAndDerivatives(final double theta,
113 final double oneMinusThetaH)
114 throws DerivativeException {
115
116 if (! vectorsInitialized) {
117
118 if (v1 == null) {
119 v1 = new double[interpolatedState.length];
120 v2 = new double[interpolatedState.length];
121 v3 = new double[interpolatedState.length];
122 v4 = new double[interpolatedState.length];
123 }
124
125
126
127
128 for (int i = 0; i < interpolatedState.length; ++i) {
129 final double yDot0 = yDotK[0][i];
130 final double yDot2 = yDotK[2][i];
131 final double yDot3 = yDotK[3][i];
132 final double yDot4 = yDotK[4][i];
133 final double yDot5 = yDotK[5][i];
134 final double yDot6 = yDotK[6][i];
135 v1[i] = a70 * yDot0 + a72 * yDot2 + a73 * yDot3 + a74 * yDot4 + a75 * yDot5;
136 v2[i] = yDot0 - v1[i];
137 v3[i] = v1[i] - v2[i] - yDot6;
138 v4[i] = d0 * yDot0 + d2 * yDot2 + d3 * yDot3 + d4 * yDot4 + d5 * yDot5 + d6 * yDot6;
139 }
140
141 vectorsInitialized = true;
142
143 }
144
145
146 final double eta = 1 - theta;
147 final double twoTheta = 2 * theta;
148 final double dot2 = 1 - twoTheta;
149 final double dot3 = theta * (2 - 3 * theta);
150 final double dot4 = twoTheta * (1 + theta * (twoTheta - 3));
151 for (int i = 0; i < interpolatedState.length; ++i) {
152 interpolatedState[i] =
153 currentState[i] - oneMinusThetaH * (v1[i] - theta * (v2[i] + theta * (v3[i] + eta * v4[i])));
154 interpolatedDerivatives[i] = v1[i] + dot2 * v2[i] + dot3 * v3[i] + dot4 * v4[i];
155 }
156
157 }
158
159
160 private double[] v1;
161
162
163 private double[] v2;
164
165
166 private double[] v3;
167
168
169 private double[] v4;
170
171
172 private boolean vectorsInitialized;
173
174
175 private static final double a70 = 35.0 / 384.0;
176
177
178
179
180 private static final double a72 = 500.0 / 1113.0;
181
182
183 private static final double a73 = 125.0 / 192.0;
184
185
186 private static final double a74 = -2187.0 / 6784.0;
187
188
189 private static final double a75 = 11.0 / 84.0;
190
191
192 private static final double d0 = -12715105075.0 / 11282082432.0;
193
194
195
196
197 private static final double d2 = 87487479700.0 / 32700410799.0;
198
199
200 private static final double d3 = -10690763975.0 / 1880347072.0;
201
202
203 private static final double d4 = 701980252875.0 / 199316789632.0;
204
205
206 private static final double d5 = -1453857185.0 / 822651844.0;
207
208
209 private static final double d6 = 69997945.0 / 29380423.0;
210
211
212 private static final long serialVersionUID = 4104157279605906956L;
213
214 }