1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17 package org.apache.commons.math4.legacy.ode.sampling;
18
19
20 import org.apache.commons.math4.legacy.core.RealFieldElement;
21 import org.apache.commons.math4.legacy.exception.DimensionMismatchException;
22 import org.apache.commons.math4.legacy.exception.MaxCountExceededException;
23 import org.apache.commons.math4.legacy.exception.NoBracketingException;
24 import org.apache.commons.math4.legacy.exception.NumberIsTooSmallException;
25 import org.apache.commons.math4.legacy.ode.FieldExpandableODE;
26 import org.apache.commons.math4.legacy.ode.FirstOrderFieldIntegrator;
27 import org.apache.commons.math4.legacy.ode.FieldODEStateAndDerivative;
28 import org.apache.commons.math4.legacy.ode.FirstOrderIntegrator;
29 import org.apache.commons.math4.legacy.ode.TestFieldProblemAbstract;
30 import org.apache.commons.math4.legacy.ode.TestProblemAbstract;
31 import org.apache.commons.math4.core.jdkmath.JdkMath;
32 import org.junit.Assert;
33
34 public final class StepInterpolatorTestUtils {
35
36
37 private StepInterpolatorTestUtils() {}
38
39 public static void checkDerivativesConsistency(final FirstOrderIntegrator integrator,
40 final TestProblemAbstract problem,
41 final double finiteDifferencesRatio,
42 final double threshold)
43 throws DimensionMismatchException, NumberIsTooSmallException,
44 MaxCountExceededException, NoBracketingException {
45 integrator.addStepHandler(new StepHandler() {
46
47 @Override
48 public void handleStep(StepInterpolator interpolator, boolean isLast)
49 throws MaxCountExceededException {
50
51 final double dt = interpolator.getCurrentTime() - interpolator.getPreviousTime();
52 final double h = finiteDifferencesRatio * dt;
53 final double t = interpolator.getCurrentTime() - 0.3 * dt;
54
55 if (JdkMath.abs(h) < 10 * JdkMath.ulp(t)) {
56 return;
57 }
58
59 interpolator.setInterpolatedTime(t - 4 * h);
60 final double[] yM4h = interpolator.getInterpolatedState().clone();
61 interpolator.setInterpolatedTime(t - 3 * h);
62 final double[] yM3h = interpolator.getInterpolatedState().clone();
63 interpolator.setInterpolatedTime(t - 2 * h);
64 final double[] yM2h = interpolator.getInterpolatedState().clone();
65 interpolator.setInterpolatedTime(t - h);
66 final double[] yM1h = interpolator.getInterpolatedState().clone();
67 interpolator.setInterpolatedTime(t + h);
68 final double[] yP1h = interpolator.getInterpolatedState().clone();
69 interpolator.setInterpolatedTime(t + 2 * h);
70 final double[] yP2h = interpolator.getInterpolatedState().clone();
71 interpolator.setInterpolatedTime(t + 3 * h);
72 final double[] yP3h = interpolator.getInterpolatedState().clone();
73 interpolator.setInterpolatedTime(t + 4 * h);
74 final double[] yP4h = interpolator.getInterpolatedState().clone();
75
76 interpolator.setInterpolatedTime(t);
77 final double[] yDot = interpolator.getInterpolatedDerivatives();
78
79 for (int i = 0; i < yDot.length; ++i) {
80 final double approYDot = ( -3 * (yP4h[i] - yM4h[i]) +
81 32 * (yP3h[i] - yM3h[i]) +
82 -168 * (yP2h[i] - yM2h[i]) +
83 672 * (yP1h[i] - yM1h[i])) / (840 * h);
84 Assert.assertEquals("" + (approYDot - yDot[i]), approYDot, yDot[i], threshold);
85 }
86 }
87
88 @Override
89 public void init(double t0, double[] y0, double t) {
90 }
91 });
92
93 integrator.integrate(problem,
94 problem.getInitialTime(), problem.getInitialState(),
95 problem.getFinalTime(), new double[problem.getDimension()]);
96 }
97
98 public static <T extends RealFieldElement<T>> void checkDerivativesConsistency(final FirstOrderFieldIntegrator<T> integrator,
99 final TestFieldProblemAbstract<T> problem,
100 final double threshold) {
101 integrator.addStepHandler(new FieldStepHandler<T>() {
102
103 @Override
104 public void handleStep(FieldStepInterpolator<T> interpolator, boolean isLast)
105 throws MaxCountExceededException {
106
107 final T h = interpolator.getCurrentState().getTime().subtract(interpolator.getPreviousState().getTime()).multiply(0.001);
108 final T t = interpolator.getCurrentState().getTime().subtract(h.multiply(300));
109
110 if (h.abs().subtract(JdkMath.ulp(t.getReal()) * 10).getReal() < 0) {
111 return;
112 }
113
114 final T[] yM4h = interpolator.getInterpolatedState(t.add(h.multiply(-4))).getState();
115 final T[] yM3h = interpolator.getInterpolatedState(t.add(h.multiply(-3))).getState();
116 final T[] yM2h = interpolator.getInterpolatedState(t.add(h.multiply(-2))).getState();
117 final T[] yM1h = interpolator.getInterpolatedState(t.add(h.multiply(-1))).getState();
118 final T[] yP1h = interpolator.getInterpolatedState(t.add(h.multiply( 1))).getState();
119 final T[] yP2h = interpolator.getInterpolatedState(t.add(h.multiply( 2))).getState();
120 final T[] yP3h = interpolator.getInterpolatedState(t.add(h.multiply( 3))).getState();
121 final T[] yP4h = interpolator.getInterpolatedState(t.add(h.multiply( 4))).getState();
122
123 final T[] yDot = interpolator.getInterpolatedState(t).getDerivative();
124
125 for (int i = 0; i < yDot.length; ++i) {
126 final T approYDot = yP4h[i].subtract(yM4h[i]).multiply( -3).
127 add(yP3h[i].subtract(yM3h[i]).multiply( 32)).
128 add(yP2h[i].subtract(yM2h[i]).multiply(-168)).
129 add(yP1h[i].subtract(yM1h[i]).multiply( 672)).
130 divide(h.multiply(840));
131 Assert.assertEquals(approYDot.getReal(), yDot[i].getReal(), threshold);
132 }
133 }
134
135 @Override
136 public void init(FieldODEStateAndDerivative<T> state0, T t) {
137 }
138 });
139
140 integrator.integrate(new FieldExpandableODE<>(problem), problem.getInitialState(), problem.getFinalTime());
141 }
142 }
143