1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23 package org.hipparchus.ode;
24
25 import org.hipparchus.CalculusFieldElement;
26 import org.hipparchus.Field;
27 import org.hipparchus.exception.MathIllegalArgumentException;
28 import org.hipparchus.ode.nonstiff.DormandPrince54FieldIntegrator;
29 import org.hipparchus.ode.nonstiff.DormandPrince853FieldIntegrator;
30 import org.hipparchus.ode.sampling.DummyFieldStepInterpolator;
31 import org.hipparchus.ode.sampling.FieldODEStateInterpolator;
32 import org.hipparchus.util.Binary64Field;
33 import org.hipparchus.util.FastMath;
34 import org.hipparchus.util.MathArrays;
35 import org.hipparchus.util.MathUtils;
36 import org.junit.jupiter.api.Test;
37
38 import java.util.Random;
39
40 import static org.junit.jupiter.api.Assertions.assertEquals;
41 import static org.junit.jupiter.api.Assertions.assertFalse;
42 import static org.junit.jupiter.api.Assertions.assertTrue;
43
44 public class FieldDenseOutputModelTest {
45
46 @Test
47 void testBoundaries() {
48 doTestBoundaries(Binary64Field.getInstance());
49 }
50
51 private <T extends CalculusFieldElement<T>> void doTestBoundaries(final Field<T> field) {
52 TestFieldProblem3<T> pb = new TestFieldProblem3<T>(field.getZero().add(0.9));
53 double minStep = 0;
54 double maxStep = pb.getFinalTime().subtract(pb.getInitialState().getTime()).getReal();
55 FieldODEIntegrator<T> integ = new DormandPrince54FieldIntegrator<T>(field, minStep, maxStep, 1.0e-8, 1.0e-8);
56 integ.addStepHandler(new FieldDenseOutputModel<T>());
57 integ.integrate(new FieldExpandableODE<T>(pb), pb.getInitialState(), pb.getFinalTime());
58 FieldDenseOutputModel<T> cm = (FieldDenseOutputModel<T>) integ.getStepHandlers().iterator().next();
59 cm.getInterpolatedState(pb.getInitialState().getTime().multiply(2).subtract(pb.getFinalTime()));
60 cm.getInterpolatedState(pb.getFinalTime().multiply(2).subtract(pb.getInitialState().getTime()));
61 cm.getInterpolatedState(pb.getInitialState().getTime().add(pb.getFinalTime()).multiply(0.5));
62 }
63
64 @Test
65 void testRandomAccess() {
66 doTestRandomAccess(Binary64Field.getInstance());
67 }
68
69 private <T extends CalculusFieldElement<T>> void doTestRandomAccess(final Field<T> field) {
70
71 TestFieldProblem3<T> pb = new TestFieldProblem3<T>(field.getZero().add(0.9));
72 double minStep = 0;
73 double maxStep = pb.getFinalTime().subtract(pb.getInitialState().getTime()).getReal();
74 FieldODEIntegrator<T> integ = new DormandPrince54FieldIntegrator<T>(field, minStep, maxStep, 1.0e-8, 1.0e-8);
75 FieldDenseOutputModel<T> cm = new FieldDenseOutputModel<T>();
76 integ.addStepHandler(cm);
77 integ.integrate(new FieldExpandableODE<T>(pb), pb.getInitialState(), pb.getFinalTime());
78
79 Random random = new Random(347588535632l);
80 T maxError = field.getZero();
81 T maxErrorDot = field.getZero();
82 for (int i = 0; i < 1000; ++i) {
83 double r = random.nextDouble();
84 T time = pb.getInitialState().getTime().multiply(r).add(pb.getFinalTime().multiply(1.0 - r));
85 FieldODEStateAndDerivative<T> interpolated = cm.getInterpolatedState(time);
86 T[] theoreticalY = pb.computeTheoreticalState(time);
87 T[] theoreticalYDot = pb.doComputeDerivatives(time, theoreticalY);
88 T dx = interpolated.getPrimaryState()[0].subtract(theoreticalY[0]);
89 T dy = interpolated.getPrimaryState()[1].subtract(theoreticalY[1]);
90 T error = dx.square().add(dy.square());
91 maxError = MathUtils.max(maxError, error);
92 T dxDot = interpolated.getPrimaryDerivative()[0].subtract(theoreticalYDot[0]);
93 T dyDot = interpolated.getPrimaryDerivative()[1].subtract(theoreticalYDot[1]);
94 T errorDot = dxDot.multiply(dxDot).add(dyDot.multiply(dyDot));
95 maxErrorDot = MathUtils.max(maxErrorDot, errorDot);
96 }
97
98 assertEquals(0.0, maxError.getReal(), 1.0e-9);
99 assertEquals(0.0, maxErrorDot.getReal(), 4.0e-7);
100
101 }
102
103 @Test
104 void testModelsMerging() {
105 doTestModelsMerging(Binary64Field.getInstance());
106 }
107
108 private <T extends CalculusFieldElement<T>> void doTestModelsMerging(final Field<T> field) {
109
110
111 FieldOrdinaryDifferentialEquation<T> problem =
112 new FieldOrdinaryDifferentialEquation<T>() {
113 public T[] computeDerivatives(T t, T[] y) {
114 T[] yDot = MathArrays.buildArray(field, 2);
115 yDot[0] = y[1].negate();
116 yDot[1] = y[0];
117 return yDot;
118 }
119 public int getDimension() {
120 return 2;
121 }
122 public void init(T t0, T[] y0, T finalTime) {
123 }
124 };
125
126
127 FieldDenseOutputModel<T> cm1 = new FieldDenseOutputModel<T>();
128 FieldODEIntegrator<T> integ1 =
129 new DormandPrince853FieldIntegrator<T>(field, 0, 1.0, 1.0e-8, 1.0e-8);
130 integ1.addStepHandler(cm1);
131 T t0 = field.getZero().add(FastMath.PI);
132 T[] y0 = MathArrays.buildArray(field, 2);
133 y0[0] = field.getOne().negate();
134 y0[1] = field.getZero();
135 integ1.integrate(new FieldExpandableODE<T>(problem),
136 new FieldODEState<T>(t0, y0),
137 field.getZero());
138
139
140 FieldDenseOutputModel<T> cm2 = new FieldDenseOutputModel<T>();
141 FieldODEIntegrator<T> integ2 =
142 new DormandPrince853FieldIntegrator<T>(field, 0, 0.1, 1.0e-12, 1.0e-12);
143 integ2.addStepHandler(cm2);
144 t0 = field.getZero().add(2.0 * FastMath.PI);
145 y0[0] = field.getOne();
146 y0[1] = field.getZero();
147 integ2.integrate(new FieldExpandableODE<T>(problem),
148 new FieldODEState<T>(t0, y0),
149 field.getZero().add(FastMath.PI));
150
151
152 FieldDenseOutputModel<T> cm = new FieldDenseOutputModel<T>();
153 cm.append(cm2);
154 cm.append(new FieldDenseOutputModel<T>());
155 cm.append(cm1);
156
157
158 assertEquals(2.0 * FastMath.PI, cm.getInitialTime().getReal(), 1.0e-12);
159 assertEquals(0, cm.getFinalTime().getReal(), 1.0e-12);
160 for (double t = 0; t < 2.0 * FastMath.PI; t += 0.1) {
161 FieldODEStateAndDerivative<T> interpolated = cm.getInterpolatedState(field.getZero().add(t));
162 assertEquals(FastMath.cos(t), interpolated.getPrimaryState()[0].getReal(), 1.0e-7);
163 assertEquals(FastMath.sin(t), interpolated.getPrimaryState()[1].getReal(), 1.0e-7);
164 }
165
166 }
167
168 @Test
169 void testErrorConditions() {
170 doTestErrorConditions(Binary64Field.getInstance());
171 }
172
173 private <T extends CalculusFieldElement<T>> void doTestErrorConditions(final Field<T> field) {
174 FieldDenseOutputModel<T> cm = new FieldDenseOutputModel<T>();
175 cm.handleStep(buildInterpolator(field, 0, 1, new double[] { 0.0, 1.0, -2.0 }));
176
177
178 assertTrue(checkAppendError(field, cm, 1.0, 2.0, new double[] { 0.0, 1.0 }));
179
180
181 assertTrue(checkAppendError(field, cm, 10.0, 20.0, new double[] { 0.0, 1.0, -2.0 }));
182
183
184 assertTrue(checkAppendError(field, cm, 1.0, 0.0, new double[] { 0.0, 1.0, -2.0 }));
185
186
187 assertFalse(checkAppendError(field, cm, 1.0, 2.0, new double[] { 0.0, 1.0, -2.0 }));
188
189 }
190
191 private <T extends CalculusFieldElement<T>> boolean checkAppendError(Field<T> field, FieldDenseOutputModel<T> cm,
192 double t0, double t1, double[] y) {
193 try {
194 FieldDenseOutputModel<T> otherCm = new FieldDenseOutputModel<T>();
195 otherCm.handleStep(buildInterpolator(field, t0, t1, y));
196 cm.append(otherCm);
197 } catch(MathIllegalArgumentException dme) {
198 return true;
199 }
200 return false;
201 }
202
203 private <T extends CalculusFieldElement<T>> FieldODEStateInterpolator<T> buildInterpolator(Field<T> field,
204 double t0, double t1, double[] y) {
205 T[] fieldY = MathArrays.buildArray(field, y.length);
206 for (int i = 0; i < y.length; ++i) {
207 fieldY[i] = field.getZero().add(y[i]);
208 }
209 final FieldODEStateAndDerivative<T> s0 = new FieldODEStateAndDerivative<T>(field.getZero().add(t0), fieldY, fieldY);
210 final FieldODEStateAndDerivative<T> s1 = new FieldODEStateAndDerivative<T>(field.getZero().add(t1), fieldY, fieldY);
211 final FieldEquationsMapper<T> mapper = new FieldExpandableODE<T>(new FieldOrdinaryDifferentialEquation<T>() {
212 public int getDimension() {
213 return s0.getPrimaryStateDimension();
214 }
215 public void init(T t0, T[] y0, T finalTime) {
216 }
217 public T[] computeDerivatives(T t, T[] y) {
218 return y;
219 }
220 }).getMapper();
221 return new DummyFieldStepInterpolator<T>(t1 >= t0, s0, s1, s0, s1, mapper);
222 }
223
224 public void checkValue(double value, double reference) {
225 assertTrue(FastMath.abs(value - reference) < 1.0e-10);
226 }
227
228 }