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