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.math.analysis.interpolation;
18  
19  import org.apache.commons.math.MathException;
20  import org.apache.commons.math.analysis.Expm1Function;
21  import org.apache.commons.math.analysis.SinFunction;
22  import org.apache.commons.math.analysis.UnivariateRealFunction;
23  
24  import junit.framework.TestCase;
25  
26  /**
27   * Testcase for Divided Difference interpolator.
28   * <p>
29   * The error of polynomial interpolation is
30   *     f(z) - p(z) = f^(n)(zeta) * (z-x[0])(z-x[1])...(z-x[n-1]) / n!
31   * where f^(n) is the n-th derivative of the approximated function and
32   * zeta is some point in the interval determined by x[] and z.
33   * <p>
34   * Since zeta is unknown, f^(n)(zeta) cannot be calculated. But we can bound
35   * it and use the absolute value upper bound for estimates. For reference,
36   * see <b>Introduction to Numerical Analysis</b>, ISBN 038795452X, chapter 2.
37   * 
38   * @version $Revision: 799857 $ $Date: 2009-08-01 09:07:12 -0400 (Sat, 01 Aug 2009) $ 
39   */
40  public final class DividedDifferenceInterpolatorTest extends TestCase {
41  
42      /**
43       * Test of interpolator for the sine function.
44       * <p>
45       * |sin^(n)(zeta)| <= 1.0, zeta in [0, 2*PI]
46       */
47      public void testSinFunction() throws MathException {
48          UnivariateRealFunction f = new SinFunction();
49          UnivariateRealInterpolator interpolator = new DividedDifferenceInterpolator();
50          double x[], y[], z, expected, result, tolerance;
51  
52          // 6 interpolating points on interval [0, 2*PI]
53          int n = 6;
54          double min = 0.0, max = 2 * Math.PI;
55          x = new double[n];
56          y = new double[n];
57          for (int i = 0; i < n; i++) {
58              x[i] = min + i * (max - min) / n;
59              y[i] = f.value(x[i]);
60          }
61          double derivativebound = 1.0;
62          UnivariateRealFunction p = interpolator.interpolate(x, y);
63  
64          z = Math.PI / 4; expected = f.value(z); result = p.value(z);
65          tolerance = Math.abs(derivativebound * partialerror(x, z));
66          assertEquals(expected, result, tolerance);
67  
68          z = Math.PI * 1.5; expected = f.value(z); result = p.value(z);
69          tolerance = Math.abs(derivativebound * partialerror(x, z));
70          assertEquals(expected, result, tolerance);
71      }
72  
73      /**
74       * Test of interpolator for the exponential function.
75       * <p>
76       * |expm1^(n)(zeta)| <= e, zeta in [-1, 1]
77       */
78      public void testExpm1Function() throws MathException {
79          UnivariateRealFunction f = new Expm1Function();
80          UnivariateRealInterpolator interpolator = new DividedDifferenceInterpolator();
81          double x[], y[], z, expected, result, tolerance;
82  
83          // 5 interpolating points on interval [-1, 1]
84          int n = 5;
85          double min = -1.0, max = 1.0;
86          x = new double[n];
87          y = new double[n];
88          for (int i = 0; i < n; i++) {
89              x[i] = min + i * (max - min) / n;
90              y[i] = f.value(x[i]);
91          }
92          double derivativebound = Math.E;
93          UnivariateRealFunction p = interpolator.interpolate(x, y);
94  
95          z = 0.0; expected = f.value(z); result = p.value(z);
96          tolerance = Math.abs(derivativebound * partialerror(x, z));
97          assertEquals(expected, result, tolerance);
98  
99          z = 0.5; expected = f.value(z); result = p.value(z);
100         tolerance = Math.abs(derivativebound * partialerror(x, z));
101         assertEquals(expected, result, tolerance);
102 
103         z = -0.5; expected = f.value(z); result = p.value(z);
104         tolerance = Math.abs(derivativebound * partialerror(x, z));
105         assertEquals(expected, result, tolerance);
106     }
107 
108     /**
109      * Test of parameters for the interpolator.
110      */
111     public void testParameters() throws Exception {
112         UnivariateRealInterpolator interpolator = new DividedDifferenceInterpolator();
113 
114         try {
115             // bad abscissas array
116             double x[] = { 1.0, 2.0, 2.0, 4.0 };
117             double y[] = { 0.0, 4.0, 4.0, 2.5 };
118             UnivariateRealFunction p = interpolator.interpolate(x, y);
119             p.value(0.0);
120             fail("Expecting MathException - bad abscissas array");
121         } catch (MathException ex) {
122             // expected
123         }
124     }
125 
126     /**
127      * Returns the partial error term (z-x[0])(z-x[1])...(z-x[n-1])/n!
128      */
129     protected double partialerror(double x[], double z) throws
130         IllegalArgumentException {
131 
132         if (x.length < 1) {
133             throw new IllegalArgumentException
134                 ("Interpolation array cannot be empty.");
135         }
136         double out = 1;
137         for (int i = 0; i < x.length; i++) {
138             out *= (z - x[i]) / (i + 1);
139         }
140         return out;
141     }
142 }