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.math3.analysis.integration;
018    
019    import org.apache.commons.math3.analysis.UnivariateFunction;
020    import org.apache.commons.math3.analysis.integration.gauss.GaussIntegratorFactory;
021    import org.apache.commons.math3.analysis.integration.gauss.GaussIntegrator;
022    import org.apache.commons.math3.exception.MaxCountExceededException;
023    import org.apache.commons.math3.exception.NotStrictlyPositiveException;
024    import org.apache.commons.math3.exception.NumberIsTooSmallException;
025    import org.apache.commons.math3.exception.TooManyEvaluationsException;
026    import org.apache.commons.math3.util.FastMath;
027    
028    /**
029     * This algorithm divides the integration interval into equally-sized
030     * sub-interval and on each of them performs a
031     * <a href="http://mathworld.wolfram.com/Legendre-GaussQuadrature.html">
032     * Legendre-Gauss</a> quadrature.
033     *
034     * @version $Id: IterativeLegendreGaussIntegrator.java 1416643 2012-12-03 19:37:14Z tn $
035     * @since 3.1
036     */
037    
038    public class IterativeLegendreGaussIntegrator
039        extends BaseAbstractUnivariateIntegrator {
040        /** Factory that computes the points and weights. */
041        private static final GaussIntegratorFactory FACTORY
042            = new GaussIntegratorFactory();
043        /** Number of integration points (per interval). */
044        private final int numberOfPoints;
045    
046        /**
047         * Builds an integrator with given accuracies and iterations counts.
048         *
049         * @param n Number of integration points.
050         * @param relativeAccuracy Relative accuracy of the result.
051         * @param absoluteAccuracy Absolute accuracy of the result.
052         * @param minimalIterationCount Minimum number of iterations.
053         * @param maximalIterationCount Maximum number of iterations.
054         * @throws NotStrictlyPositiveException if minimal number of iterations
055         * is not strictly positive.
056         * @throws NumberIsTooSmallException if maximal number of iterations
057         * is smaller than or equal to the minimal number of iterations.
058         */
059        public IterativeLegendreGaussIntegrator(final int n,
060                                                final double relativeAccuracy,
061                                                final double absoluteAccuracy,
062                                                final int minimalIterationCount,
063                                                final int maximalIterationCount)
064            throws NotStrictlyPositiveException, NumberIsTooSmallException {
065            super(relativeAccuracy, absoluteAccuracy, minimalIterationCount, maximalIterationCount);
066            numberOfPoints = n;
067        }
068    
069        /**
070         * Builds an integrator with given accuracies.
071         *
072         * @param n Number of integration points.
073         * @param relativeAccuracy Relative accuracy of the result.
074         * @param absoluteAccuracy Absolute accuracy of the result.
075         */
076        public IterativeLegendreGaussIntegrator(final int n,
077                                                final double relativeAccuracy,
078                                                final double absoluteAccuracy) {
079            this(n, relativeAccuracy, absoluteAccuracy,
080                 DEFAULT_MIN_ITERATIONS_COUNT, DEFAULT_MAX_ITERATIONS_COUNT);
081        }
082    
083        /**
084         * Builds an integrator with given iteration counts.
085         *
086         * @param n Number of integration points.
087         * @param minimalIterationCount Minimum number of iterations.
088         * @param maximalIterationCount Maximum number of iterations.
089         * @throws NotStrictlyPositiveException if minimal number of iterations
090         * is not strictly positive.
091         * @throws NumberIsTooSmallException if maximal number of iterations
092         * is smaller than or equal to the minimal number of iterations.
093         */
094        public IterativeLegendreGaussIntegrator(final int n,
095                                                final int minimalIterationCount,
096                                                final int maximalIterationCount) {
097            this(n, DEFAULT_RELATIVE_ACCURACY, DEFAULT_ABSOLUTE_ACCURACY,
098                 minimalIterationCount, maximalIterationCount);
099        }
100    
101        /** {@inheritDoc} */
102        @Override
103        protected double doIntegrate()
104            throws TooManyEvaluationsException, MaxCountExceededException {
105            // Compute first estimate with a single step.
106            double oldt = stage(1);
107    
108            int n = 2;
109            while (true) {
110                // Improve integral with a larger number of steps.
111                final double t = stage(n);
112    
113                // Estimate the error.
114                final double delta = FastMath.abs(t - oldt);
115                final double limit =
116                    FastMath.max(getAbsoluteAccuracy(),
117                                 getRelativeAccuracy() * (FastMath.abs(oldt) + FastMath.abs(t)) * 0.5);
118    
119                // check convergence
120                if (iterations.getCount() + 1 >= getMinimalIterationCount() &&
121                    delta <= limit) {
122                    return t;
123                }
124    
125                // Prepare next iteration.
126                final double ratio = FastMath.min(4, FastMath.pow(delta / limit, 0.5 / numberOfPoints));
127                n = FastMath.max((int) (ratio * n), n + 1);
128                oldt = t;
129                iterations.incrementCount();
130            }
131        }
132    
133        /**
134         * Compute the n-th stage integral.
135         *
136         * @param n Number of steps.
137         * @return the value of n-th stage integral.
138         * @throws TooManyEvaluationsException if the maximum number of evaluations
139         * is exceeded.
140         */
141        private double stage(final int n)
142            throws TooManyEvaluationsException {
143            // Function to be integrated is stored in the base class.
144            final UnivariateFunction f = new UnivariateFunction() {
145                    public double value(double x) {
146                        return computeObjectiveValue(x);
147                    }
148                };
149    
150            final double min = getMin();
151            final double max = getMax();
152            final double step = (max - min) / n;
153    
154            double sum = 0;
155            for (int i = 0; i < n; i++) {
156                // Integrate over each sub-interval [a, b].
157                final double a = min + i * step;
158                final double b = a + step;
159                final GaussIntegrator g = FACTORY.legendreHighPrecision(numberOfPoints, a, b);
160                sum += g.integrate(f);
161            }
162    
163            return sum;
164        }
165    }