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    *      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  /*
19   * This is not the original file distributed by the Apache Software Foundation
20   * It has been modified by the Hipparchus project
21   */
22  package org.hipparchus.analysis.interpolation;
23  
24  import org.hipparchus.analysis.UnivariateFunction;
25  import org.hipparchus.analysis.function.Expm1;
26  import org.hipparchus.analysis.function.Sin;
27  import org.hipparchus.exception.MathIllegalArgumentException;
28  import org.hipparchus.util.FastMath;
29  import org.junit.Assert;
30  import org.junit.Test;
31  
32  
33  /**
34   * Test case for Neville interpolator.
35   * <p>
36   * The error of polynomial interpolation is
37   *     f(z) - p(z) = f^(n)(zeta) * (z-x[0])(z-x[1])...(z-x[n-1]) / n!
38   * where f^(n) is the n-th derivative of the approximated function and
39   * zeta is some point in the interval determined by x[] and z.
40   * <p>
41   * Since zeta is unknown, f^(n)(zeta) cannot be calculated. But we can bound
42   * it and use the absolute value upper bound for estimates. For reference,
43   * see <b>Introduction to Numerical Analysis</b>, ISBN 038795452X, chapter 2.
44   *
45   */
46  public final class NevilleInterpolatorTest {
47  
48      /**
49       * Test of interpolator for the sine function.
50       * <p>
51       * |sin^(n)(zeta)| &lt;= 1.0, zeta in [0, 2*PI]
52       */
53      @Test
54      public void testSinFunction() {
55          UnivariateFunction f = new Sin();
56          UnivariateInterpolator interpolator = new NevilleInterpolator();
57          double[] x;
58          double[] y;
59          double z;
60          double expected;
61          double result;
62          double tolerance;
63  
64          // 6 interpolating points on interval [0, 2*PI]
65          int n = 6;
66          double min = 0.0, max = 2 * FastMath.PI;
67          x = new double[n];
68          y = new double[n];
69          for (int i = 0; i < n; i++) {
70              x[i] = min + i * (max - min) / n;
71              y[i] = f.value(x[i]);
72          }
73          double derivativebound = 1.0;
74          UnivariateFunction p = interpolator.interpolate(x, y);
75  
76          z = FastMath.PI / 4; expected = f.value(z); result = p.value(z);
77          tolerance = FastMath.abs(derivativebound * partialerror(x, z));
78          Assert.assertEquals(expected, result, tolerance);
79  
80          z = FastMath.PI * 1.5; expected = f.value(z); result = p.value(z);
81          tolerance = FastMath.abs(derivativebound * partialerror(x, z));
82          Assert.assertEquals(expected, result, tolerance);
83      }
84  
85      /**
86       * Test of interpolator for the exponential function.
87       * <p>
88       * |expm1^(n)(zeta)| &lt;= e, zeta in [-1, 1]
89       */
90      @Test
91      public void testExpm1Function() {
92          UnivariateFunction f = new Expm1();
93          UnivariateInterpolator interpolator = new NevilleInterpolator();
94          double[] x;
95          double[] y;
96          double z;
97          double expected;
98          double result;
99          double tolerance;
100 
101         // 5 interpolating points on interval [-1, 1]
102         int n = 5;
103         double min = -1.0, max = 1.0;
104         x = new double[n];
105         y = new double[n];
106         for (int i = 0; i < n; i++) {
107             x[i] = min + i * (max - min) / n;
108             y[i] = f.value(x[i]);
109         }
110         double derivativebound = FastMath.E;
111         UnivariateFunction p = interpolator.interpolate(x, y);
112 
113         z = 0.0; expected = f.value(z); result = p.value(z);
114         tolerance = FastMath.abs(derivativebound * partialerror(x, z));
115         Assert.assertEquals(expected, result, tolerance);
116 
117         z = 0.5; expected = f.value(z); result = p.value(z);
118         tolerance = FastMath.abs(derivativebound * partialerror(x, z));
119         Assert.assertEquals(expected, result, tolerance);
120 
121         z = -0.5; expected = f.value(z); result = p.value(z);
122         tolerance = FastMath.abs(derivativebound * partialerror(x, z));
123         Assert.assertEquals(expected, result, tolerance);
124     }
125 
126     /**
127      * Test of parameters for the interpolator.
128      */
129     @Test
130     public void testParameters() {
131         UnivariateInterpolator interpolator = new NevilleInterpolator();
132 
133         try {
134             // bad abscissas array
135             double[] x = { 1.0, 2.0, 2.0, 4.0 };
136             double[] y = { 0.0, 4.0, 4.0, 2.5 };
137             UnivariateFunction p = interpolator.interpolate(x, y);
138             p.value(0.0);
139             Assert.fail("Expecting MathIllegalArgumentException - bad abscissas array");
140         } catch (MathIllegalArgumentException ex) {
141             // expected
142         }
143     }
144 
145     /**
146      * Returns the partial error term (z-x[0])(z-x[1])...(z-x[n-1])/n!
147      */
148     protected double partialerror(double[] x, double z) throws
149         IllegalArgumentException {
150 
151         if (x.length < 1) {
152             throw new IllegalArgumentException
153                 ("Interpolation array cannot be empty.");
154         }
155         double out = 1;
156         for (int i = 0; i < x.length; i++) {
157             out *= (z - x[i]) / (i + 1);
158         }
159         return out;
160     }
161 }