Classes in this File | Line Coverage | Branch Coverage | Complexity | ||||||||
PolynomialFunctionLagrangeForm |
|
| 3.3333333333333335;3.333 |
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.FunctionEvaluationException; |
|
20 | ||
21 | /** |
|
22 | * Implements the representation of a real polynomial function in |
|
23 | * <a href="http://mathworld.wolfram.com/LagrangeInterpolatingPolynomial.html"> |
|
24 | * Lagrange Form</a>. For reference, see <b>Introduction to Numerical |
|
25 | * Analysis</b>, ISBN 038795452X, chapter 2. |
|
26 | * <p> |
|
27 | * The approximated function should be smooth enough for Lagrange polynomial |
|
28 | * to work well. Otherwise, consider using splines instead. |
|
29 | * |
|
30 | * @version $Revision$ $Date$ |
|
31 | */ |
|
32 | public class PolynomialFunctionLagrangeForm implements UnivariateRealFunction, |
|
33 | Serializable { |
|
34 | ||
35 | /** serializable version identifier */ |
|
36 | static final long serialVersionUID = -3965199246151093920L; |
|
37 | ||
38 | /** |
|
39 | * The coefficients of the polynomial, ordered by degree -- i.e. |
|
40 | * coefficients[0] is the constant term and coefficients[n] is the |
|
41 | * coefficient of x^n where n is the degree of the polynomial. |
|
42 | */ |
|
43 | private double coefficients[]; |
|
44 | ||
45 | /** |
|
46 | * Interpolating points (abscissas) and the function values at these points. |
|
47 | */ |
|
48 | private double x[], y[]; |
|
49 | ||
50 | /** |
|
51 | * Whether the polynomial coefficients are available. |
|
52 | */ |
|
53 | private boolean coefficientsComputed; |
|
54 | ||
55 | /** |
|
56 | * Construct a Lagrange polynomial with the given abscissas and function |
|
57 | * values. The order of interpolating points are not important. |
|
58 | * <p> |
|
59 | * The constructor makes copy of the input arrays and assigns them. |
|
60 | * |
|
61 | * @param x interpolating points |
|
62 | * @param y function values at interpolating points |
|
63 | * @throws IllegalArgumentException if input arrays are not valid |
|
64 | */ |
|
65 | PolynomialFunctionLagrangeForm(double x[], double y[]) throws |
|
66 | 16 | IllegalArgumentException { |
67 | ||
68 | 16 | verifyInterpolationArray(x, y); |
69 | 12 | this.x = new double[x.length]; |
70 | 12 | this.y = new double[y.length]; |
71 | 12 | System.arraycopy(x, 0, this.x, 0, x.length); |
72 | 12 | System.arraycopy(y, 0, this.y, 0, y.length); |
73 | 12 | coefficientsComputed = false; |
74 | 12 | } |
75 | ||
76 | /** |
|
77 | * Calculate the function value at the given point. |
|
78 | * |
|
79 | * @param z the point at which the function value is to be computed |
|
80 | * @return the function value |
|
81 | * @throws FunctionEvaluationException if a runtime error occurs |
|
82 | * @see UnivariateRealFunction#value(double) |
|
83 | */ |
|
84 | public double value(double z) throws FunctionEvaluationException { |
|
85 | 30 | return evaluate(x, y, z); |
86 | } |
|
87 | ||
88 | /** |
|
89 | * Returns the degree of the polynomial. |
|
90 | * |
|
91 | * @return the degree of the polynomial |
|
92 | */ |
|
93 | public int degree() { |
|
94 | 12 | return x.length - 1; |
95 | } |
|
96 | ||
97 | /** |
|
98 | * Returns a copy of the interpolating points array. |
|
99 | * <p> |
|
100 | * Changes made to the returned copy will not affect the polynomial. |
|
101 | * |
|
102 | * @return a fresh copy of the interpolating points array |
|
103 | */ |
|
104 | public double[] getInterpolatingPoints() { |
|
105 | 0 | double[] out = new double[x.length]; |
106 | 0 | System.arraycopy(x, 0, out, 0, x.length); |
107 | 0 | return out; |
108 | } |
|
109 | ||
110 | /** |
|
111 | * Returns a copy of the interpolating values array. |
|
112 | * <p> |
|
113 | * Changes made to the returned copy will not affect the polynomial. |
|
114 | * |
|
115 | * @return a fresh copy of the interpolating values array |
|
116 | */ |
|
117 | public double[] getInterpolatingValues() { |
|
118 | 0 | double[] out = new double[y.length]; |
119 | 0 | System.arraycopy(y, 0, out, 0, y.length); |
120 | 0 | return out; |
121 | } |
|
122 | ||
123 | /** |
|
124 | * Returns a copy of the coefficients array. |
|
125 | * <p> |
|
126 | * Changes made to the returned copy will not affect the polynomial. |
|
127 | * |
|
128 | * @return a fresh copy of the coefficients array |
|
129 | */ |
|
130 | public double[] getCoefficients() { |
|
131 | 6 | if (!coefficientsComputed) { |
132 | 6 | computeCoefficients(); |
133 | } |
|
134 | 6 | double[] out = new double[coefficients.length]; |
135 | 6 | System.arraycopy(coefficients, 0, out, 0, coefficients.length); |
136 | 6 | return out; |
137 | } |
|
138 | ||
139 | /** |
|
140 | * Evaluate the Lagrange polynomial using |
|
141 | * <a href="http://mathworld.wolfram.com/NevillesAlgorithm.html"> |
|
142 | * Neville's Algorithm</a>. It takes O(N^2) time. |
|
143 | * <p> |
|
144 | * This function is made public static so that users can call it directly |
|
145 | * without instantiating PolynomialFunctionLagrangeForm object. |
|
146 | * |
|
147 | * @param x the interpolating points array |
|
148 | * @param y the interpolating values array |
|
149 | * @param z the point at which the function value is to be computed |
|
150 | * @return the function value |
|
151 | * @throws FunctionEvaluationException if a runtime error occurs |
|
152 | * @throws IllegalArgumentException if inputs are not valid |
|
153 | */ |
|
154 | public static double evaluate(double x[], double y[], double z) throws |
|
155 | FunctionEvaluationException, IllegalArgumentException { |
|
156 | ||
157 | 30 | int i, j, n, nearest = 0; |
158 | double value, c[], d[], tc, td, divider, w, dist, min_dist; |
|
159 | ||
160 | 30 | verifyInterpolationArray(x, y); |
161 | ||
162 | 30 | n = x.length; |
163 | 30 | c = new double[n]; |
164 | 30 | d = new double[n]; |
165 | 30 | min_dist = Double.POSITIVE_INFINITY; |
166 | 158 | for (i = 0; i < n; i++) { |
167 | // initialize the difference arrays |
|
168 | 128 | c[i] = y[i]; |
169 | 128 | d[i] = y[i]; |
170 | // find out the abscissa closest to z |
|
171 | 128 | dist = Math.abs(z - x[i]); |
172 | 128 | if (dist < min_dist) { |
173 | 74 | nearest = i; |
174 | 74 | min_dist = dist; |
175 | } |
|
176 | } |
|
177 | ||
178 | // initial approximation to the function value at z |
|
179 | 30 | value = y[nearest]; |
180 | ||
181 | 122 | for (i = 1; i < n; i++) { |
182 | 330 | for (j = 0; j < n-i; j++) { |
183 | 238 | tc = x[j] - z; |
184 | 238 | td = x[i+j] - z; |
185 | 238 | divider = x[j] - x[i+j]; |
186 | 238 | if (divider == 0.0) { |
187 | // This happens only when two abscissas are identical. |
|
188 | 2 | throw new FunctionEvaluationException(z, |
189 | "Identical abscissas cause division by zero: x[" + |
|
190 | i + "] = x[" + (i+j) + "] = " + x[i]); |
|
191 | } |
|
192 | // update the difference arrays |
|
193 | 236 | w = (c[j+1] - d[j]) / divider; |
194 | 236 | c[j] = tc * w; |
195 | 236 | d[j] = td * w; |
196 | } |
|
197 | // sum up the difference terms to get the final value |
|
198 | 92 | if (nearest < 0.5*(n-i+1)) { |
199 | 34 | value += c[nearest]; // fork down |
200 | } else { |
|
201 | 58 | nearest--; |
202 | 58 | value += d[nearest]; // fork up |
203 | } |
|
204 | } |
|
205 | ||
206 | 28 | return value; |
207 | } |
|
208 | ||
209 | /** |
|
210 | * Calculate the coefficients of Lagrange polynomial from the |
|
211 | * interpolation data. It takes O(N^2) time. |
|
212 | * <p> |
|
213 | * Note this computation can be ill-conditioned. Use with caution |
|
214 | * and only when it is necessary. |
|
215 | * |
|
216 | * @throws ArithmeticException if any abscissas coincide |
|
217 | */ |
|
218 | protected void computeCoefficients() throws ArithmeticException { |
|
219 | int i, j, n; |
|
220 | double c[], tc[], d, t; |
|
221 | ||
222 | 6 | n = degree() + 1; |
223 | 6 | coefficients = new double[n]; |
224 | 28 | for (i = 0; i < n; i++) { |
225 | 22 | coefficients[i] = 0.0; |
226 | } |
|
227 | ||
228 | // c[] are the coefficients of P(x) = (x-x[0])(x-x[1])...(x-x[n-1]) |
|
229 | 6 | c = new double[n+1]; |
230 | 6 | c[0] = 1.0; |
231 | 28 | for (i = 0; i < n; i++) { |
232 | 60 | for (j = i; j > 0; j--) { |
233 | 38 | c[j] = c[j-1] - c[j] * x[i]; |
234 | } |
|
235 | 22 | c[0] *= (-x[i]); |
236 | 22 | c[i+1] = 1; |
237 | } |
|
238 | ||
239 | 6 | tc = new double[n]; |
240 | 28 | for (i = 0; i < n; i++) { |
241 | // d = (x[i]-x[0])...(x[i]-x[i-1])(x[i]-x[i+1])...(x[i]-x[n-1]) |
|
242 | 22 | d = 1; |
243 | 120 | for (j = 0; j < n; j++) { |
244 | 98 | if (i != j) { |
245 | 76 | d *= (x[i] - x[j]); |
246 | } |
|
247 | } |
|
248 | 22 | if (d == 0.0) { |
249 | // This happens only when two abscissas are identical. |
|
250 | 0 | throw new ArithmeticException |
251 | ("Identical abscissas cause division by zero."); |
|
252 | } |
|
253 | 22 | t = y[i] / d; |
254 | // Lagrange polynomial is the sum of n terms, each of which is a |
|
255 | // polynomial of degree n-1. tc[] are the coefficients of the i-th |
|
256 | // numerator Pi(x) = (x-x[0])...(x-x[i-1])(x-x[i+1])...(x-x[n-1]). |
|
257 | 22 | tc[n-1] = c[n]; // actually c[n] = 1 |
258 | 22 | coefficients[n-1] += t * tc[n-1]; |
259 | 98 | for (j = n-2; j >= 0; j--) { |
260 | 76 | tc[j] = c[j+1] + tc[j+1] * x[i]; |
261 | 76 | coefficients[j] += t * tc[j]; |
262 | } |
|
263 | } |
|
264 | ||
265 | 6 | coefficientsComputed = true; |
266 | 6 | } |
267 | ||
268 | /** |
|
269 | * Verifies that the interpolation arrays are valid. |
|
270 | * <p> |
|
271 | * The interpolating points must be distinct. However it is not |
|
272 | * verified here, it is checked in evaluate() and computeCoefficients(). |
|
273 | * |
|
274 | * @throws IllegalArgumentException if not valid |
|
275 | * @see #evaluate(double[], double[], double) |
|
276 | * @see #computeCoefficients() |
|
277 | */ |
|
278 | protected static void verifyInterpolationArray(double x[], double y[]) throws |
|
279 | IllegalArgumentException { |
|
280 | ||
281 | 58 | if (x.length < 2 || y.length < 2) { |
282 | 2 | throw new IllegalArgumentException |
283 | ("Interpolation requires at least two points."); |
|
284 | } |
|
285 | 56 | if (x.length != y.length) { |
286 | 2 | throw new IllegalArgumentException |
287 | ("Abscissa and value arrays must have the same length."); |
|
288 | } |
|
289 | 54 | } |
290 | } |