← Back to branch summary

~ubuntu-branches/ubuntu/quantal/commons-math/quantal

« back to all changes in this revision

Viewing changes to src/main/java/org/apache/commons/math/analysis/polynomials/PolynomialFunctionLagrangeForm.java

  • Committer: Bazaar Package Importer
  • Author(s): Damien Raude-Morvan
  • Date: 2010-04-05 23:33:02 UTC
  • mfrom: (1.1.2 upstream)
  • Revision ID: james.westby@ubuntu.com-20100405233302-gpqlceked76nw28a
Tags: 2.1-1
* New upstream release.
* Bump Standards-Version to 3.8.4: no changes needed
* Bump debhelper to >= 7
* Switch to 3.0 (quilt) source format:
  - Remove B-D on quilt
  - Add d/source/format
  - Remove d/README.source

Show diffs side-by-side

added added

removed removed

Lines of Context:
src/main/java/org/apache/commons/math/analysis/polynomials/PolynomialFunctionLagrangeForm.java
30
30
 * The approximated function should be smooth enough for Lagrange polynomial
31
31
 * to work well. Otherwise, consider using splines instead.</p>
32
32
 *
33
 
 * @version $Revision: 799857 $ $Date: 2009-08-01 09:07:12 -0400 (Sat, 01 Aug 2009) $
 
33
 * @version $Revision: 922708 $ $Date: 2010-03-13 20:15:47 -0500 (Sat, 13 Mar 2010) $
34
34
 * @since 1.2
35
35
 */
36
36
public class PolynomialFunctionLagrangeForm implements UnivariateRealFunction {
37
37
 
38
38
    /**
39
39
     * The coefficients of the polynomial, ordered by degree -- i.e.
40
 
     * coefficients[0] is the constant term and coefficients[n] is the 
 
40
     * coefficients[0] is the constant term and coefficients[n] is the
41
41
     * coefficient of x^n where n is the degree of the polynomial.
42
42
     */
43
43
    private double coefficients[];
44
44
 
45
45
    /**
46
 
     * Interpolating points (abscissas) and the function values at these points.
47
 
     */
48
 
    private double x[], y[];
 
46
     * Interpolating points (abscissas).
 
47
     */
 
48
    private final double x[];
 
49
 
 
50
    /**
 
51
     * Function values at interpolating points.
 
52
     */
 
53
    private final double y[];
49
54
 
50
55
    /**
51
56
     * Whether the polynomial coefficients are available.
57
62
     * values. The order of interpolating points are not important.
58
63
     * <p>
59
64
     * The constructor makes copy of the input arrays and assigns them.</p>
60
 
     * 
 
65
     *
61
66
     * @param x interpolating points
62
67
     * @param y function values at interpolating points
63
68
     * @throws IllegalArgumentException if input arrays are not valid
91
96
 
92
97
    /**
93
98
     * Returns the degree of the polynomial.
94
 
     * 
 
99
     *
95
100
     * @return the degree of the polynomial
96
101
     */
97
102
    public int degree() {
102
107
     * Returns a copy of the interpolating points array.
103
108
     * <p>
104
109
     * Changes made to the returned copy will not affect the polynomial.</p>
105
 
     * 
 
110
     *
106
111
     * @return a fresh copy of the interpolating points array
107
112
     */
108
113
    public double[] getInterpolatingPoints() {
115
120
     * Returns a copy of the interpolating values array.
116
121
     * <p>
117
122
     * Changes made to the returned copy will not affect the polynomial.</p>
118
 
     * 
 
123
     *
119
124
     * @return a fresh copy of the interpolating values array
120
125
     */
121
126
    public double[] getInterpolatingValues() {
128
133
     * Returns a copy of the coefficients array.
129
134
     * <p>
130
135
     * Changes made to the returned copy will not affect the polynomial.</p>
131
 
     * 
 
136
     * <p>
 
137
     * Note that coefficients computation can be ill-conditioned. Use with caution
 
138
     * and only when it is necessary.</p>
 
139
     *
132
140
     * @return a fresh copy of the coefficients array
133
141
     */
134
142
    public double[] getCoefficients() {
141
149
    }
142
150
 
143
151
    /**
144
 
     * Evaluate the Lagrange polynomial using 
 
152
     * Evaluate the Lagrange polynomial using
145
153
     * <a href="http://mathworld.wolfram.com/NevillesAlgorithm.html">
146
154
     * Neville's Algorithm</a>. It takes O(N^2) time.
147
155
     * <p>
158
166
    public static double evaluate(double x[], double y[], double z) throws
159
167
        DuplicateSampleAbscissaException, IllegalArgumentException {
160
168
 
161
 
        int i, j, n, nearest = 0;
162
 
        double value, c[], d[], tc, td, divider, w, dist, min_dist;
163
 
 
164
169
        verifyInterpolationArray(x, y);
165
170
 
166
 
        n = x.length;
167
 
        c = new double[n];
168
 
        d = new double[n];
169
 
        min_dist = Double.POSITIVE_INFINITY;
170
 
        for (i = 0; i < n; i++) {
 
171
        int nearest = 0;
 
172
        final int n = x.length;
 
173
        final double[] c = new double[n];
 
174
        final double[] d = new double[n];
 
175
        double min_dist = Double.POSITIVE_INFINITY;
 
176
        for (int i = 0; i < n; i++) {
171
177
            // initialize the difference arrays
172
178
            c[i] = y[i];
173
179
            d[i] = y[i];
174
180
            // find out the abscissa closest to z
175
 
            dist = Math.abs(z - x[i]);
 
181
            final double dist = Math.abs(z - x[i]);
176
182
            if (dist < min_dist) {
177
183
                nearest = i;
178
184
                min_dist = dist;
180
186
        }
181
187
 
182
188
        // initial approximation to the function value at z
183
 
        value = y[nearest];
 
189
        double value = y[nearest];
184
190
 
185
 
        for (i = 1; i < n; i++) {
186
 
            for (j = 0; j < n-i; j++) {
187
 
                tc = x[j] - z;
188
 
                td = x[i+j] - z;
189
 
                divider = x[j] - x[i+j];
 
191
        for (int i = 1; i < n; i++) {
 
192
            for (int j = 0; j < n-i; j++) {
 
193
                final double tc = x[j] - z;
 
194
                final double td = x[i+j] - z;
 
195
                final double divider = x[j] - x[i+j];
190
196
                if (divider == 0.0) {
191
197
                    // This happens only when two abscissas are identical.
192
198
                    throw new DuplicateSampleAbscissaException(x[i], i, i+j);
193
199
                }
194
200
                // update the difference arrays
195
 
                w = (c[j+1] - d[j]) / divider;
 
201
                final double w = (c[j+1] - d[j]) / divider;
196
202
                c[j] = tc * w;
197
203
                d[j] = td * w;
198
204
            }
218
224
     * @throws ArithmeticException if any abscissas coincide
219
225
     */
220
226
    protected void computeCoefficients() throws ArithmeticException {
221
 
        int i, j, n;
222
 
        double c[], tc[], d, t;
223
227
 
224
 
        n = degree() + 1;
 
228
        final int n = degree() + 1;
225
229
        coefficients = new double[n];
226
 
        for (i = 0; i < n; i++) {
 
230
        for (int i = 0; i < n; i++) {
227
231
            coefficients[i] = 0.0;
228
232
        }
229
233
 
230
234
        // c[] are the coefficients of P(x) = (x-x[0])(x-x[1])...(x-x[n-1])
231
 
        c = new double[n+1];
 
235
        final double[] c = new double[n+1];
232
236
        c[0] = 1.0;
233
 
        for (i = 0; i < n; i++) {
234
 
            for (j = i; j > 0; j--) {
 
237
        for (int i = 0; i < n; i++) {
 
238
            for (int j = i; j > 0; j--) {
235
239
                c[j] = c[j-1] - c[j] * x[i];
236
240
            }
237
 
            c[0] *= (-x[i]);
 
241
            c[0] *= -x[i];
238
242
            c[i+1] = 1;
239
243
        }
240
244
 
241
 
        tc = new double[n];
242
 
        for (i = 0; i < n; i++) {
 
245
        final double[] tc = new double[n];
 
246
        for (int i = 0; i < n; i++) {
243
247
            // d = (x[i]-x[0])...(x[i]-x[i-1])(x[i]-x[i+1])...(x[i]-x[n-1])
244
 
            d = 1;
245
 
            for (j = 0; j < n; j++) {
 
248
            double d = 1;
 
249
            for (int j = 0; j < n; j++) {
246
250
                if (i != j) {
247
 
                    d *= (x[i] - x[j]);
 
251
                    d *= x[i] - x[j];
248
252
                }
249
253
            }
250
254
            if (d == 0.0) {
256
260
                    }
257
261
                }
258
262
            }
259
 
            t = y[i] / d;
 
263
            final double t = y[i] / d;
260
264
            // Lagrange polynomial is the sum of n terms, each of which is a
261
265
            // polynomial of degree n-1. tc[] are the coefficients of the i-th
262
266
            // numerator Pi(x) = (x-x[0])...(x-x[i-1])(x-x[i+1])...(x-x[n-1]).
263
267
            tc[n-1] = c[n];     // actually c[n] = 1
264
268
            coefficients[n-1] += t * tc[n-1];
265
 
            for (j = n-2; j >= 0; j--) {
 
269
            for (int j = n-2; j >= 0; j--) {
266
270
                tc[j] = c[j+1] + tc[j+1] * x[i];
267
271
                coefficients[j] += t * tc[j];
268
272
            }
274
278
    /**
275
279
     * Verifies that the interpolation arrays are valid.
276
280
     * <p>
 
281
     * The arrays features checked by this method are that both arrays have the
 
282
     * same length and this length is at least 2.
 
283
     * </p>
 
284
     * <p>
277
285
     * The interpolating points must be distinct. However it is not
278
 
     * verified here, it is checked in evaluate() and computeCoefficients().</p>
279
 
     * 
 
286
     * verified here, it is checked in evaluate() and computeCoefficients().
 
287
     * </p>
 
288
     *
280
289
     * @param x the interpolating points array
281
290
     * @param y the interpolating values array
282
291
     * @throws IllegalArgumentException if not valid
283
292
     * @see #evaluate(double[], double[], double)
284
293
     * @see #computeCoefficients()
285
294
     */
286
 
    public static void verifyInterpolationArray(double x[], double y[]) throws
287
 
        IllegalArgumentException {
 
295
    public static void verifyInterpolationArray(double x[], double y[])
 
296
        throws IllegalArgumentException {
288
297
 
289
 
        if (Math.min(x.length, y.length) < 2) {
290
 
            throw MathRuntimeException.createIllegalArgumentException(
291
 
                  "{0} points are required, got only {1}",
292
 
                  2, Math.min(x.length, y.length));
293
 
        }
294
298
        if (x.length != y.length) {
295
299
            throw MathRuntimeException.createIllegalArgumentException(
296
300
                  "dimension mismatch {0} != {1}", x.length, y.length);
297
301
        }
 
302
 
 
303
        if (x.length < 2) {
 
304
            throw MathRuntimeException.createIllegalArgumentException(
 
305
                  "{0} points are required, got only {1}", 2, x.length);
 
306
        }
 
307
 
298
308
    }
299
309
}