001    /*
002     * Licensed to the Apache Software Foundation (ASF) under one or more
003     * contributor license agreements.  See the NOTICE file distributed with
004     * this work for additional information regarding copyright ownership.
005     * The ASF licenses this file to You under the Apache License, Version 2.0
006     * (the "License"); you may not use this file except in compliance with
007     * the License.  You may obtain a copy of the License at
008     *
009     *      http://www.apache.org/licenses/LICENSE-2.0
010     *
011     * Unless required by applicable law or agreed to in writing, software
012     * distributed under the License is distributed on an "AS IS" BASIS,
013     * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
014     * See the License for the specific language governing permissions and
015     * limitations under the License.
016     */
017    package org.apache.commons.math.analysis.integration;
018    
019    import org.apache.commons.math.ConvergenceException;
020    import org.apache.commons.math.FunctionEvaluationException;
021    import org.apache.commons.math.MathRuntimeException;
022    import org.apache.commons.math.MaxIterationsExceededException;
023    import org.apache.commons.math.analysis.UnivariateRealFunction;
024    
025    /**
026     * Implements the <a href="http://mathworld.wolfram.com/Legendre-GaussQuadrature.html">
027     * Legendre-Gauss</a> quadrature formula.
028     * <p>
029     * Legendre-Gauss integrators are efficient integrators that can
030     * accurately integrate functions with few functions evaluations. A
031     * Legendre-Gauss integrator using an n-points quadrature formula can
032     * integrate exactly 2n-1 degree polynomialss.
033     * </p>
034     * <p>
035     * These integrators evaluate the function on n carefully chosen
036     * abscissas in each step interval (mapped to the canonical [-1  1] interval).
037     * The evaluation abscissas are not evenly spaced and none of them are
038     * at the interval endpoints. This implies the function integrated can be
039     * undefined at integration interval endpoints.
040     * </p>
041     * <p>
042     * The evaluation abscissas x<sub>i</sub> are the roots of the degree n
043     * Legendre polynomial. The weights a<sub>i</sub> of the quadrature formula
044     * integrals from -1 to +1 &int; Li<sup>2</sup> where Li (x) =
045     * &prod; (x-x<sub>k</sub>)/(x<sub>i</sub>-x<sub>k</sub>) for k != i.
046     * </p>
047     * <p>
048     * @version $Revision: 799857 $ $Date: 2009-08-01 09:07:12 -0400 (Sat, 01 Aug 2009) $
049     * @since 1.2
050     */
051    
052    public class LegendreGaussIntegrator extends UnivariateRealIntegratorImpl {
053    
054        /** Abscissas for the 2 points method. */
055        private static final double[] ABSCISSAS_2 = {
056            -1.0 / Math.sqrt(3.0),
057             1.0 / Math.sqrt(3.0)
058        };
059    
060        /** Weights for the 2 points method. */
061        private static final double[] WEIGHTS_2 = {
062            1.0,
063            1.0
064        };
065    
066        /** Abscissas for the 3 points method. */
067        private static final double[] ABSCISSAS_3 = {
068            -Math.sqrt(0.6),
069             0.0,
070             Math.sqrt(0.6)
071        };
072    
073        /** Weights for the 3 points method. */
074        private static final double[] WEIGHTS_3 = {
075            5.0 / 9.0,
076            8.0 / 9.0,
077            5.0 / 9.0
078        };
079    
080        /** Abscissas for the 4 points method. */
081        private static final double[] ABSCISSAS_4 = {
082            -Math.sqrt((15.0 + 2.0 * Math.sqrt(30.0)) / 35.0),
083            -Math.sqrt((15.0 - 2.0 * Math.sqrt(30.0)) / 35.0),
084             Math.sqrt((15.0 - 2.0 * Math.sqrt(30.0)) / 35.0),
085             Math.sqrt((15.0 + 2.0 * Math.sqrt(30.0)) / 35.0)
086        };
087    
088        /** Weights for the 4 points method. */
089        private static final double[] WEIGHTS_4 = {
090            (90.0 - 5.0 * Math.sqrt(30.0)) / 180.0,
091            (90.0 + 5.0 * Math.sqrt(30.0)) / 180.0,
092            (90.0 + 5.0 * Math.sqrt(30.0)) / 180.0,
093            (90.0 - 5.0 * Math.sqrt(30.0)) / 180.0
094        };
095    
096        /** Abscissas for the 5 points method. */
097        private static final double[] ABSCISSAS_5 = {
098            -Math.sqrt((35.0 + 2.0 * Math.sqrt(70.0)) / 63.0),
099            -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    }