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