1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22 package org.hipparchus.ode.sampling;
23
24
25 import org.hipparchus.CalculusFieldElement;
26 import org.hipparchus.exception.MathIllegalArgumentException;
27 import org.hipparchus.exception.MathIllegalStateException;
28 import org.hipparchus.ode.FieldExpandableODE;
29 import org.hipparchus.ode.FieldODEIntegrator;
30 import org.hipparchus.ode.ODEIntegrator;
31 import org.hipparchus.ode.TestFieldProblemAbstract;
32 import org.hipparchus.ode.TestProblemAbstract;
33 import org.hipparchus.util.FastMath;
34
35 import static org.junit.jupiter.api.Assertions.assertEquals;
36
37 public class StepInterpolatorTestUtils {
38
39 public static void checkDerivativesConsistency(final ODEIntegrator integrator,
40 final TestProblemAbstract problem,
41 final double finiteDifferencesRatio,
42 final double threshold)
43 throws MathIllegalArgumentException, MathIllegalStateException {
44 integrator.addStepHandler(new ODEStepHandler() {
45
46 public void handleStep(ODEStateInterpolator interpolator) {
47
48 final double dt = interpolator.getCurrentState().getTime() - interpolator.getPreviousState().getTime();
49 final double h = finiteDifferencesRatio * dt;
50 final double t = interpolator.getCurrentState().getTime() - 0.3 * dt;
51
52 if (FastMath.abs(h) < 10 * FastMath.ulp(t)) {
53 return;
54 }
55
56 final double[] yM4h = interpolator.getInterpolatedState(t - 4 * h).getPrimaryState();
57 final double[] yM3h = interpolator.getInterpolatedState(t - 3 * h).getPrimaryState();
58 final double[] yM2h = interpolator.getInterpolatedState(t - 2 * h).getPrimaryState();
59 final double[] yM1h = interpolator.getInterpolatedState(t - h).getPrimaryState();
60 final double[] yP1h = interpolator.getInterpolatedState(t + h).getPrimaryState();
61 final double[] yP2h = interpolator.getInterpolatedState(t + 2 * h).getPrimaryState();
62 final double[] yP3h = interpolator.getInterpolatedState(t + 3 * h).getPrimaryState();
63 final double[] yP4h = interpolator.getInterpolatedState(t + 4 * h).getPrimaryState();
64
65 final double[] yDot = interpolator.getInterpolatedState(t).getPrimaryDerivative();
66
67 for (int i = 0; i < yDot.length; ++i) {
68 final double approYDot = ( -3 * (yP4h[i] - yM4h[i]) +
69 32 * (yP3h[i] - yM3h[i]) +
70 -168 * (yP2h[i] - yM2h[i]) +
71 672 * (yP1h[i] - yM1h[i])) / (840 * h);
72 assertEquals(approYDot, yDot[i], threshold, "" + (approYDot - yDot[i]));
73 }
74
75 }
76
77 });
78
79 integrator.integrate(problem, problem.getInitialState(), problem.getFinalTime());
80
81 }
82
83 public static <T extends CalculusFieldElement<T>> void checkDerivativesConsistency(final FieldODEIntegrator<T> integrator,
84 final TestFieldProblemAbstract<T> problem,
85 final double threshold) {
86 integrator.addStepHandler(new FieldODEStepHandler<T>() {
87
88 public void handleStep(FieldODEStateInterpolator<T> interpolator) {
89
90 final T h = interpolator.getCurrentState().getTime().subtract(interpolator.getPreviousState().getTime()).multiply(0.001);
91 final T t = interpolator.getCurrentState().getTime().subtract(h.multiply(300));
92
93 if (h.abs().subtract(FastMath.ulp(t.getReal()) * 10).getReal() < 0) {
94 return;
95 }
96
97 final T[] yM4h = interpolator.getInterpolatedState(t.add(h.multiply(-4))).getPrimaryState();
98 final T[] yM3h = interpolator.getInterpolatedState(t.add(h.multiply(-3))).getPrimaryState();
99 final T[] yM2h = interpolator.getInterpolatedState(t.add(h.multiply(-2))).getPrimaryState();
100 final T[] yM1h = interpolator.getInterpolatedState(t.add(h.multiply(-1))).getPrimaryState();
101 final T[] yP1h = interpolator.getInterpolatedState(t.add(h.multiply( 1))).getPrimaryState();
102 final T[] yP2h = interpolator.getInterpolatedState(t.add(h.multiply( 2))).getPrimaryState();
103 final T[] yP3h = interpolator.getInterpolatedState(t.add(h.multiply( 3))).getPrimaryState();
104 final T[] yP4h = interpolator.getInterpolatedState(t.add(h.multiply( 4))).getPrimaryState();
105
106 final T[] yDot = interpolator.getInterpolatedState(t).getPrimaryDerivative();
107
108 for (int i = 0; i < yDot.length; ++i) {
109 final T approYDot = yP4h[i].subtract(yM4h[i]).multiply( -3).
110 add(yP3h[i].subtract(yM3h[i]).multiply( 32)).
111 add(yP2h[i].subtract(yM2h[i]).multiply(-168)).
112 add(yP1h[i].subtract(yM1h[i]).multiply( 672)).
113 divide(h.multiply(840));
114 assertEquals(approYDot.getReal(), yDot[i].getReal(), threshold);
115 }
116
117 }
118
119 });
120
121 integrator.integrate(new FieldExpandableODE<T>(problem), problem.getInitialState(), problem.getFinalTime());
122
123 }
124 }
125