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 }