← Back to branch summary

~ubuntu-branches/ubuntu/maverick/commons-math/maverick

« back to all changes in this revision

Viewing changes to src/main/java/org/apache/commons/math/analysis/integration/LegendreGaussIntegrator.java

  • Committer: Bazaar Package Importer
  • Author(s): Damien Raude-Morvan
  • Date: 2009-08-22 01:13:25 UTC
  • mfrom: (1.1.1 upstream)
  • Revision ID: james.westby@ubuntu.com-20090822011325-hi4peq1ua5weguwn
Tags: 2.0-1
* New upstream release.
* Set Maintainer field to Debian Java Team
* Add myself as Uploaders
* Switch to Quilt patch system:
  - Refresh all patchs
  - Remove B-D on dpatch, Add B-D on quilt
  - Include patchsys-quilt.mk in debian/rules
* Bump Standards-Version to 3.8.3:
  - Add a README.source to describe patch system
* Maven POMs:
  - Add a Build-Depends-Indep dependency on maven-repo-helper
  - Use mh_installpom and mh_installjar to install the POM and the jar to the
    Maven repository
* Use default-jdk/jre:
  - Depends on java5-runtime-headless
  - Build-Depends on default-jdk
  - Use /usr/lib/jvm/default-java as JAVA_HOME
* Move api documentation to /usr/share/doc/libcommons-math-java/api
* Build-Depends on junit4 instead of junit

Show diffs side-by-side

added added

removed removed

Lines of Context:
src/main/java/org/apache/commons/math/analysis/integration/LegendreGaussIntegrator.java
 
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.math.analysis.integration;
 
18
 
 
19
import org.apache.commons.math.ConvergenceException;
 
20
import org.apache.commons.math.FunctionEvaluationException;
 
21
import org.apache.commons.math.MathRuntimeException;
 
22
import org.apache.commons.math.MaxIterationsExceededException;
 
23
import org.apache.commons.math.analysis.UnivariateRealFunction;
 
24
 
 
25
/**
 
26
 * Implements the <a href="http://mathworld.wolfram.com/Legendre-GaussQuadrature.html">
 
27
 * Legendre-Gauss</a> quadrature formula.
 
28
 * <p>
 
29
 * Legendre-Gauss integrators are efficient integrators that can
 
30
 * accurately integrate functions with few functions evaluations. A
 
31
 * Legendre-Gauss integrator using an n-points quadrature formula can
 
32
 * integrate exactly 2n-1 degree polynomialss.
 
33
 * </p>
 
34
 * <p>
 
35
 * These integrators evaluate the function on n carefully chosen
 
36
 * abscissas in each step interval (mapped to the canonical [-1  1] interval).
 
37
 * The evaluation abscissas are not evenly spaced and none of them are
 
38
 * at the interval endpoints. This implies the function integrated can be
 
39
 * undefined at integration interval endpoints.
 
40
 * </p>
 
41
 * <p>
 
42
 * The evaluation abscissas x<sub>i</sub> are the roots of the degree n
 
43
 * Legendre polynomial. The weights a<sub>i</sub> of the quadrature formula
 
44
 * integrals from -1 to +1 &int; Li<sup>2</sup> where Li (x) =
 
45
 * &prod; (x-x<sub>k</sub>)/(x<sub>i</sub>-x<sub>k</sub>) for k != i.
 
46
 * </p>
 
47
 * <p>
 
48
 * @version $Revision: 799857 $ $Date: 2009-08-01 09:07:12 -0400 (Sat, 01 Aug 2009) $
 
49
 * @since 1.2
 
50
 */
 
51
 
 
52
public class LegendreGaussIntegrator extends UnivariateRealIntegratorImpl {
 
53
 
 
54
    /** Abscissas for the 2 points method. */
 
55
    private static final double[] ABSCISSAS_2 = {
 
56
        -1.0 / Math.sqrt(3.0),
 
57
         1.0 / Math.sqrt(3.0)
 
58
    };
 
59
 
 
60
    /** Weights for the 2 points method. */
 
61
    private static final double[] WEIGHTS_2 = {
 
62
        1.0,
 
63
        1.0
 
64
    };
 
65
 
 
66
    /** Abscissas for the 3 points method. */
 
67
    private static final double[] ABSCISSAS_3 = {
 
68
        -Math.sqrt(0.6),
 
69
         0.0,
 
70
         Math.sqrt(0.6)
 
71
    };
 
72
 
 
73
    /** Weights for the 3 points method. */
 
74
    private static final double[] WEIGHTS_3 = {
 
75
        5.0 / 9.0,
 
76
        8.0 / 9.0,
 
77
        5.0 / 9.0
 
78
    };
 
79
 
 
80
    /** Abscissas for the 4 points method. */
 
81
    private static final double[] ABSCISSAS_4 = {
 
82
        -Math.sqrt((15.0 + 2.0 * Math.sqrt(30.0)) / 35.0),
 
83
        -Math.sqrt((15.0 - 2.0 * Math.sqrt(30.0)) / 35.0),
 
84
         Math.sqrt((15.0 - 2.0 * Math.sqrt(30.0)) / 35.0),
 
85
         Math.sqrt((15.0 + 2.0 * Math.sqrt(30.0)) / 35.0)
 
86
    };
 
87
 
 
88
    /** Weights for the 4 points method. */
 
89
    private static final double[] WEIGHTS_4 = {
 
90
        (90.0 - 5.0 * Math.sqrt(30.0)) / 180.0,
 
91
        (90.0 + 5.0 * Math.sqrt(30.0)) / 180.0,
 
92
        (90.0 + 5.0 * Math.sqrt(30.0)) / 180.0,
 
93
        (90.0 - 5.0 * Math.sqrt(30.0)) / 180.0
 
94
    };
 
95
 
 
96
    /** Abscissas for the 5 points method. */
 
97
    private static final double[] ABSCISSAS_5 = {
 
98
        -Math.sqrt((35.0 + 2.0 * Math.sqrt(70.0)) / 63.0),
 
99
        -Math.sqrt((35.0 - 2.0 * Math.sqrt(70.0)) / 63.0),
 
100
         0.0,
 
101
         Math.sqrt((35.0 - 2.0 * Math.sqrt(70.0)) / 63.0),
 
102
         Math.sqrt((35.0 + 2.0 * Math.sqrt(70.0)) / 63.0)
 
103
    };
 
104
 
 
105
    /** Weights for the 5 points method. */
 
106
    private static final double[] WEIGHTS_5 = {
 
107
        (322.0 - 13.0 * Math.sqrt(70.0)) / 900.0,
 
108
        (322.0 + 13.0 * Math.sqrt(70.0)) / 900.0,
 
109
        128.0 / 225.0,
 
110
        (322.0 + 13.0 * Math.sqrt(70.0)) / 900.0,
 
111
        (322.0 - 13.0 * Math.sqrt(70.0)) / 900.0
 
112
    };
 
113
 
 
114
    /** Abscissas for the current method. */
 
115
    private final double[] abscissas;
 
116
 
 
117
    /** Weights for the current method. */
 
118
    private final double[] weights;
 
119
 
 
120
    /** Build a Legendre-Gauss integrator.
 
121
     * @param n number of points desired (must be between 2 and 5 inclusive)
 
122
     * @param defaultMaximalIterationCount maximum number of iterations
 
123
     * @exception IllegalArgumentException if the number of points is not
 
124
     * in the supported range
 
125
     */
 
126
    public LegendreGaussIntegrator(final int n, final int defaultMaximalIterationCount)
 
127
        throws IllegalArgumentException {
 
128
        super(defaultMaximalIterationCount);
 
129
        switch(n) {
 
130
        case 2 :
 
131
            abscissas = ABSCISSAS_2;
 
132
            weights   = WEIGHTS_2;
 
133
            break;
 
134
        case 3 :
 
135
            abscissas = ABSCISSAS_3;
 
136
            weights   = WEIGHTS_3;
 
137
            break;
 
138
        case 4 :
 
139
            abscissas = ABSCISSAS_4;
 
140
            weights   = WEIGHTS_4;
 
141
            break;
 
142
        case 5 :
 
143
            abscissas = ABSCISSAS_5;
 
144
            weights   = WEIGHTS_5;
 
145
            break;
 
146
        default :
 
147
            throw MathRuntimeException.createIllegalArgumentException(
 
148
                    "{0} points Legendre-Gauss integrator not supported, " +
 
149
                    "number of points must be in the {1}-{2} range",
 
150
                    n, 2, 5);
 
151
        }
 
152
 
 
153
    }
 
154
 
 
155
    /** {@inheritDoc} */
 
156
    @Deprecated
 
157
    public double integrate(final double min, final double max)
 
158
        throws ConvergenceException,  FunctionEvaluationException, IllegalArgumentException {
 
159
        return integrate(f, min, max);
 
160
    }
 
161
 
 
162
    /** {@inheritDoc} */
 
163
    public double integrate(final UnivariateRealFunction f,
 
164
            final double min, final double max)
 
165
        throws ConvergenceException,  FunctionEvaluationException, IllegalArgumentException {
 
166
        
 
167
        clearResult();
 
168
        verifyInterval(min, max);
 
169
        verifyIterationCount();
 
170
 
 
171
        // compute first estimate with a single step
 
172
        double oldt = stage(f, min, max, 1);
 
173
 
 
174
        int n = 2;
 
175
        for (int i = 0; i < maximalIterationCount; ++i) {
 
176
 
 
177
            // improve integral with a larger number of steps
 
178
            final double t = stage(f, min, max, n);
 
179
 
 
180
            // estimate error
 
181
            final double delta = Math.abs(t - oldt);
 
182
            final double limit =
 
183
                Math.max(absoluteAccuracy,
 
184
                         relativeAccuracy * (Math.abs(oldt) + Math.abs(t)) * 0.5);
 
185
 
 
186
            // check convergence
 
187
            if ((i + 1 >= minimalIterationCount) && (delta <= limit)) {
 
188
                setResult(t, i);
 
189
                return result;
 
190
            }
 
191
 
 
192
            // prepare next iteration
 
193
            double ratio = Math.min(4, Math.pow(delta / limit, 0.5 / abscissas.length));
 
194
            n = Math.max((int) (ratio * n), n + 1);
 
195
            oldt = t;
 
196
 
 
197
        }
 
198
 
 
199
        throw new MaxIterationsExceededException(maximalIterationCount);
 
200
 
 
201
    }
 
202
 
 
203
    /**
 
204
     * Compute the n-th stage integral.
 
205
     * @param f the integrand function
 
206
     * @param min the lower bound for the interval
 
207
     * @param max the upper bound for the interval
 
208
     * @param n number of steps
 
209
     * @return the value of n-th stage integral
 
210
     * @throws FunctionEvaluationException if an error occurs evaluating the
 
211
     * function
 
212
     */
 
213
    private double stage(final UnivariateRealFunction f,
 
214
                         final double min, final double max, final int n)
 
215
        throws FunctionEvaluationException {
 
216
 
 
217
        // set up the step for the current stage
 
218
        final double step     = (max - min) / n;
 
219
        final double halfStep = step / 2.0;
 
220
 
 
221
        // integrate over all elementary steps
 
222
        double midPoint = min + halfStep;
 
223
        double sum = 0.0;
 
224
        for (int i = 0; i < n; ++i) {
 
225
            for (int j = 0; j < abscissas.length; ++j) {
 
226
                sum += weights[j] * f.value(midPoint + halfStep * abscissas[j]);
 
227
            }
 
228
            midPoint += step;
 
229
        }
 
230
 
 
231
        return halfStep * sum;
 
232
 
 
233
    }
 
234
 
 
235
}