View Javadoc
1   /*
2    * Licensed to the Hipparchus project under one or more
3    * contributor license agreements.  See the NOTICE file distributed with
4    * this work for additional information regarding copyright ownership.
5    * The Hipparchus project licenses this file to You under the Apache License, Version 2.0
6    * (the "License"); you may not use this file except in compliance with
7    * the License.  You may obtain a copy of the License at
8    *
9    *      https://www.apache.org/licenses/LICENSE-2.0
10   *
11   * Unless required by applicable law or agreed to in writing, software
12   * distributed under the License is distributed on an "AS IS" BASIS,
13   * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
14   * See the License for the specific language governing permissions and
15   * limitations under the License.
16   */
17  
18  package org.hipparchus.ode.nonstiff;
19  
20  import org.hipparchus.CalculusFieldElement;
21  import org.hipparchus.Field;
22  import org.hipparchus.ode.EquationsMapper;
23  import org.hipparchus.ode.FieldEquationsMapper;
24  import org.hipparchus.ode.FieldODEStateAndDerivative;
25  import org.hipparchus.ode.FieldOrdinaryDifferentialEquation;
26  import org.hipparchus.ode.ODEStateAndDerivative;
27  import org.hipparchus.ode.sampling.AbstractFieldODEStateInterpolator;
28  import org.hipparchus.ode.sampling.FieldODEStateInterpolator;
29  import org.hipparchus.ode.sampling.ODEStateInterpolator;
30  import org.hipparchus.util.Binary64Field;
31  import org.hipparchus.util.FastMath;
32  import org.hipparchus.util.MathArrays;
33  import org.junit.jupiter.api.Test;
34  
35  import java.lang.reflect.InvocationTargetException;
36  
37  import static org.junit.jupiter.api.Assertions.assertEquals;
38  import static org.junit.jupiter.api.Assertions.assertFalse;
39  import static org.junit.jupiter.api.Assertions.assertNotSame;
40  import static org.junit.jupiter.api.Assertions.assertSame;
41  import static org.junit.jupiter.api.Assertions.assertTrue;
42  
43  public abstract class FieldODEStateInterpolatorAbstractTest {
44  
45      @Test
46      public abstract void interpolationAtBounds();
47  
48      protected <T extends CalculusFieldElement<T>> void doInterpolationAtBounds(final Field<T> field, double epsilon) {
49  
50          FieldODEStateInterpolator<T> interpolator = setUpInterpolator(field,
51                                                                        new SinCos<T>(field),
52                                                                        0.0, new double[] { 0.0, 1.0 }, 0.125);
53  
54          assertEquals(0.0, interpolator.getPreviousState().getTime().getReal(), 1.0e-15);
55          for (int i = 0; i < 2; ++i) {
56              assertEquals(interpolator.getPreviousState().getPrimaryState()[i].getReal(),
57                                  interpolator.getInterpolatedState(interpolator.getPreviousState().getTime()).getPrimaryState()[i].getReal(),
58                                  epsilon);
59          }
60          assertEquals(0.125, interpolator.getCurrentState().getTime().getReal(), 1.0e-15);
61          for (int i = 0; i < 2; ++i) {
62              assertEquals(interpolator.getCurrentState().getPrimaryState()[i].getReal(),
63                                  interpolator.getInterpolatedState(interpolator.getCurrentState().getTime()).getPrimaryState()[i].getReal(),
64                                  epsilon);
65          }
66  
67          assertFalse(interpolator.isPreviousStateInterpolated());
68          assertFalse(interpolator.isCurrentStateInterpolated());
69      }
70  
71      @Test
72      public abstract void interpolationInside();
73  
74      protected <T extends CalculusFieldElement<T>> void doInterpolationInside(final Field<T> field,
75                                                                           double epsilonSin, double epsilonCos) {
76  
77          ReferenceFieldODE<T> sinCos =  new SinCos<T>(field);
78          FieldODEStateInterpolator<T> interpolator = setUpInterpolator(field, sinCos,
79                                                                        0.0, new double[] { 0.0, 1.0 }, 0.0125);
80  
81          int n = 100;
82          double maxErrorSin = 0;
83          double maxErrorCos = 0;
84          for (int i = 0; i <= n; ++i) {
85              T t =     interpolator.getPreviousState().getTime().multiply(n - i).
86                    add(interpolator.getCurrentState().getTime().multiply(i)).
87                    divide(n);
88              FieldODEStateAndDerivative<T> state = interpolator.getInterpolatedState(t);
89              T[] ref = sinCos.theoreticalState(t);
90              maxErrorSin = FastMath.max(maxErrorSin, state.getPrimaryState()[0].subtract(ref[0]).norm());
91              maxErrorCos = FastMath.max(maxErrorCos, state.getPrimaryState()[1].subtract(ref[1]).norm());
92          }
93          assertEquals(0.0, maxErrorSin, epsilonSin);
94          assertEquals(0.0, maxErrorCos, epsilonCos);
95  
96          assertFalse(interpolator.isPreviousStateInterpolated());
97          assertFalse(interpolator.isCurrentStateInterpolated());
98      }
99  
100     @Test
101     public void restrictPrevious() {
102         doRestrictPrevious(Binary64Field.getInstance(), 1e-15, 1e-15);
103     }
104 
105     protected <T extends CalculusFieldElement<T>> void doRestrictPrevious(
106             Field<T> field,
107             double epsilon,
108             double epsilonDot) {
109 
110         AbstractFieldODEStateInterpolator<T> original = setUpInterpolator(
111                 field, new SinCos<>(field), 0.0, new double[]{0.0, 1.0}, 0.125);
112 
113         assertFalse(original.isPreviousStateInterpolated());
114         assertFalse(original.isCurrentStateInterpolated());
115 
116         AbstractFieldODEStateInterpolator<T> restricted = original.restrictStep(
117                 original.getInterpolatedState(field.getZero().add(1.0 / 32)),
118                 original.getCurrentState());
119 
120         assertSame(original.getPreviousState(),       original.getGlobalPreviousState());
121         assertSame(original.getCurrentState(),        original.getGlobalCurrentState());
122         assertSame(original.getGlobalPreviousState(), restricted.getGlobalPreviousState());
123         assertSame(original.getGlobalCurrentState(),  restricted.getGlobalCurrentState());
124         assertNotSame(restricted.getPreviousState(),  restricted.getGlobalPreviousState());
125         assertSame(restricted.getCurrentState(),      restricted.getGlobalCurrentState());
126         assertEquals(1.0 / 32, restricted.getPreviousState().getTime().getReal(), 1.0e-15);
127         assertTrue(restricted.isPreviousStateInterpolated());
128         assertFalse(restricted.isCurrentStateInterpolated());
129 
130         checkRestricted(original, restricted, epsilon, epsilonDot);
131 
132     }
133 
134     @Test
135     public void restrictCurrent() {
136         doRestrictCurrent(Binary64Field.getInstance(), 1e-15, 1e-15);
137     }
138 
139     protected <T extends CalculusFieldElement<T>> void doRestrictCurrent(Field<T> field,
140                                                                      double epsilon,
141                                                                      double epsilonDot) {
142 
143         AbstractFieldODEStateInterpolator<T> original = setUpInterpolator(
144                 field, new SinCos<>(field), 0.0, new double[]{0.0, 1.0}, 0.125);
145 
146         assertFalse(original.isPreviousStateInterpolated());
147         assertFalse(original.isCurrentStateInterpolated());
148 
149         AbstractFieldODEStateInterpolator<T> restricted = original.restrictStep(
150                 original.getPreviousState(),
151                 original.getInterpolatedState(field.getZero().add(3.0 / 32)));
152 
153         assertSame(original.getPreviousState(),       original.getGlobalPreviousState());
154         assertSame(original.getCurrentState(),        original.getGlobalCurrentState());
155         assertSame(original.getGlobalPreviousState(), restricted.getGlobalPreviousState());
156         assertSame(original.getGlobalCurrentState(),  restricted.getGlobalCurrentState());
157         assertSame(restricted.getPreviousState(),     restricted.getGlobalPreviousState());
158         assertNotSame(restricted.getCurrentState(),   restricted.getGlobalCurrentState());
159         assertEquals(3.0 / 32, restricted.getCurrentState().getTime().getReal(), 1.0e-15);
160         assertFalse(restricted.isPreviousStateInterpolated());
161         assertTrue(restricted.isCurrentStateInterpolated());
162 
163         checkRestricted(original, restricted, epsilon, epsilonDot);
164 
165     }
166 
167     @Test
168     public void restrictBothEnds() {
169         doRestrictBothEnds(Binary64Field.getInstance(), 1e-15, 1e-15);
170     }
171 
172     protected <T extends CalculusFieldElement<T>> void doRestrictBothEnds(Field<T> field,
173                                                                       double epsilon,
174                                                                       double epsilonDot) {
175 
176         AbstractFieldODEStateInterpolator<T> original = setUpInterpolator(
177                 field, new SinCos<>(field), 0.0, new double[]{0.0, 1.0}, 0.125);
178 
179         assertFalse(original.isPreviousStateInterpolated());
180         assertFalse(original.isCurrentStateInterpolated());
181 
182         AbstractFieldODEStateInterpolator<T> restricted = original.restrictStep(
183                 original.getInterpolatedState(field.getZero().add(1.0 / 32)),
184                 original.getInterpolatedState(field.getZero().add(3.0 / 32)));
185 
186         assertSame(original.getPreviousState(),       original.getGlobalPreviousState());
187         assertSame(original.getCurrentState(),        original.getGlobalCurrentState());
188         assertSame(original.getGlobalPreviousState(), restricted.getGlobalPreviousState());
189         assertSame(original.getGlobalCurrentState(),  restricted.getGlobalCurrentState());
190         assertNotSame(restricted.getPreviousState(),  restricted.getGlobalPreviousState());
191         assertNotSame(restricted.getCurrentState(),   restricted.getGlobalCurrentState());
192         assertEquals(1.0 / 32, restricted.getPreviousState().getTime().getReal(), 1.0e-15);
193         assertEquals(3.0 / 32, restricted.getCurrentState().getTime().getReal(), 1.0e-15);
194         assertTrue(restricted.isPreviousStateInterpolated());
195         assertTrue(restricted.isCurrentStateInterpolated());
196 
197         checkRestricted(original, restricted, epsilon, epsilonDot);
198 
199     }
200 
201     @Test
202     public void degenerateInterpolation() {
203         doDegenerateInterpolation(Binary64Field.getInstance());
204     }
205 
206     protected <T extends CalculusFieldElement<T>> void doDegenerateInterpolation(Field<T> field) {
207         AbstractFieldODEStateInterpolator<T> interpolator = setUpInterpolator(
208                 field, new SinCos<>(field), 0.0, new double[] { 0.0, 1.0 }, 0.0);
209         FieldODEStateAndDerivative<T> interpolatedState = interpolator.getInterpolatedState(field.getZero());
210         assertEquals(0.0, interpolatedState.getTime().getReal(), 0.0);
211         assertEquals(0.0, interpolatedState.getPrimaryState()[0].getReal(), 0.0);
212         assertEquals(1.0, interpolatedState.getPrimaryState()[1].getReal(), 0.0);
213         assertEquals(1.0, interpolatedState.getPrimaryDerivative()[0].getReal(), 0.0);
214         assertEquals(0.0, interpolatedState.getPrimaryDerivative()[1].getReal(), 0.0);
215     }
216 
217     private <T extends CalculusFieldElement<T>> void checkRestricted(
218             AbstractFieldODEStateInterpolator<T> original,
219             AbstractFieldODEStateInterpolator<T> restricted,
220             double epsilon,
221             double epsilonDot) {
222 
223         for (T t = restricted.getPreviousState().getTime();
224              t.getReal() <= restricted.getCurrentState().getTime().getReal();
225              t = t.add(1.0 / 256)) {
226             FieldODEStateAndDerivative<T> originalInterpolated   = original.getInterpolatedState(t);
227             FieldODEStateAndDerivative<T> restrictedInterpolated = restricted.getInterpolatedState(t);
228             assertEquals(t.getReal(), originalInterpolated.getTime().getReal(), 1.0e-15);
229             assertEquals(t.getReal(), restrictedInterpolated.getTime().getReal(), 1.0e-15);
230             assertEquals(originalInterpolated.getPrimaryState()[0].getReal(),
231                                 restrictedInterpolated.getPrimaryState()[0].getReal(),
232                                 epsilon);
233             assertEquals(originalInterpolated.getPrimaryState()[1].getReal(),
234                                 restrictedInterpolated.getPrimaryState()[1].getReal(),
235                                 epsilon);
236             assertEquals(originalInterpolated.getPrimaryDerivative()[0].getReal(),
237                                 restrictedInterpolated.getPrimaryDerivative()[0].getReal(),
238                                 epsilonDot);
239             assertEquals(originalInterpolated.getPrimaryDerivative()[1].getReal(),
240                                 restrictedInterpolated.getPrimaryDerivative()[1].getReal(),
241                                 epsilonDot);
242         }
243 
244     }
245 
246     @Test
247     public abstract void nonFieldInterpolatorConsistency();
248 
249     protected <T extends CalculusFieldElement<T>> void doNonFieldInterpolatorConsistency(final Field<T> field,
250                                                                                      double epsilonSin, double epsilonCos,
251                                                                                      double epsilonSinDot, double epsilonCosDot) {
252 
253         ReferenceFieldODE<T> eqn = new SinCos<T>(field);
254         FieldODEStateInterpolator<T> fieldInterpolator =
255                         setUpInterpolator(field, eqn, 0.0, new double[] { 0.0, 1.0 }, 0.125);
256         ODEStateInterpolator regularInterpolator = convertInterpolator(fieldInterpolator, eqn);
257 
258         int n = 100;
259         double maxErrorSin    = 0;
260         double maxErrorCos    = 0;
261         double maxErrorSinDot = 0;
262         double maxErrorCosDot = 0;
263         for (int i = 0; i <= n; ++i) {
264 
265             T t =     fieldInterpolator.getPreviousState().getTime().multiply(n - i).
266                   add(fieldInterpolator.getCurrentState().getTime().multiply(i)).
267                   divide(n);
268 
269             FieldODEStateAndDerivative<T> state = fieldInterpolator.getInterpolatedState(t);
270             T[] fieldY    = state.getPrimaryState();
271             T[] fieldYDot = state.getPrimaryDerivative();
272 
273             ODEStateAndDerivative regularState = regularInterpolator.getInterpolatedState(t.getReal());
274             double[] regularY     = regularState.getPrimaryState();
275             double[] regularYDot  = regularState.getPrimaryDerivative();
276 
277             maxErrorSin    = FastMath.max(maxErrorSin,    fieldY[0].subtract(regularY[0]).norm());
278             maxErrorCos    = FastMath.max(maxErrorCos,    fieldY[1].subtract(regularY[1]).norm());
279             maxErrorSinDot = FastMath.max(maxErrorSinDot, fieldYDot[0].subtract(regularYDot[0]).norm());
280             maxErrorCosDot = FastMath.max(maxErrorCosDot, fieldYDot[1].subtract(regularYDot[1]).norm());
281 
282         }
283         assertEquals(0.0, maxErrorSin,    epsilonSin);
284         assertEquals(0.0, maxErrorCos,    epsilonCos);
285         assertEquals(0.0, maxErrorSinDot, epsilonSinDot);
286         assertEquals(0.0, maxErrorCosDot, epsilonCosDot);
287 
288     }
289 
290     public interface ReferenceFieldODE<T extends CalculusFieldElement<T>> extends FieldOrdinaryDifferentialEquation<T> {
291         T[] theoreticalState(T t);
292     }
293 
294     protected abstract <T extends CalculusFieldElement<T>>
295     AbstractFieldODEStateInterpolator<T> setUpInterpolator(final Field<T> field,
296                                                            final ReferenceFieldODE<T> eqn,
297                                                            final double t0,
298                                                            final double[] y0,
299                                                            final double t1);
300 
301     protected abstract <T extends CalculusFieldElement<T>>
302     ODEStateInterpolator convertInterpolator(final FieldODEStateInterpolator<T> fieldInterpolator,
303                                              final FieldOrdinaryDifferentialEquation<T> eqn);
304 
305     protected <T extends CalculusFieldElement<T>>
306     ODEStateAndDerivative convertODEStateAndDerivative(final FieldODEStateAndDerivative<T> s) {
307         final double[][] secondaryStates;
308         final double[][] secondaryDerivatives;
309         if (s.getNumberOfSecondaryStates() == 0) {
310             secondaryStates      = null;
311             secondaryDerivatives = null;
312         } else {
313             secondaryStates      = new double[s.getNumberOfSecondaryStates()][];
314             secondaryDerivatives = new double[s.getNumberOfSecondaryStates()][];
315             for (int i = 0; i < secondaryStates.length; ++i) {
316                 secondaryStates[i]      = convertArray(s.getSecondaryState(i));
317                 secondaryDerivatives[i] = convertArray(s.getSecondaryDerivative(i));
318             }
319         }
320         return new ODEStateAndDerivative(s.getTime().getReal(),
321                                          convertArray(s.getPrimaryState()),
322                                          convertArray(s.getPrimaryDerivative()),
323                                          secondaryStates,
324                                          secondaryDerivatives);
325     }
326 
327     protected <T extends CalculusFieldElement<T>> double[][] convertArray(final T[][] fieldArray) {
328         if (fieldArray == null) {
329             return null;
330         }
331         double[][] array = new double[fieldArray.length][];
332         for (int i = 0; i < array.length; ++i) {
333             array[i] = convertArray(fieldArray[i]);
334         }
335         return array;
336     }
337 
338     protected <T extends CalculusFieldElement<T>> double[] convertArray(final T[] fieldArray) {
339         if (fieldArray == null) {
340             return null;
341         }
342         double[] array = new double[fieldArray.length];
343         for (int i = 0; i < array.length; ++i) {
344             array[i] = fieldArray[i].getReal();
345         }
346         return array;
347     }
348 
349     protected <T extends CalculusFieldElement<T>>
350     EquationsMapper convertMapper(final FieldEquationsMapper<T> fieldmapper)
351         throws NoSuchMethodException, SecurityException, NoSuchFieldException,
352                IllegalArgumentException, IllegalAccessException,
353                InstantiationException, InvocationTargetException {
354         java.lang.reflect.Field fStart = FieldEquationsMapper.class.getDeclaredField("start");
355         fStart.setAccessible(true);
356         int[] start = (int[]) fStart.get(fieldmapper);
357 
358         java.lang.reflect.Constructor<EquationsMapper> regularMapperConstructor =
359                         EquationsMapper.class.getDeclaredConstructor(EquationsMapper.class,
360                                                                      Integer.TYPE);
361         regularMapperConstructor.setAccessible(true);
362         EquationsMapper regularMapper = regularMapperConstructor.newInstance(null, start[1]);
363         for (int k = 1; k < fieldmapper.getNumberOfEquations(); ++k) {
364             regularMapper = regularMapperConstructor.newInstance(regularMapper,
365                                                                  start[k + 1] - start[k]);
366         }
367         return regularMapper;
368     }
369 
370     private static class SinCos<T extends CalculusFieldElement<T>> implements ReferenceFieldODE<T> {
371         private final Field<T> field;
372         protected SinCos(final Field<T> field) {
373             this.field = field;
374         }
375         public int getDimension() {
376             return 2;
377         }
378         public void init(final T t0, final T[] y0, final T finalTime) {
379         }
380         public T[] computeDerivatives(final T t, final T[] y) {
381             T[] yDot = MathArrays.buildArray(field, 2);
382             yDot[0] = y[1];
383             yDot[1] = y[0].negate();
384             return yDot;
385         }
386         public T[] theoreticalState(final T t) {
387             T[] state = MathArrays.buildArray(field, 2);
388             state[0] = t.sin();
389             state[1] = t.cos();
390             return state;
391         }
392     }
393 
394 }