    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    import org.apache.commons.math.exception.util.LocalizedFormats;
    025    import org.apache.commons.math.util.FastMath;
    026
    027    /**
    028     * Implements the <a href="http://mathworld.wolfram.com/Legendre-GaussQuadrature.html">
    029     * Legendre-Gauss</a> quadrature formula.
    030     * <p>
    031     * Legendre-Gauss integrators are efficient integrators that can
    032     * accurately integrate functions with few functions evaluations. A
    033     * Legendre-Gauss integrator using an n-points quadrature formula can
    034     * integrate exactly 2n-1 degree polynomials.
    035     * </p>
    036     * <p>
    037     * These integrators evaluate the function on n carefully chosen
    038     * abscissas in each step interval (mapped to the canonical [-1  1] interval).
    039     * The evaluation abscissas are not evenly spaced and none of them are
    040     * at the interval endpoints. This implies the function integrated can be
    041     * undefined at integration interval endpoints.
    042     * </p>
    043     * <p>
    044     * The evaluation abscissas x<sub>i</sub> are the roots of the degree n
    045     * Legendre polynomial. The weights a<sub>i</sub> of the quadrature formula
    046     * integrals from -1 to +1 &int; Li<sup>2</sup> where Li (x) =
    047     * &prod; (x-x<sub>k</sub>)/(x<sub>i</sub>-x<sub>k</sub>) for k != i.
    048     * </p>
    049     * <p>
    050     * @version $Revision: 1070725 $ $Date: 2011-02-15 02:31:12 +0100 (mar. 15 f??vr. 2011) $
    051     * @since 1.2
    052     */
    053
    054    public class LegendreGaussIntegrator extends UnivariateRealIntegratorImpl {
    055
    056        /** Abscissas for the 2 points method. */
    057        private static final double[] ABSCISSAS_2 = {
    058            -1.0 / FastMath.sqrt(3.0),
    059             1.0 / FastMath.sqrt(3.0)
    060        };
    061
    062        /** Weights for the 2 points method. */
    063        private static final double[] WEIGHTS_2 = {
    064            1.0,
    065            1.0
    066        };
    067
    068        /** Abscissas for the 3 points method. */
    069        private static final double[] ABSCISSAS_3 = {
    070            -FastMath.sqrt(0.6),
    071             0.0,
    072             FastMath.sqrt(0.6)
    073        };
    074
    075        /** Weights for the 3 points method. */
    076        private static final double[] WEIGHTS_3 = {
    077            5.0 / 9.0,
    078            8.0 / 9.0,
    079            5.0 / 9.0
    080        };
    081
    082        /** Abscissas for the 4 points method. */
    083        private static final double[] ABSCISSAS_4 = {
    084            -FastMath.sqrt((15.0 + 2.0 * FastMath.sqrt(30.0)) / 35.0),
    085            -FastMath.sqrt((15.0 - 2.0 * FastMath.sqrt(30.0)) / 35.0),
    086             FastMath.sqrt((15.0 - 2.0 * FastMath.sqrt(30.0)) / 35.0),
    087             FastMath.sqrt((15.0 + 2.0 * FastMath.sqrt(30.0)) / 35.0)
    088        };
    089
    090        /** Weights for the 4 points method. */
    091        private static final double[] WEIGHTS_4 = {
    092            (90.0 - 5.0 * FastMath.sqrt(30.0)) / 180.0,
    093            (90.0 + 5.0 * FastMath.sqrt(30.0)) / 180.0,
    094            (90.0 + 5.0 * FastMath.sqrt(30.0)) / 180.0,
    095            (90.0 - 5.0 * FastMath.sqrt(30.0)) / 180.0
    096        };
    097
    098        /** Abscissas for the 5 points method. */
    099        private static final double[] ABSCISSAS_5 = {
    100            -FastMath.sqrt((35.0 + 2.0 * FastMath.sqrt(70.0)) / 63.0),
    101            -FastMath.sqrt((35.0 - 2.0 * FastMath.sqrt(70.0)) / 63.0),
    102             0.0,
    103             FastMath.sqrt((35.0 - 2.0 * FastMath.sqrt(70.0)) / 63.0),
    104             FastMath.sqrt((35.0 + 2.0 * FastMath.sqrt(70.0)) / 63.0)
    105        };
    106
    107        /** Weights for the 5 points method. */
    108        private static final double[] WEIGHTS_5 = {
    109            (322.0 - 13.0 * FastMath.sqrt(70.0)) / 900.0,
    110            (322.0 + 13.0 * FastMath.sqrt(70.0)) / 900.0,
    111            128.0 / 225.0,
    112            (322.0 + 13.0 * FastMath.sqrt(70.0)) / 900.0,
    113            (322.0 - 13.0 * FastMath.sqrt(70.0)) / 900.0
    114        };
    115
    116        /** Abscissas for the current method. */
    117        private final double[] abscissas;
    118
    119        /** Weights for the current method. */
    120        private final double[] weights;
    121
    122        /**
    123         * Build a Legendre-Gauss integrator.
    124         * @param n number of points desired (must be between 2 and 5 inclusive)
    125         * @param defaultMaximalIterationCount maximum number of iterations
    126         * @exception IllegalArgumentException if the number of points is not
    127         * in the supported range
    128         */
    129        public LegendreGaussIntegrator(final int n, final int defaultMaximalIterationCount)
    130            throws IllegalArgumentException {
    131            super(defaultMaximalIterationCount);
    132            switch(n) {
    133            case 2 :
    134                abscissas = ABSCISSAS_2;
    135                weights   = WEIGHTS_2;
    136                break;
    137            case 3 :
    138                abscissas = ABSCISSAS_3;
    139                weights   = WEIGHTS_3;
    140                break;
    141            case 4 :
    142                abscissas = ABSCISSAS_4;
    143                weights   = WEIGHTS_4;
    144                break;
    145            case 5 :
    146                abscissas = ABSCISSAS_5;
    147                weights   = WEIGHTS_5;
    148                break;
    149            default :
    150                throw MathRuntimeException.createIllegalArgumentException(
    151                        LocalizedFormats.N_POINTS_GAUSS_LEGENDRE_INTEGRATOR_NOT_SUPPORTED,
    152                        n, 2, 5);
    153            }
    154
    155        }
    156
    157        /** {@inheritDoc} */
    158        @Deprecated
    159        public double integrate(final double min, final double max)
    160            throws ConvergenceException,  FunctionEvaluationException, IllegalArgumentException {
    161            return integrate(f, min, max);
    162        }
    163
    164        /** {@inheritDoc} */
    165        public double integrate(final UnivariateRealFunction f, final double min, final double max)
    166            throws ConvergenceException,  FunctionEvaluationException, IllegalArgumentException {
    167
    168            clearResult();
    169            verifyInterval(min, max);
    170            verifyIterationCount();
    171
    172            // compute first estimate with a single step
    173            double oldt = stage(f, min, max, 1);
    174
    175            int n = 2;
    176            for (int i = 0; i < maximalIterationCount; ++i) {
    177
    178                // improve integral with a larger number of steps
    179                final double t = stage(f, min, max, n);
    180
    181                // estimate error
    182                final double delta = FastMath.abs(t - oldt);
    183                final double limit =
    184                    FastMath.max(absoluteAccuracy,
    185                             relativeAccuracy * (FastMath.abs(oldt) + FastMath.abs(t)) * 0.5);
    186
    187                // check convergence
    188                if ((i + 1 >= minimalIterationCount) && (delta <= limit)) {
    189                    setResult(t, i);
    190                    return result;
    191                }
    192
    193                // prepare next iteration
    194                double ratio = FastMath.min(4, FastMath.pow(delta / limit, 0.5 / abscissas.length));
    195                n = FastMath.max((int) (ratio * n), n + 1);
    196                oldt = t;
    197
    198            }
    199
    200            throw new MaxIterationsExceededException(maximalIterationCount);
    201
    202        }
    203
    204        /**
    205         * Compute the n-th stage integral.
    206         * @param f the integrand function
    207         * @param min the lower bound for the interval
    208         * @param max the upper bound for the interval
    209         * @param n number of steps
    210         * @return the value of n-th stage integral
    211         * @throws FunctionEvaluationException if an error occurs evaluating the
    212         * function
    213         */
    214        private double stage(final UnivariateRealFunction f,
    215                             final double min, final double max, final int n)
    216            throws FunctionEvaluationException {
    217
    218            // set up the step for the current stage
    219            final double step     = (max - min) / n;
    220            final double halfStep = step / 2.0;
    221
    222            // integrate over all elementary steps
    223            double midPoint = min + halfStep;
    224            double sum = 0.0;
    225            for (int i = 0; i < n; ++i) {
    226                for (int j = 0; j < abscissas.length; ++j) {
    227                    sum += weights[j] * f.value(midPoint + halfStep * abscissas[j]);
    228                }
    229                midPoint += step;
    230            }
    231
    232            return halfStep * sum;
    233
    234        }
    235
    236    }