001    /*
002     * Licensed to the Apache Software Foundation (ASF) under one or more
003     * contributor license agreements.  See the NOTICE file distributed with
004     * this work for additional information regarding copyright ownership.
005     * The ASF licenses this file to You under the Apache License, Version 2.0
006     * (the "License"); you may not use this file except in compliance with
007     * the License.  You may obtain a copy of the License at
008     *
009     *      http://www.apache.org/licenses/LICENSE-2.0
010     *
011     * Unless required by applicable law or agreed to in writing, software
012     * distributed under the License is distributed on an "AS IS" BASIS,
013     * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
014     * See the License for the specific language governing permissions and
015     * limitations under the License.
016     */
017    package org.apache.commons.math.analysis.interpolation;
018    
019    import org.apache.commons.math.MathException;
020    import org.apache.commons.math.analysis.Expm1Function;
021    import org.apache.commons.math.analysis.SinFunction;
022    import org.apache.commons.math.analysis.UnivariateRealFunction;
023    
024    import junit.framework.TestCase;
025    
026    /**
027     * Testcase for Divided Difference interpolator.
028     * <p>
029     * The error of polynomial interpolation is
030     *     f(z) - p(z) = f^(n)(zeta) * (z-x[0])(z-x[1])...(z-x[n-1]) / n!
031     * where f^(n) is the n-th derivative of the approximated function and
032     * zeta is some point in the interval determined by x[] and z.
033     * <p>
034     * Since zeta is unknown, f^(n)(zeta) cannot be calculated. But we can bound
035     * it and use the absolute value upper bound for estimates. For reference,
036     * see <b>Introduction to Numerical Analysis</b>, ISBN 038795452X, chapter 2.
037     * 
038     * @version $Revision: 799857 $ $Date: 2009-08-01 09:07:12 -0400 (Sat, 01 Aug 2009) $ 
039     */
040    public final class DividedDifferenceInterpolatorTest extends TestCase {
041    
042        /**
043         * Test of interpolator for the sine function.
044         * <p>
045         * |sin^(n)(zeta)| <= 1.0, zeta in [0, 2*PI]
046         */
047        public void testSinFunction() throws MathException {
048            UnivariateRealFunction f = new SinFunction();
049            UnivariateRealInterpolator interpolator = new DividedDifferenceInterpolator();
050            double x[], y[], z, expected, result, tolerance;
051    
052            // 6 interpolating points on interval [0, 2*PI]
053            int n = 6;
054            double min = 0.0, max = 2 * Math.PI;
055            x = new double[n];
056            y = new double[n];
057            for (int i = 0; i < n; i++) {
058                x[i] = min + i * (max - min) / n;
059                y[i] = f.value(x[i]);
060            }
061            double derivativebound = 1.0;
062            UnivariateRealFunction p = interpolator.interpolate(x, y);
063    
064            z = Math.PI / 4; expected = f.value(z); result = p.value(z);
065            tolerance = Math.abs(derivativebound * partialerror(x, z));
066            assertEquals(expected, result, tolerance);
067    
068            z = Math.PI * 1.5; expected = f.value(z); result = p.value(z);
069            tolerance = Math.abs(derivativebound * partialerror(x, z));
070            assertEquals(expected, result, tolerance);
071        }
072    
073        /**
074         * Test of interpolator for the exponential function.
075         * <p>
076         * |expm1^(n)(zeta)| <= e, zeta in [-1, 1]
077         */
078        public void testExpm1Function() throws MathException {
079            UnivariateRealFunction f = new Expm1Function();
080            UnivariateRealInterpolator interpolator = new DividedDifferenceInterpolator();
081            double x[], y[], z, expected, result, tolerance;
082    
083            // 5 interpolating points on interval [-1, 1]
084            int n = 5;
085            double min = -1.0, max = 1.0;
086            x = new double[n];
087            y = new double[n];
088            for (int i = 0; i < n; i++) {
089                x[i] = min + i * (max - min) / n;
090                y[i] = f.value(x[i]);
091            }
092            double derivativebound = Math.E;
093            UnivariateRealFunction p = interpolator.interpolate(x, y);
094    
095            z = 0.0; expected = f.value(z); result = p.value(z);
096            tolerance = Math.abs(derivativebound * partialerror(x, z));
097            assertEquals(expected, result, tolerance);
098    
099            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    }