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