Classes in this File | Line Coverage | Branch Coverage | Complexity | ||||||||
DividedDifferenceInterpolator |
|
| 4.0;4 |
1 | /* |
|
2 | * Copyright 2003-2005 The Apache Software Foundation. |
|
3 | * |
|
4 | * Licensed under the Apache License, Version 2.0 (the "License"); |
|
5 | * you may not use this file except in compliance with the License. |
|
6 | * You may obtain a copy of the License at |
|
7 | * |
|
8 | * http://www.apache.org/licenses/LICENSE-2.0 |
|
9 | * |
|
10 | * Unless required by applicable law or agreed to in writing, software |
|
11 | * distributed under the License is distributed on an "AS IS" BASIS, |
|
12 | * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. |
|
13 | * See the License for the specific language governing permissions and |
|
14 | * limitations under the License. |
|
15 | */ |
|
16 | package org.apache.commons.math.analysis; |
|
17 | ||
18 | import java.io.Serializable; |
|
19 | import org.apache.commons.math.MathException; |
|
20 | ||
21 | /** |
|
22 | * Implements the <a href=" |
|
23 | * "http://mathworld.wolfram.com/NewtonsDividedDifferenceInterpolationFormula.html"> |
|
24 | * Divided Difference Algorithm</a> for interpolation of real univariate |
|
25 | * functions. For reference, see <b>Introduction to Numerical Analysis</b>, |
|
26 | * ISBN 038795452X, chapter 2. |
|
27 | * <p> |
|
28 | * The actual code of Neville's evalution is in PolynomialFunctionLagrangeForm, |
|
29 | * this class provides an easy-to-use interface to it. |
|
30 | * |
|
31 | * @version $Revision$ $Date$ |
|
32 | */ |
|
33 | 6 | public class DividedDifferenceInterpolator implements UnivariateRealInterpolator, |
34 | Serializable { |
|
35 | ||
36 | /** serializable version identifier */ |
|
37 | static final long serialVersionUID = 107049519551235069L; |
|
38 | ||
39 | /** |
|
40 | * Computes an interpolating function for the data set. |
|
41 | * |
|
42 | * @param x the interpolating points array |
|
43 | * @param y the interpolating values array |
|
44 | * @return a function which interpolates the data set |
|
45 | * @throws MathException if arguments are invalid |
|
46 | */ |
|
47 | public UnivariateRealFunction interpolate(double x[], double y[]) throws |
|
48 | MathException { |
|
49 | ||
50 | /** |
|
51 | * a[] and c[] are defined in the general formula of Newton form: |
|
52 | * p(x) = a[0] + a[1](x-c[0]) + a[2](x-c[0])(x-c[1]) + ... + |
|
53 | * a[n](x-c[0])(x-c[1])...(x-c[n-1]) |
|
54 | */ |
|
55 | double a[], c[]; |
|
56 | ||
57 | 6 | PolynomialFunctionLagrangeForm.verifyInterpolationArray(x, y); |
58 | ||
59 | /** |
|
60 | * When used for interpolation, the Newton form formula becomes |
|
61 | * p(x) = f[x0] + f[x0,x1](x-x0) + f[x0,x1,x2](x-x0)(x-x1) + ... + |
|
62 | * f[x0,x1,...,x[n-1]](x-x0)(x-x1)...(x-x[n-2]) |
|
63 | * Therefore, a[k] = f[x0,x1,...,xk], c[k] = x[k]. |
|
64 | * <p> |
|
65 | * Note x[], y[], a[] have the same length but c[]'s size is one less. |
|
66 | */ |
|
67 | 6 | c = new double[x.length-1]; |
68 | 30 | for (int i = 0; i < c.length; i++) { |
69 | 24 | c[i] = x[i]; |
70 | } |
|
71 | 6 | a = computeDividedDifference(x, y); |
72 | ||
73 | PolynomialFunctionNewtonForm p; |
|
74 | 4 | p = new PolynomialFunctionNewtonForm(a, c); |
75 | 4 | return p; |
76 | } |
|
77 | ||
78 | /** |
|
79 | * Returns a copy of the divided difference array. |
|
80 | * <p> |
|
81 | * The divided difference array is defined recursively by <pre> |
|
82 | * f[x0] = f(x0) |
|
83 | * f[x0,x1,...,xk] = (f(x1,...,xk) - f(x0,...,x[k-1])) / (xk - x0) |
|
84 | * </pre><p> |
|
85 | * The computational complexity is O(N^2). |
|
86 | * |
|
87 | * @return a fresh copy of the divided difference array |
|
88 | * @throws MathException if any abscissas coincide |
|
89 | */ |
|
90 | protected static double[] computeDividedDifference(double x[], double y[]) |
|
91 | throws MathException { |
|
92 | ||
93 | int i, j, n; |
|
94 | double divdiff[], a[], denominator; |
|
95 | ||
96 | 6 | PolynomialFunctionLagrangeForm.verifyInterpolationArray(x, y); |
97 | ||
98 | 6 | n = x.length; |
99 | 6 | divdiff = new double[n]; |
100 | 36 | for (i = 0; i < n; i++) { |
101 | 30 | divdiff[i] = y[i]; // initialization |
102 | } |
|
103 | ||
104 | 6 | a = new double [n]; |
105 | 6 | a[0] = divdiff[0]; |
106 | 24 | for (i = 1; i < n; i++) { |
107 | 72 | for (j = 0; j < n-i; j++) { |
108 | 54 | denominator = x[j+i] - x[j]; |
109 | 54 | if (denominator == 0.0) { |
110 | // This happens only when two abscissas are identical. |
|
111 | 2 | throw new MathException |
112 | ("Identical abscissas cause division by zero."); |
|
113 | } |
|
114 | 52 | divdiff[j] = (divdiff[j+1] - divdiff[j]) / denominator; |
115 | } |
|
116 | 18 | a[i] = divdiff[0]; |
117 | } |
|
118 | ||
119 | 4 | return a; |
120 | } |
|
121 | } |