View Javadoc
1   /*
2    * Licensed to the Apache Software Foundation (ASF) 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 ASF 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    *      http://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  package org.apache.commons.math4.legacy.analysis.interpolation;
18  
19  import java.util.Random;
20  
21  import org.apache.commons.math4.legacy.analysis.differentiation.DerivativeStructure;
22  import org.apache.commons.math4.legacy.analysis.polynomials.PolynomialFunction;
23  import org.apache.commons.math4.legacy.exception.MathIllegalArgumentException;
24  import org.apache.commons.math4.legacy.exception.NoDataException;
25  import org.apache.commons.math4.core.jdkmath.JdkMath;
26  import org.junit.Assert;
27  import org.junit.Test;
28  
29  public class HermiteInterpolatorTest {
30  
31      @Test
32      public void testZero() {
33          HermiteInterpolator interpolator = new HermiteInterpolator();
34          interpolator.addSamplePoint(0.0, new double[] { 0.0 });
35          for (double x = -10; x < 10; x += 1.0) {
36              DerivativeStructure y = interpolator.value(new DerivativeStructure(1, 1, 0, x))[0];
37              Assert.assertEquals(0.0, y.getValue(), 1.0e-15);
38              Assert.assertEquals(0.0, y.getPartialDerivative(1), 1.0e-15);
39          }
40          checkPolynomial(new PolynomialFunction(new double[] { 0.0 }),
41                          interpolator.getPolynomials()[0]);
42      }
43  
44      @Test
45      public void testQuadratic() {
46          HermiteInterpolator interpolator = new HermiteInterpolator();
47          interpolator.addSamplePoint(0.0, new double[] { 2.0 });
48          interpolator.addSamplePoint(1.0, new double[] { 0.0 });
49          interpolator.addSamplePoint(2.0, new double[] { 0.0 });
50          for (double x = -10; x < 10; x += 1.0) {
51              DerivativeStructure y = interpolator.value(new DerivativeStructure(1, 1, 0, x))[0];
52              Assert.assertEquals((x - 1.0) * (x - 2.0), y.getValue(), 1.0e-15);
53              Assert.assertEquals(2 * x - 3.0, y.getPartialDerivative(1), 1.0e-15);
54          }
55          checkPolynomial(new PolynomialFunction(new double[] { 2.0, -3.0, 1.0 }),
56                          interpolator.getPolynomials()[0]);
57      }
58  
59      @Test
60      public void testMixedDerivatives() {
61          HermiteInterpolator interpolator = new HermiteInterpolator();
62          interpolator.addSamplePoint(0.0, new double[] { 1.0 }, new double[] { 2.0 });
63          interpolator.addSamplePoint(1.0, new double[] { 4.0 });
64          interpolator.addSamplePoint(2.0, new double[] { 5.0 }, new double[] { 2.0 });
65          Assert.assertEquals(4, interpolator.getPolynomials()[0].degree());
66          DerivativeStructure y0 = interpolator.value(new DerivativeStructure(1, 1, 0, 0.0))[0];
67          Assert.assertEquals(1.0, y0.getValue(), 1.0e-15);
68          Assert.assertEquals(2.0, y0.getPartialDerivative(1), 1.0e-15);
69          Assert.assertEquals(4.0, interpolator.value(1.0)[0], 1.0e-15);
70          DerivativeStructure y2 = interpolator.value(new DerivativeStructure(1, 1, 0, 2.0))[0];
71          Assert.assertEquals(5.0, y2.getValue(), 1.0e-15);
72          Assert.assertEquals(2.0, y2.getPartialDerivative(1), 1.0e-15);
73          checkPolynomial(new PolynomialFunction(new double[] { 1.0, 2.0, 4.0, -4.0, 1.0 }),
74                          interpolator.getPolynomials()[0]);
75      }
76  
77      @Test
78      public void testRandomPolynomialsValuesOnly() {
79  
80          Random random = new Random(0x42b1e7dbd361a932L);
81  
82          for (int i = 0; i < 100; ++i) {
83  
84              int maxDegree = 0;
85              PolynomialFunction[] p = new PolynomialFunction[5];
86              for (int k = 0; k < p.length; ++k) {
87                  int degree = random.nextInt(7);
88                  p[k] = randomPolynomial(degree, random);
89                  maxDegree = JdkMath.max(maxDegree, degree);
90              }
91  
92              HermiteInterpolator interpolator = new HermiteInterpolator();
93              for (int j = 0; j < 1 + maxDegree; ++j) {
94                  double x = 0.1 * j;
95                  double[] values = new double[p.length];
96                  for (int k = 0; k < p.length; ++k) {
97                      values[k] = p[k].value(x);
98                  }
99                  interpolator.addSamplePoint(x, values);
100             }
101 
102             for (double x = 0; x < 2; x += 0.1) {
103                 double[] values = interpolator.value(x);
104                 Assert.assertEquals(p.length, values.length);
105                 for (int k = 0; k < p.length; ++k) {
106                     Assert.assertEquals(p[k].value(x), values[k], 1.0e-8 * JdkMath.abs(p[k].value(x)));
107                 }
108             }
109 
110             PolynomialFunction[] result = interpolator.getPolynomials();
111             for (int k = 0; k < p.length; ++k) {
112                 checkPolynomial(p[k], result[k]);
113             }
114         }
115     }
116 
117     @Test
118     public void testRandomPolynomialsFirstDerivative() {
119 
120         Random random = new Random(0x570803c982ca5d3bL);
121 
122         for (int i = 0; i < 100; ++i) {
123 
124             int maxDegree = 0;
125             PolynomialFunction[] p      = new PolynomialFunction[5];
126             PolynomialFunction[] pPrime = new PolynomialFunction[5];
127             for (int k = 0; k < p.length; ++k) {
128                 int degree = random.nextInt(7);
129                 p[k]      = randomPolynomial(degree, random);
130                 pPrime[k] = p[k].polynomialDerivative();
131                 maxDegree = JdkMath.max(maxDegree, degree);
132             }
133 
134             HermiteInterpolator interpolator = new HermiteInterpolator();
135             for (int j = 0; j < 1 + maxDegree / 2; ++j) {
136                 double x = 0.1 * j;
137                 double[] values      = new double[p.length];
138                 double[] derivatives = new double[p.length];
139                 for (int k = 0; k < p.length; ++k) {
140                     values[k]      = p[k].value(x);
141                     derivatives[k] = pPrime[k].value(x);
142                 }
143                 interpolator.addSamplePoint(x, values, derivatives);
144             }
145 
146             for (double x = 0; x < 2; x += 0.1) {
147                 DerivativeStructure[] y = interpolator.value(new DerivativeStructure(1, 1, 0, x));
148                 Assert.assertEquals(p.length, y.length);
149                 for (int k = 0; k < p.length; ++k) {
150                     Assert.assertEquals(p[k].value(x), y[k].getValue(), 1.0e-8 * JdkMath.abs(p[k].value(x)));
151                     Assert.assertEquals(pPrime[k].value(x), y[k].getPartialDerivative(1), 4.0e-8 * JdkMath.abs(p[k].value(x)));
152                 }
153             }
154 
155             PolynomialFunction[] result = interpolator.getPolynomials();
156             for (int k = 0; k < p.length; ++k) {
157                 checkPolynomial(p[k], result[k]);
158             }
159         }
160     }
161 
162     @Test
163     public void testSine() {
164         HermiteInterpolator interpolator = new HermiteInterpolator();
165         for (double x = 0; x < JdkMath.PI; x += 0.5) {
166             interpolator.addSamplePoint(x, new double[] { JdkMath.sin(x) });
167         }
168         for (double x = 0.1; x <= 2.9; x += 0.01) {
169             DerivativeStructure y = interpolator.value(new DerivativeStructure(1, 2, 0, x))[0];
170             Assert.assertEquals( JdkMath.sin(x), y.getValue(), 3.5e-5);
171             Assert.assertEquals( JdkMath.cos(x), y.getPartialDerivative(1), 1.3e-4);
172             Assert.assertEquals(-JdkMath.sin(x), y.getPartialDerivative(2), 2.9e-3);
173         }
174     }
175 
176     @Test
177     public void testSquareRoot() {
178         HermiteInterpolator interpolator = new HermiteInterpolator();
179         for (double x = 1.0; x < 3.6; x += 0.5) {
180             interpolator.addSamplePoint(x, new double[] { JdkMath.sqrt(x) });
181         }
182         for (double x = 1.1; x < 3.5; x += 0.01) {
183             DerivativeStructure y = interpolator.value(new DerivativeStructure(1, 1, 0, x))[0];
184             Assert.assertEquals(JdkMath.sqrt(x), y.getValue(), 1.5e-4);
185             Assert.assertEquals(0.5 / JdkMath.sqrt(x), y.getPartialDerivative(1), 8.5e-4);
186         }
187     }
188 
189     @Test
190     public void testWikipedia() {
191         // this test corresponds to the example from Wikipedia page:
192         // http://en.wikipedia.org/wiki/Hermite_interpolation
193         HermiteInterpolator interpolator = new HermiteInterpolator();
194         interpolator.addSamplePoint(-1, new double[] { 2 }, new double[] { -8 }, new double[] { 56 });
195         interpolator.addSamplePoint( 0, new double[] { 1 }, new double[] {  0 }, new double[] {  0 });
196         interpolator.addSamplePoint( 1, new double[] { 2 }, new double[] {  8 }, new double[] { 56 });
197         for (double x = -1.0; x <= 1.0; x += 0.125) {
198             DerivativeStructure y = interpolator.value(new DerivativeStructure(1, 1, 0, x))[0];
199             double x2 = x * x;
200             double x4 = x2 * x2;
201             double x8 = x4 * x4;
202             Assert.assertEquals(x8 + 1, y.getValue(), 1.0e-15);
203             Assert.assertEquals(8 * x4 * x2 * x, y.getPartialDerivative(1), 1.0e-15);
204         }
205         checkPolynomial(new PolynomialFunction(new double[] { 1, 0, 0, 0, 0, 0, 0, 0, 1 }),
206                         interpolator.getPolynomials()[0]);
207     }
208 
209     @Test
210     public void testOnePointParabola() {
211         HermiteInterpolator interpolator = new HermiteInterpolator();
212         interpolator.addSamplePoint(0, new double[] { 1 }, new double[] { 1 }, new double[] { 2 });
213         for (double x = -1.0; x <= 1.0; x += 0.125) {
214             DerivativeStructure y = interpolator.value(new DerivativeStructure(1, 1, 0, x))[0];
215             Assert.assertEquals(1 + x * (1 + x), y.getValue(), 1.0e-15);
216             Assert.assertEquals(1 + 2 * x, y.getPartialDerivative(1), 1.0e-15);
217         }
218         checkPolynomial(new PolynomialFunction(new double[] { 1, 1, 1 }),
219                         interpolator.getPolynomials()[0]);
220     }
221 
222     private PolynomialFunction randomPolynomial(int degree, Random random) {
223         double[] coeff = new double[ 1 + degree];
224         for (int j = 0; j < degree; ++j) {
225             coeff[j] = random.nextDouble();
226         }
227         return new PolynomialFunction(coeff);
228     }
229 
230     @Test(expected=NoDataException.class)
231     public void testEmptySample() {
232         new HermiteInterpolator().value(0.0);
233     }
234 
235     @Test(expected=MathIllegalArgumentException.class)
236     public void testDuplicatedAbscissa() {
237         HermiteInterpolator interpolator = new HermiteInterpolator();
238         interpolator.addSamplePoint(1.0, new double[] { 0.0 });
239         interpolator.addSamplePoint(1.0, new double[] { 1.0 });
240     }
241 
242     private void checkPolynomial(PolynomialFunction expected, PolynomialFunction result) {
243         Assert.assertTrue(result.degree() >= expected.degree());
244         double[] cE = expected.getCoefficients();
245         double[] cR = result.getCoefficients();
246         for (int i = 0; i < cE.length; ++i) {
247             Assert.assertEquals(cE[i], cR[i], 1.0e-8 * JdkMath.abs(cE[i]));
248         }
249         for (int i = cE.length; i < cR.length; ++i) {
250             Assert.assertEquals(0.0, cR[i], 1.0e-9);
251         }
252     }
253 }
254