1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18 package org.hipparchus.ode.nonstiff;
19
20 import org.hipparchus.CalculusFieldElement;
21 import org.hipparchus.Field;
22 import org.hipparchus.ode.EquationsMapper;
23 import org.hipparchus.ode.FieldEquationsMapper;
24 import org.hipparchus.ode.FieldODEStateAndDerivative;
25 import org.hipparchus.ode.FieldOrdinaryDifferentialEquation;
26 import org.hipparchus.ode.ODEStateAndDerivative;
27 import org.hipparchus.ode.sampling.AbstractFieldODEStateInterpolator;
28 import org.hipparchus.ode.sampling.FieldODEStateInterpolator;
29 import org.hipparchus.ode.sampling.ODEStateInterpolator;
30 import org.hipparchus.util.Binary64Field;
31 import org.hipparchus.util.FastMath;
32 import org.hipparchus.util.MathArrays;
33 import org.junit.jupiter.api.Test;
34
35 import java.lang.reflect.InvocationTargetException;
36
37 import static org.junit.jupiter.api.Assertions.assertEquals;
38 import static org.junit.jupiter.api.Assertions.assertFalse;
39 import static org.junit.jupiter.api.Assertions.assertNotSame;
40 import static org.junit.jupiter.api.Assertions.assertSame;
41 import static org.junit.jupiter.api.Assertions.assertTrue;
42
43 public abstract class FieldODEStateInterpolatorAbstractTest {
44
45 @Test
46 public abstract void interpolationAtBounds();
47
48 protected <T extends CalculusFieldElement<T>> void doInterpolationAtBounds(final Field<T> field, double epsilon) {
49
50 FieldODEStateInterpolator<T> interpolator = setUpInterpolator(field,
51 new SinCos<T>(field),
52 0.0, new double[] { 0.0, 1.0 }, 0.125);
53
54 assertEquals(0.0, interpolator.getPreviousState().getTime().getReal(), 1.0e-15);
55 for (int i = 0; i < 2; ++i) {
56 assertEquals(interpolator.getPreviousState().getPrimaryState()[i].getReal(),
57 interpolator.getInterpolatedState(interpolator.getPreviousState().getTime()).getPrimaryState()[i].getReal(),
58 epsilon);
59 }
60 assertEquals(0.125, interpolator.getCurrentState().getTime().getReal(), 1.0e-15);
61 for (int i = 0; i < 2; ++i) {
62 assertEquals(interpolator.getCurrentState().getPrimaryState()[i].getReal(),
63 interpolator.getInterpolatedState(interpolator.getCurrentState().getTime()).getPrimaryState()[i].getReal(),
64 epsilon);
65 }
66
67 assertFalse(interpolator.isPreviousStateInterpolated());
68 assertFalse(interpolator.isCurrentStateInterpolated());
69 }
70
71 @Test
72 public abstract void interpolationInside();
73
74 protected <T extends CalculusFieldElement<T>> void doInterpolationInside(final Field<T> field,
75 double epsilonSin, double epsilonCos) {
76
77 ReferenceFieldODE<T> sinCos = new SinCos<T>(field);
78 FieldODEStateInterpolator<T> interpolator = setUpInterpolator(field, sinCos,
79 0.0, new double[] { 0.0, 1.0 }, 0.0125);
80
81 int n = 100;
82 double maxErrorSin = 0;
83 double maxErrorCos = 0;
84 for (int i = 0; i <= n; ++i) {
85 T t = interpolator.getPreviousState().getTime().multiply(n - i).
86 add(interpolator.getCurrentState().getTime().multiply(i)).
87 divide(n);
88 FieldODEStateAndDerivative<T> state = interpolator.getInterpolatedState(t);
89 T[] ref = sinCos.theoreticalState(t);
90 maxErrorSin = FastMath.max(maxErrorSin, state.getPrimaryState()[0].subtract(ref[0]).norm());
91 maxErrorCos = FastMath.max(maxErrorCos, state.getPrimaryState()[1].subtract(ref[1]).norm());
92 }
93 assertEquals(0.0, maxErrorSin, epsilonSin);
94 assertEquals(0.0, maxErrorCos, epsilonCos);
95
96 assertFalse(interpolator.isPreviousStateInterpolated());
97 assertFalse(interpolator.isCurrentStateInterpolated());
98 }
99
100 @Test
101 public void restrictPrevious() {
102 doRestrictPrevious(Binary64Field.getInstance(), 1e-15, 1e-15);
103 }
104
105 protected <T extends CalculusFieldElement<T>> void doRestrictPrevious(
106 Field<T> field,
107 double epsilon,
108 double epsilonDot) {
109
110 AbstractFieldODEStateInterpolator<T> original = setUpInterpolator(
111 field, new SinCos<>(field), 0.0, new double[]{0.0, 1.0}, 0.125);
112
113 assertFalse(original.isPreviousStateInterpolated());
114 assertFalse(original.isCurrentStateInterpolated());
115
116 AbstractFieldODEStateInterpolator<T> restricted = original.restrictStep(
117 original.getInterpolatedState(field.getZero().add(1.0 / 32)),
118 original.getCurrentState());
119
120 assertSame(original.getPreviousState(), original.getGlobalPreviousState());
121 assertSame(original.getCurrentState(), original.getGlobalCurrentState());
122 assertSame(original.getGlobalPreviousState(), restricted.getGlobalPreviousState());
123 assertSame(original.getGlobalCurrentState(), restricted.getGlobalCurrentState());
124 assertNotSame(restricted.getPreviousState(), restricted.getGlobalPreviousState());
125 assertSame(restricted.getCurrentState(), restricted.getGlobalCurrentState());
126 assertEquals(1.0 / 32, restricted.getPreviousState().getTime().getReal(), 1.0e-15);
127 assertTrue(restricted.isPreviousStateInterpolated());
128 assertFalse(restricted.isCurrentStateInterpolated());
129
130 checkRestricted(original, restricted, epsilon, epsilonDot);
131
132 }
133
134 @Test
135 public void restrictCurrent() {
136 doRestrictCurrent(Binary64Field.getInstance(), 1e-15, 1e-15);
137 }
138
139 protected <T extends CalculusFieldElement<T>> void doRestrictCurrent(Field<T> field,
140 double epsilon,
141 double epsilonDot) {
142
143 AbstractFieldODEStateInterpolator<T> original = setUpInterpolator(
144 field, new SinCos<>(field), 0.0, new double[]{0.0, 1.0}, 0.125);
145
146 assertFalse(original.isPreviousStateInterpolated());
147 assertFalse(original.isCurrentStateInterpolated());
148
149 AbstractFieldODEStateInterpolator<T> restricted = original.restrictStep(
150 original.getPreviousState(),
151 original.getInterpolatedState(field.getZero().add(3.0 / 32)));
152
153 assertSame(original.getPreviousState(), original.getGlobalPreviousState());
154 assertSame(original.getCurrentState(), original.getGlobalCurrentState());
155 assertSame(original.getGlobalPreviousState(), restricted.getGlobalPreviousState());
156 assertSame(original.getGlobalCurrentState(), restricted.getGlobalCurrentState());
157 assertSame(restricted.getPreviousState(), restricted.getGlobalPreviousState());
158 assertNotSame(restricted.getCurrentState(), restricted.getGlobalCurrentState());
159 assertEquals(3.0 / 32, restricted.getCurrentState().getTime().getReal(), 1.0e-15);
160 assertFalse(restricted.isPreviousStateInterpolated());
161 assertTrue(restricted.isCurrentStateInterpolated());
162
163 checkRestricted(original, restricted, epsilon, epsilonDot);
164
165 }
166
167 @Test
168 public void restrictBothEnds() {
169 doRestrictBothEnds(Binary64Field.getInstance(), 1e-15, 1e-15);
170 }
171
172 protected <T extends CalculusFieldElement<T>> void doRestrictBothEnds(Field<T> field,
173 double epsilon,
174 double epsilonDot) {
175
176 AbstractFieldODEStateInterpolator<T> original = setUpInterpolator(
177 field, new SinCos<>(field), 0.0, new double[]{0.0, 1.0}, 0.125);
178
179 assertFalse(original.isPreviousStateInterpolated());
180 assertFalse(original.isCurrentStateInterpolated());
181
182 AbstractFieldODEStateInterpolator<T> restricted = original.restrictStep(
183 original.getInterpolatedState(field.getZero().add(1.0 / 32)),
184 original.getInterpolatedState(field.getZero().add(3.0 / 32)));
185
186 assertSame(original.getPreviousState(), original.getGlobalPreviousState());
187 assertSame(original.getCurrentState(), original.getGlobalCurrentState());
188 assertSame(original.getGlobalPreviousState(), restricted.getGlobalPreviousState());
189 assertSame(original.getGlobalCurrentState(), restricted.getGlobalCurrentState());
190 assertNotSame(restricted.getPreviousState(), restricted.getGlobalPreviousState());
191 assertNotSame(restricted.getCurrentState(), restricted.getGlobalCurrentState());
192 assertEquals(1.0 / 32, restricted.getPreviousState().getTime().getReal(), 1.0e-15);
193 assertEquals(3.0 / 32, restricted.getCurrentState().getTime().getReal(), 1.0e-15);
194 assertTrue(restricted.isPreviousStateInterpolated());
195 assertTrue(restricted.isCurrentStateInterpolated());
196
197 checkRestricted(original, restricted, epsilon, epsilonDot);
198
199 }
200
201 @Test
202 public void degenerateInterpolation() {
203 doDegenerateInterpolation(Binary64Field.getInstance());
204 }
205
206 protected <T extends CalculusFieldElement<T>> void doDegenerateInterpolation(Field<T> field) {
207 AbstractFieldODEStateInterpolator<T> interpolator = setUpInterpolator(
208 field, new SinCos<>(field), 0.0, new double[] { 0.0, 1.0 }, 0.0);
209 FieldODEStateAndDerivative<T> interpolatedState = interpolator.getInterpolatedState(field.getZero());
210 assertEquals(0.0, interpolatedState.getTime().getReal(), 0.0);
211 assertEquals(0.0, interpolatedState.getPrimaryState()[0].getReal(), 0.0);
212 assertEquals(1.0, interpolatedState.getPrimaryState()[1].getReal(), 0.0);
213 assertEquals(1.0, interpolatedState.getPrimaryDerivative()[0].getReal(), 0.0);
214 assertEquals(0.0, interpolatedState.getPrimaryDerivative()[1].getReal(), 0.0);
215 }
216
217 private <T extends CalculusFieldElement<T>> void checkRestricted(
218 AbstractFieldODEStateInterpolator<T> original,
219 AbstractFieldODEStateInterpolator<T> restricted,
220 double epsilon,
221 double epsilonDot) {
222
223 for (T t = restricted.getPreviousState().getTime();
224 t.getReal() <= restricted.getCurrentState().getTime().getReal();
225 t = t.add(1.0 / 256)) {
226 FieldODEStateAndDerivative<T> originalInterpolated = original.getInterpolatedState(t);
227 FieldODEStateAndDerivative<T> restrictedInterpolated = restricted.getInterpolatedState(t);
228 assertEquals(t.getReal(), originalInterpolated.getTime().getReal(), 1.0e-15);
229 assertEquals(t.getReal(), restrictedInterpolated.getTime().getReal(), 1.0e-15);
230 assertEquals(originalInterpolated.getPrimaryState()[0].getReal(),
231 restrictedInterpolated.getPrimaryState()[0].getReal(),
232 epsilon);
233 assertEquals(originalInterpolated.getPrimaryState()[1].getReal(),
234 restrictedInterpolated.getPrimaryState()[1].getReal(),
235 epsilon);
236 assertEquals(originalInterpolated.getPrimaryDerivative()[0].getReal(),
237 restrictedInterpolated.getPrimaryDerivative()[0].getReal(),
238 epsilonDot);
239 assertEquals(originalInterpolated.getPrimaryDerivative()[1].getReal(),
240 restrictedInterpolated.getPrimaryDerivative()[1].getReal(),
241 epsilonDot);
242 }
243
244 }
245
246 @Test
247 public abstract void nonFieldInterpolatorConsistency();
248
249 protected <T extends CalculusFieldElement<T>> void doNonFieldInterpolatorConsistency(final Field<T> field,
250 double epsilonSin, double epsilonCos,
251 double epsilonSinDot, double epsilonCosDot) {
252
253 ReferenceFieldODE<T> eqn = new SinCos<T>(field);
254 FieldODEStateInterpolator<T> fieldInterpolator =
255 setUpInterpolator(field, eqn, 0.0, new double[] { 0.0, 1.0 }, 0.125);
256 ODEStateInterpolator regularInterpolator = convertInterpolator(fieldInterpolator, eqn);
257
258 int n = 100;
259 double maxErrorSin = 0;
260 double maxErrorCos = 0;
261 double maxErrorSinDot = 0;
262 double maxErrorCosDot = 0;
263 for (int i = 0; i <= n; ++i) {
264
265 T t = fieldInterpolator.getPreviousState().getTime().multiply(n - i).
266 add(fieldInterpolator.getCurrentState().getTime().multiply(i)).
267 divide(n);
268
269 FieldODEStateAndDerivative<T> state = fieldInterpolator.getInterpolatedState(t);
270 T[] fieldY = state.getPrimaryState();
271 T[] fieldYDot = state.getPrimaryDerivative();
272
273 ODEStateAndDerivative regularState = regularInterpolator.getInterpolatedState(t.getReal());
274 double[] regularY = regularState.getPrimaryState();
275 double[] regularYDot = regularState.getPrimaryDerivative();
276
277 maxErrorSin = FastMath.max(maxErrorSin, fieldY[0].subtract(regularY[0]).norm());
278 maxErrorCos = FastMath.max(maxErrorCos, fieldY[1].subtract(regularY[1]).norm());
279 maxErrorSinDot = FastMath.max(maxErrorSinDot, fieldYDot[0].subtract(regularYDot[0]).norm());
280 maxErrorCosDot = FastMath.max(maxErrorCosDot, fieldYDot[1].subtract(regularYDot[1]).norm());
281
282 }
283 assertEquals(0.0, maxErrorSin, epsilonSin);
284 assertEquals(0.0, maxErrorCos, epsilonCos);
285 assertEquals(0.0, maxErrorSinDot, epsilonSinDot);
286 assertEquals(0.0, maxErrorCosDot, epsilonCosDot);
287
288 }
289
290 public interface ReferenceFieldODE<T extends CalculusFieldElement<T>> extends FieldOrdinaryDifferentialEquation<T> {
291 T[] theoreticalState(T t);
292 }
293
294 protected abstract <T extends CalculusFieldElement<T>>
295 AbstractFieldODEStateInterpolator<T> setUpInterpolator(final Field<T> field,
296 final ReferenceFieldODE<T> eqn,
297 final double t0,
298 final double[] y0,
299 final double t1);
300
301 protected abstract <T extends CalculusFieldElement<T>>
302 ODEStateInterpolator convertInterpolator(final FieldODEStateInterpolator<T> fieldInterpolator,
303 final FieldOrdinaryDifferentialEquation<T> eqn);
304
305 protected <T extends CalculusFieldElement<T>>
306 ODEStateAndDerivative convertODEStateAndDerivative(final FieldODEStateAndDerivative<T> s) {
307 final double[][] secondaryStates;
308 final double[][] secondaryDerivatives;
309 if (s.getNumberOfSecondaryStates() == 0) {
310 secondaryStates = null;
311 secondaryDerivatives = null;
312 } else {
313 secondaryStates = new double[s.getNumberOfSecondaryStates()][];
314 secondaryDerivatives = new double[s.getNumberOfSecondaryStates()][];
315 for (int i = 0; i < secondaryStates.length; ++i) {
316 secondaryStates[i] = convertArray(s.getSecondaryState(i));
317 secondaryDerivatives[i] = convertArray(s.getSecondaryDerivative(i));
318 }
319 }
320 return new ODEStateAndDerivative(s.getTime().getReal(),
321 convertArray(s.getPrimaryState()),
322 convertArray(s.getPrimaryDerivative()),
323 secondaryStates,
324 secondaryDerivatives);
325 }
326
327 protected <T extends CalculusFieldElement<T>> double[][] convertArray(final T[][] fieldArray) {
328 if (fieldArray == null) {
329 return null;
330 }
331 double[][] array = new double[fieldArray.length][];
332 for (int i = 0; i < array.length; ++i) {
333 array[i] = convertArray(fieldArray[i]);
334 }
335 return array;
336 }
337
338 protected <T extends CalculusFieldElement<T>> double[] convertArray(final T[] fieldArray) {
339 if (fieldArray == null) {
340 return null;
341 }
342 double[] array = new double[fieldArray.length];
343 for (int i = 0; i < array.length; ++i) {
344 array[i] = fieldArray[i].getReal();
345 }
346 return array;
347 }
348
349 protected <T extends CalculusFieldElement<T>>
350 EquationsMapper convertMapper(final FieldEquationsMapper<T> fieldmapper)
351 throws NoSuchMethodException, SecurityException, NoSuchFieldException,
352 IllegalArgumentException, IllegalAccessException,
353 InstantiationException, InvocationTargetException {
354 java.lang.reflect.Field fStart = FieldEquationsMapper.class.getDeclaredField("start");
355 fStart.setAccessible(true);
356 int[] start = (int[]) fStart.get(fieldmapper);
357
358 java.lang.reflect.Constructor<EquationsMapper> regularMapperConstructor =
359 EquationsMapper.class.getDeclaredConstructor(EquationsMapper.class,
360 Integer.TYPE);
361 regularMapperConstructor.setAccessible(true);
362 EquationsMapper regularMapper = regularMapperConstructor.newInstance(null, start[1]);
363 for (int k = 1; k < fieldmapper.getNumberOfEquations(); ++k) {
364 regularMapper = regularMapperConstructor.newInstance(regularMapper,
365 start[k + 1] - start[k]);
366 }
367 return regularMapper;
368 }
369
370 private static class SinCos<T extends CalculusFieldElement<T>> implements ReferenceFieldODE<T> {
371 private final Field<T> field;
372 protected SinCos(final Field<T> field) {
373 this.field = field;
374 }
375 public int getDimension() {
376 return 2;
377 }
378 public void init(final T t0, final T[] y0, final T finalTime) {
379 }
380 public T[] computeDerivatives(final T t, final T[] y) {
381 T[] yDot = MathArrays.buildArray(field, 2);
382 yDot[0] = y[1];
383 yDot[1] = y[0].negate();
384 return yDot;
385 }
386 public T[] theoreticalState(final T t) {
387 T[] state = MathArrays.buildArray(field, 2);
388 state[0] = t.sin();
389 state[1] = t.cos();
390 return state;
391 }
392 }
393
394 }