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