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