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 org.apache.commons.math4.legacy.analysis.polynomials.PolynomialFunctionLagrangeForm;
20 import org.apache.commons.math4.legacy.analysis.polynomials.PolynomialFunctionNewtonForm;
21 import org.apache.commons.math4.legacy.exception.DimensionMismatchException;
22 import org.apache.commons.math4.legacy.exception.NonMonotonicSequenceException;
23 import org.apache.commons.math4.legacy.exception.NumberIsTooSmallException;
24
25 /**
26 * Implements the <a href=
27 * "http://mathworld.wolfram.com/NewtonsDividedDifferenceInterpolationFormula.html">
28 * Divided Difference Algorithm</a> for interpolation of real univariate
29 * functions. For reference, see <b>Introduction to Numerical Analysis</b>,
30 * ISBN 038795452X, chapter 2.
31 * <p>
32 * The actual code of Neville's evaluation is in PolynomialFunctionLagrangeForm,
33 * this class provides an easy-to-use interface to it.</p>
34 *
35 * @since 1.2
36 */
37 public class DividedDifferenceInterpolator
38 implements UnivariateInterpolator {
39 /**
40 * Compute an interpolating function for the dataset.
41 *
42 * @param x Interpolating points array.
43 * @param y Interpolating values array.
44 * @return a function which interpolates the dataset.
45 * @throws DimensionMismatchException if the array lengths are different.
46 * @throws NumberIsTooSmallException if the number of points is less than 2.
47 * @throws NonMonotonicSequenceException if {@code x} is not sorted in
48 * strictly increasing order.
49 */
50 @Override
51 public PolynomialFunctionNewtonForm interpolate(double[] x, double[] y)
52 throws DimensionMismatchException,
53 NumberIsTooSmallException,
54 NonMonotonicSequenceException {
55 /*
56 * a[] and c[] are defined in the general formula of Newton form:
57 * p(x) = a[0] + a[1](x-c[0]) + a[2](x-c[0])(x-c[1]) + ... +
58 * a[n](x-c[0])(x-c[1])...(x-c[n-1])
59 */
60 PolynomialFunctionLagrangeForm.verifyInterpolationArray(x, y, true);
61
62 /*
63 * When used for interpolation, the Newton form formula becomes
64 * p(x) = f[x0] + f[x0,x1](x-x0) + f[x0,x1,x2](x-x0)(x-x1) + ... +
65 * f[x0,x1,...,x[n-1]](x-x0)(x-x1)...(x-x[n-2])
66 * Therefore, a[k] = f[x0,x1,...,xk], c[k] = x[k].
67 * <p>
68 * Note x[], y[], a[] have the same length but c[]'s size is one less.</p>
69 */
70 final double[] c = new double[x.length-1];
71 System.arraycopy(x, 0, c, 0, c.length);
72
73 final double[] a = computeDividedDifference(x, y);
74 return new PolynomialFunctionNewtonForm(a, c);
75 }
76
77 /**
78 * Return a copy of the divided difference array.
79 * <p>
80 * The divided difference array is defined recursively by <pre>
81 * f[x0] = f(x0)
82 * f[x0,x1,...,xk] = (f[x1,...,xk] - f[x0,...,x[k-1]]) / (xk - x0)
83 * </pre>
84 * <p>
85 * The computational complexity is \(O(n^2)\) where \(n\) is the common
86 * length of {@code x} and {@code y}.</p>
87 *
88 * @param x Interpolating points array.
89 * @param y Interpolating values array.
90 * @return a fresh copy of the divided difference array.
91 * @throws DimensionMismatchException if the array lengths are different.
92 * @throws NumberIsTooSmallException if the number of points is less than 2.
93 * @throws NonMonotonicSequenceException
94 * if {@code x} is not sorted in strictly increasing order.
95 */
96 protected static double[] computeDividedDifference(final double[] x, final double[] y)
97 throws DimensionMismatchException,
98 NumberIsTooSmallException,
99 NonMonotonicSequenceException {
100 PolynomialFunctionLagrangeForm.verifyInterpolationArray(x, y, true);
101
102 final double[] divdiff = y.clone(); // initialization
103
104 final int n = x.length;
105 final double[] a = new double [n];
106 a[0] = divdiff[0];
107 for (int i = 1; i < n; i++) {
108 for (int j = 0; j < n-i; j++) {
109 final double denominator = x[j+i] - x[j];
110 divdiff[j] = (divdiff[j+1] - divdiff[j]) / denominator;
111 }
112 a[i] = divdiff[0];
113 }
114
115 return a;
116 }
117 }