1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18 package org.apache.commons.math4.legacy.ode.sampling;
19
20 import java.io.IOException;
21 import java.io.ObjectInput;
22 import java.io.ObjectOutput;
23 import java.util.Arrays;
24
25 import org.apache.commons.math4.legacy.exception.MaxCountExceededException;
26 import org.apache.commons.math4.legacy.linear.Array2DRowRealMatrix;
27 import org.apache.commons.math4.legacy.ode.EquationsMapper;
28 import org.apache.commons.math4.core.jdkmath.JdkMath;
29
30
31
32
33
34
35
36
37
38
39
40
41 public class NordsieckStepInterpolator extends AbstractStepInterpolator {
42
43
44 private static final long serialVersionUID = -7179861704951334960L;
45
46
47 protected double[] stateVariation;
48
49
50 private double scalingH;
51
52
53
54
55
56
57
58 private double referenceTime;
59
60
61 private double[] scaled;
62
63
64 private Array2DRowRealMatrix nordsieck;
65
66
67
68
69
70
71
72
73 public NordsieckStepInterpolator() {
74 }
75
76
77
78
79
80
81 public NordsieckStepInterpolator(final NordsieckStepInterpolator interpolator) {
82 super(interpolator);
83 scalingH = interpolator.scalingH;
84 referenceTime = interpolator.referenceTime;
85 if (interpolator.scaled != null) {
86 scaled = interpolator.scaled.clone();
87 }
88 if (interpolator.nordsieck != null) {
89 nordsieck = new Array2DRowRealMatrix(interpolator.nordsieck.getDataRef(), true);
90 }
91 if (interpolator.stateVariation != null) {
92 stateVariation = interpolator.stateVariation.clone();
93 }
94 }
95
96
97 @Override
98 protected StepInterpolator doCopy() {
99 return new NordsieckStepInterpolator(this);
100 }
101
102
103
104
105
106
107
108
109
110
111 @Override
112 public void reinitialize(final double[] y, final boolean forward,
113 final EquationsMapper primaryMapper,
114 final EquationsMapper[] secondaryMappers) {
115 super.reinitialize(y, forward, primaryMapper, secondaryMappers);
116 stateVariation = new double[y.length];
117 }
118
119
120
121
122
123
124
125
126
127
128
129 public void reinitialize(final double time, final double stepSize,
130 final double[] scaledDerivative,
131 final Array2DRowRealMatrix nordsieckVector) {
132 this.referenceTime = time;
133 this.scalingH = stepSize;
134 this.scaled = scaledDerivative;
135 this.nordsieck = nordsieckVector;
136
137
138 setInterpolatedTime(getInterpolatedTime());
139 }
140
141
142
143
144
145
146 public void rescale(final double stepSize) {
147
148 final double ratio = stepSize / scalingH;
149 for (int i = 0; i < scaled.length; ++i) {
150 scaled[i] *= ratio;
151 }
152
153 final double[][] nData = nordsieck.getDataRef();
154 double power = ratio;
155 for (int i = 0; i < nData.length; ++i) {
156 power *= ratio;
157 final double[] nDataI = nData[i];
158 for (int j = 0; j < nDataI.length; ++j) {
159 nDataI[j] *= power;
160 }
161 }
162
163 scalingH = stepSize;
164 }
165
166
167
168
169
170
171
172
173
174
175
176
177
178 public double[] getInterpolatedStateVariation() throws MaxCountExceededException {
179
180
181 getInterpolatedState();
182 return stateVariation;
183 }
184
185
186 @Override
187 protected void computeInterpolatedStateAndDerivatives(final double theta, final double oneMinusThetaH) {
188
189 final double x = interpolatedTime - referenceTime;
190 final double normalizedAbscissa = x / scalingH;
191
192 Arrays.fill(stateVariation, 0.0);
193 Arrays.fill(interpolatedDerivatives, 0.0);
194
195
196
197 final double[][] nData = nordsieck.getDataRef();
198 for (int i = nData.length - 1; i >= 0; --i) {
199 final int order = i + 2;
200 final double[] nDataI = nData[i];
201 final double power = JdkMath.pow(normalizedAbscissa, order);
202 for (int j = 0; j < nDataI.length; ++j) {
203 final double d = nDataI[j] * power;
204 stateVariation[j] += d;
205 interpolatedDerivatives[j] += order * d;
206 }
207 }
208
209 for (int j = 0; j < currentState.length; ++j) {
210 stateVariation[j] += scaled[j] * normalizedAbscissa;
211 interpolatedState[j] = currentState[j] + stateVariation[j];
212 interpolatedDerivatives[j] =
213 (interpolatedDerivatives[j] + scaled[j] * normalizedAbscissa) / x;
214 }
215 }
216
217
218 @Override
219 public void writeExternal(final ObjectOutput out)
220 throws IOException {
221
222
223 writeBaseExternal(out);
224
225
226 out.writeDouble(scalingH);
227 out.writeDouble(referenceTime);
228
229 final int n = (currentState == null) ? -1 : currentState.length;
230 if (scaled == null) {
231 out.writeBoolean(false);
232 } else {
233 out.writeBoolean(true);
234 for (int j = 0; j < n; ++j) {
235 out.writeDouble(scaled[j]);
236 }
237 }
238
239 if (nordsieck == null) {
240 out.writeBoolean(false);
241 } else {
242 out.writeBoolean(true);
243 out.writeObject(nordsieck);
244 }
245
246
247 }
248
249
250 @Override
251 public void readExternal(final ObjectInput in)
252 throws IOException, ClassNotFoundException {
253
254
255 final double t = readBaseExternal(in);
256
257
258 scalingH = in.readDouble();
259 referenceTime = in.readDouble();
260
261 final int n = (currentState == null) ? -1 : currentState.length;
262 final boolean hasScaled = in.readBoolean();
263 if (hasScaled) {
264 scaled = new double[n];
265 for (int j = 0; j < n; ++j) {
266 scaled[j] = in.readDouble();
267 }
268 } else {
269 scaled = null;
270 }
271
272 final boolean hasNordsieck = in.readBoolean();
273 if (hasNordsieck) {
274 nordsieck = (Array2DRowRealMatrix) in.readObject();
275 } else {
276 nordsieck = null;
277 }
278
279 if (hasScaled && hasNordsieck) {
280
281 stateVariation = new double[n];
282 setInterpolatedTime(t);
283 } else {
284 stateVariation = null;
285 }
286 }
287 }