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.interpolators;
19  
20  
21  import org.hipparchus.analysis.differentiation.DerivativeStructure;
22  import org.hipparchus.ode.ODEStateAndDerivative;
23  import org.hipparchus.ode.OrdinaryDifferentialEquation;
24  import org.hipparchus.ode.sampling.AbstractODEStateInterpolator;
25  import org.hipparchus.ode.sampling.ODEStateInterpolator;
26  import org.hipparchus.util.FastMath;
27  import org.junit.jupiter.api.Test;
28  
29  import static org.junit.jupiter.api.Assertions.assertEquals;
30  import static org.junit.jupiter.api.Assertions.assertFalse;
31  import static org.junit.jupiter.api.Assertions.assertNotSame;
32  import static org.junit.jupiter.api.Assertions.assertSame;
33  import static org.junit.jupiter.api.Assertions.assertTrue;
34  
35  public abstract class ODEStateInterpolatorAbstractTest {
36  
37      @Test
38      public abstract void interpolationAtBounds();
39  
40      protected void doInterpolationAtBounds(double epsilon) {
41          ODEStateInterpolator interpolator = setUpInterpolator(new SinCos(), 0.0, new double[] { 0.0, 1.0 }, 0.125);
42  
43          assertEquals(0.0, interpolator.getPreviousState().getTime(), 1.0e-15);
44          for (int i = 0; i < 2; ++i) {
45              assertEquals(interpolator.getPreviousState().getPrimaryState()[i],
46                                  interpolator.getInterpolatedState(interpolator.getPreviousState().getTime()).getPrimaryState()[i],
47                                  epsilon);
48          }
49          assertEquals(0.125, interpolator.getCurrentState().getTime(), 1.0e-15);
50          for (int i = 0; i < 2; ++i) {
51              assertEquals(interpolator.getCurrentState().getPrimaryState()[i],
52                                  interpolator.getInterpolatedState(interpolator.getCurrentState().getTime()).getPrimaryState()[i],
53                                  epsilon);
54          }
55          assertFalse(interpolator.isPreviousStateInterpolated());
56          assertFalse(interpolator.isCurrentStateInterpolated());
57      }
58  
59      @Test
60      public abstract void interpolationInside();
61  
62      protected void doInterpolationInside(double epsilonSin, double epsilonCos) {
63  
64          SinCos sinCos = new SinCos();
65          ODEStateInterpolator interpolator = setUpInterpolator(sinCos, 0.0, new double[] { 0.0, 1.0 }, 0.0125);
66  
67          int n = 100;
68          double maxErrorSin = 0;
69          double maxErrorCos = 0;
70          for (int i = 0; i <= n; ++i) {
71              double t =     ((n - i) * interpolator.getPreviousState().getTime() +
72                                   i  * interpolator.getCurrentState().getTime()) / n;
73              final double[] interpolated = interpolator.getInterpolatedState(t).getPrimaryState();
74              double[] reference = sinCos.theoreticalState(t);
75              maxErrorSin = FastMath.max(maxErrorSin, FastMath.abs(interpolated[0] - reference[0]));
76              maxErrorCos = FastMath.max(maxErrorCos, FastMath.abs(interpolated[1] - reference[1]));
77          }
78  
79          assertEquals(0.0, maxErrorSin, epsilonSin);
80          assertEquals(0.0, maxErrorCos, epsilonCos);
81  
82          assertFalse(interpolator.isPreviousStateInterpolated());
83          assertFalse(interpolator.isCurrentStateInterpolated());
84      }
85  
86      @Test
87      public abstract void restrictPrevious();
88  
89      protected void doRestrictPrevious(double epsilon, double epsilonDot) {
90  
91          AbstractODEStateInterpolator original   = setUpInterpolator(new SinCos(), 0.0, new double[] { 0.0, 1.0 }, 0.125);
92  
93          assertFalse(original.isPreviousStateInterpolated());
94          assertFalse(original.isCurrentStateInterpolated());
95  
96          AbstractODEStateInterpolator restricted = original.restrictStep(original.getInterpolatedState(1.0 / 32),
97                                                                          original.getCurrentState());
98  
99          assertSame(original.getPreviousState(),       original.getGlobalPreviousState());
100         assertSame(original.getCurrentState(),        original.getGlobalCurrentState());
101         assertSame(original.getGlobalPreviousState(), restricted.getGlobalPreviousState());
102         assertSame(original.getGlobalCurrentState(),  restricted.getGlobalCurrentState());
103         assertNotSame(restricted.getPreviousState(),  restricted.getGlobalPreviousState());
104         assertSame(restricted.getCurrentState(),      restricted.getGlobalCurrentState());
105         assertEquals(1.0 / 32, restricted.getPreviousState().getTime(), 1.0e-15);
106         assertTrue(restricted.isPreviousStateInterpolated());
107         assertFalse(restricted.isCurrentStateInterpolated());
108 
109         checkRestricted(original, restricted, epsilon, epsilonDot);
110 
111     }
112 
113     @Test
114     public abstract void restrictCurrent();
115 
116     protected void doRestrictCurrent(double epsilon, double epsilonDot) {
117 
118         AbstractODEStateInterpolator original   = setUpInterpolator(new SinCos(), 0.0, new double[] { 0.0, 1.0 }, 0.125);
119 
120         assertFalse(original.isPreviousStateInterpolated());
121         assertFalse(original.isCurrentStateInterpolated());
122 
123         AbstractODEStateInterpolator restricted = original.restrictStep(original.getPreviousState(),
124                                                                         original.getInterpolatedState(3.0 / 32));
125 
126         assertSame(original.getPreviousState(),       original.getGlobalPreviousState());
127         assertSame(original.getCurrentState(),        original.getGlobalCurrentState());
128         assertSame(original.getGlobalPreviousState(), restricted.getGlobalPreviousState());
129         assertSame(original.getGlobalCurrentState(),  restricted.getGlobalCurrentState());
130         assertSame(restricted.getPreviousState(),     restricted.getGlobalPreviousState());
131         assertNotSame(restricted.getCurrentState(),   restricted.getGlobalCurrentState());
132         assertEquals(3.0 / 32, restricted.getCurrentState().getTime(), 1.0e-15);
133         assertFalse(restricted.isPreviousStateInterpolated());
134         assertTrue(restricted.isCurrentStateInterpolated());
135 
136         checkRestricted(original, restricted, epsilon, epsilonDot);
137 
138     }
139 
140     @Test
141     public abstract void restrictBothEnds();
142 
143     protected void doRestrictBothEnds(double epsilon, double epsilonDot) {
144 
145         AbstractODEStateInterpolator original   = setUpInterpolator(new SinCos(), 0.0, new double[] { 0.0, 1.0 }, 0.125);
146 
147         assertFalse(original.isPreviousStateInterpolated());
148         assertFalse(original.isCurrentStateInterpolated());
149 
150         AbstractODEStateInterpolator restricted = original.restrictStep(original.getInterpolatedState(1.0 / 32),
151                                                                         original.getInterpolatedState(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         assertNotSame(restricted.getPreviousState(),  restricted.getGlobalPreviousState());
158         assertNotSame(restricted.getCurrentState(),   restricted.getGlobalCurrentState());
159         assertEquals(1.0 / 32, restricted.getPreviousState().getTime(), 1.0e-15);
160         assertEquals(3.0 / 32, restricted.getCurrentState().getTime(), 1.0e-15);
161         assertTrue(restricted.isPreviousStateInterpolated());
162         assertTrue(restricted.isCurrentStateInterpolated());
163 
164         checkRestricted(original, restricted, epsilon, epsilonDot);
165 
166     }
167 
168     @Test
169     public abstract void degenerateInterpolation();
170 
171     protected void doDegenerateInterpolation() {
172         AbstractODEStateInterpolator interpolator = setUpInterpolator(new SinCos(), 0.0, new double[] { 0.0, 1.0 }, 0.0);
173         ODEStateAndDerivative interpolatedState = interpolator.getInterpolatedState(0.0);
174         assertEquals(0.0, interpolatedState.getTime(), 0.0);
175         assertEquals(0.0, interpolatedState.getPrimaryState()[0], 0.0);
176         assertEquals(1.0, interpolatedState.getPrimaryState()[1], 0.0);
177         assertEquals(1.0, interpolatedState.getPrimaryDerivative()[0], 0.0);
178         assertEquals(0.0, interpolatedState.getPrimaryDerivative()[1], 0.0);
179     }
180 
181     private void checkRestricted(AbstractODEStateInterpolator original, AbstractODEStateInterpolator restricted,
182                                  double epsilon, double epsilonDot) {
183         for (double t = restricted.getPreviousState().getTime();
184              t <= restricted.getCurrentState().getTime();
185              t += 1.0 / 256) {
186             ODEStateAndDerivative originalInterpolated   = original.getInterpolatedState(t);
187             ODEStateAndDerivative restrictedInterpolated = restricted.getInterpolatedState(t);
188             assertEquals(t, originalInterpolated.getTime(), 1.0e-15);
189             assertEquals(t, restrictedInterpolated.getTime(), 1.0e-15);
190             assertEquals(originalInterpolated.getPrimaryState()[0],
191                                 restrictedInterpolated.getPrimaryState()[0],
192                                 epsilon);
193             assertEquals(originalInterpolated.getPrimaryState()[1],
194                                 restrictedInterpolated.getPrimaryState()[1],
195                                 epsilon);
196             assertEquals(originalInterpolated.getPrimaryDerivative()[0],
197                                 restrictedInterpolated.getPrimaryDerivative()[0],
198                                 epsilonDot);
199             assertEquals(originalInterpolated.getPrimaryDerivative()[1],
200                                 restrictedInterpolated.getPrimaryDerivative()[1],
201                                 epsilonDot);
202         }
203 
204     }
205 
206     public interface ReferenceODE extends OrdinaryDifferentialEquation {
207         double[]              theoreticalState(double t);
208         DerivativeStructure[] theoreticalState(DerivativeStructure t);
209     }
210 
211     protected abstract AbstractODEStateInterpolator setUpInterpolator(final ReferenceODE eqn,
212                                                                       final double t0, final double[] y0,
213                                                                       final double t1);
214 
215     private static class SinCos implements ReferenceODE {
216         public int getDimension() {
217             return 2;
218         }
219         public double[] computeDerivatives(final double t, final double[] y) {
220             return new double[] {
221                 y[1], -y[0]
222             };
223         }
224         public double[] theoreticalState(final double t) {
225             return new double[] {
226                 FastMath.sin(t), FastMath.cos(t)
227             };
228         }
229         public DerivativeStructure[] theoreticalState(final DerivativeStructure t) {
230             return new DerivativeStructure[] {
231                 t.sin(), t.cos()
232             };
233         }
234     }
235 
236 }