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.exception.MathIllegalArgumentException;
020    import org.apache.commons.math3.exception.MaxCountExceededException;
021    import org.apache.commons.math3.exception.NotStrictlyPositiveException;
022    import org.apache.commons.math3.exception.NumberIsTooSmallException;
023    import org.apache.commons.math3.exception.TooManyEvaluationsException;
024    import org.apache.commons.math3.exception.util.LocalizedFormats;
025    import org.apache.commons.math3.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 function evaluations. A
033     * Legendre-Gauss integrator using an n-points quadrature formula can
034     * integrate 2n-1 degree polynomials exactly.
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 $Id: LegendreGaussIntegrator.java 1364452 2012-07-22 22:30:01Z erans $
051     * @since 1.2
052     * @deprecated As of 3.1 (to be removed in 4.0). Please use
053     * {@link IterativeLegendreGaussIntegrator} instead.
054     */
055    @Deprecated
056    public class LegendreGaussIntegrator extends BaseAbstractUnivariateIntegrator {
057    
058        /** Abscissas for the 2 points method. */
059        private static final double[] ABSCISSAS_2 = {
060            -1.0 / FastMath.sqrt(3.0),
061             1.0 / FastMath.sqrt(3.0)
062        };
063    
064        /** Weights for the 2 points method. */
065        private static final double[] WEIGHTS_2 = {
066            1.0,
067            1.0
068        };
069    
070        /** Abscissas for the 3 points method. */
071        private static final double[] ABSCISSAS_3 = {
072            -FastMath.sqrt(0.6),
073             0.0,
074             FastMath.sqrt(0.6)
075        };
076    
077        /** Weights for the 3 points method. */
078        private static final double[] WEIGHTS_3 = {
079            5.0 / 9.0,
080            8.0 / 9.0,
081            5.0 / 9.0
082        };
083    
084        /** Abscissas for the 4 points method. */
085        private static final double[] ABSCISSAS_4 = {
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             FastMath.sqrt((15.0 - 2.0 * FastMath.sqrt(30.0)) / 35.0),
089             FastMath.sqrt((15.0 + 2.0 * FastMath.sqrt(30.0)) / 35.0)
090        };
091    
092        /** Weights for the 4 points method. */
093        private static final double[] WEIGHTS_4 = {
094            (90.0 - 5.0 * FastMath.sqrt(30.0)) / 180.0,
095            (90.0 + 5.0 * FastMath.sqrt(30.0)) / 180.0,
096            (90.0 + 5.0 * FastMath.sqrt(30.0)) / 180.0,
097            (90.0 - 5.0 * FastMath.sqrt(30.0)) / 180.0
098        };
099    
100        /** Abscissas for the 5 points method. */
101        private static final double[] ABSCISSAS_5 = {
102            -FastMath.sqrt((35.0 + 2.0 * FastMath.sqrt(70.0)) / 63.0),
103            -FastMath.sqrt((35.0 - 2.0 * FastMath.sqrt(70.0)) / 63.0),
104             0.0,
105             FastMath.sqrt((35.0 - 2.0 * FastMath.sqrt(70.0)) / 63.0),
106             FastMath.sqrt((35.0 + 2.0 * FastMath.sqrt(70.0)) / 63.0)
107        };
108    
109        /** Weights for the 5 points method. */
110        private static final double[] WEIGHTS_5 = {
111            (322.0 - 13.0 * FastMath.sqrt(70.0)) / 900.0,
112            (322.0 + 13.0 * FastMath.sqrt(70.0)) / 900.0,
113            128.0 / 225.0,
114            (322.0 + 13.0 * FastMath.sqrt(70.0)) / 900.0,
115            (322.0 - 13.0 * FastMath.sqrt(70.0)) / 900.0
116        };
117    
118        /** Abscissas for the current method. */
119        private final double[] abscissas;
120    
121        /** Weights for the current method. */
122        private final double[] weights;
123    
124        /**
125         * Build a Legendre-Gauss integrator with given accuracies and iterations counts.
126         * @param n number of points desired (must be between 2 and 5 inclusive)
127         * @param relativeAccuracy relative accuracy of the result
128         * @param absoluteAccuracy absolute accuracy of the result
129         * @param minimalIterationCount minimum number of iterations
130         * @param maximalIterationCount maximum number of iterations
131         * @exception NotStrictlyPositiveException if minimal number of iterations
132         * is not strictly positive
133         * @exception NumberIsTooSmallException if maximal number of iterations
134         * is lesser than or equal to the minimal number of iterations
135         */
136        public LegendreGaussIntegrator(final int n,
137                                       final double relativeAccuracy,
138                                       final double absoluteAccuracy,
139                                       final int minimalIterationCount,
140                                       final int maximalIterationCount)
141            throws NotStrictlyPositiveException, NumberIsTooSmallException {
142            super(relativeAccuracy, absoluteAccuracy, minimalIterationCount, maximalIterationCount);
143            switch(n) {
144            case 2 :
145                abscissas = ABSCISSAS_2;
146                weights   = WEIGHTS_2;
147                break;
148            case 3 :
149                abscissas = ABSCISSAS_3;
150                weights   = WEIGHTS_3;
151                break;
152            case 4 :
153                abscissas = ABSCISSAS_4;
154                weights   = WEIGHTS_4;
155                break;
156            case 5 :
157                abscissas = ABSCISSAS_5;
158                weights   = WEIGHTS_5;
159                break;
160            default :
161                throw new MathIllegalArgumentException(
162                        LocalizedFormats.N_POINTS_GAUSS_LEGENDRE_INTEGRATOR_NOT_SUPPORTED,
163                        n, 2, 5);
164            }
165    
166        }
167    
168        /**
169         * Build a Legendre-Gauss integrator with given accuracies.
170         * @param n number of points desired (must be between 2 and 5 inclusive)
171         * @param relativeAccuracy relative accuracy of the result
172         * @param absoluteAccuracy absolute accuracy of the result
173         */
174        public LegendreGaussIntegrator(final int n,
175                                       final double relativeAccuracy,
176                                       final double absoluteAccuracy) {
177            this(n, relativeAccuracy, absoluteAccuracy,
178                 DEFAULT_MIN_ITERATIONS_COUNT, DEFAULT_MAX_ITERATIONS_COUNT);
179        }
180    
181        /**
182         * Build a Legendre-Gauss integrator with given iteration counts.
183         * @param n number of points desired (must be between 2 and 5 inclusive)
184         * @param minimalIterationCount minimum number of iterations
185         * @param maximalIterationCount maximum number of iterations
186         * @exception NotStrictlyPositiveException if minimal number of iterations
187         * is not strictly positive
188         * @exception NumberIsTooSmallException if maximal number of iterations
189         * is lesser than or equal to the minimal number of iterations
190         */
191        public LegendreGaussIntegrator(final int n,
192                                       final int minimalIterationCount,
193                                       final int maximalIterationCount) {
194            this(n, DEFAULT_RELATIVE_ACCURACY, DEFAULT_ABSOLUTE_ACCURACY,
195                 minimalIterationCount, maximalIterationCount);
196        }
197    
198        /** {@inheritDoc} */
199        @Override
200        protected double doIntegrate()
201            throws TooManyEvaluationsException, MaxCountExceededException {
202    
203            // compute first estimate with a single step
204            double oldt = stage(1);
205    
206            int n = 2;
207            while (true) {
208    
209                // improve integral with a larger number of steps
210                final double t = stage(n);
211    
212                // estimate error
213                final double delta = FastMath.abs(t - oldt);
214                final double limit =
215                    FastMath.max(getAbsoluteAccuracy(),
216                                 getRelativeAccuracy() * (FastMath.abs(oldt) + FastMath.abs(t)) * 0.5);
217    
218                // check convergence
219                if ((iterations.getCount() + 1 >= getMinimalIterationCount()) && (delta <= limit)) {
220                    return t;
221                }
222    
223                // prepare next iteration
224                double ratio = FastMath.min(4, FastMath.pow(delta / limit, 0.5 / abscissas.length));
225                n = FastMath.max((int) (ratio * n), n + 1);
226                oldt = t;
227                iterations.incrementCount();
228    
229            }
230    
231        }
232    
233        /**
234         * Compute the n-th stage integral.
235         * @param n number of steps
236         * @return the value of n-th stage integral
237         * @throws TooManyEvaluationsException if the maximum number of evaluations
238         * is exceeded.
239         */
240        private double stage(final int n)
241            throws TooManyEvaluationsException {
242    
243            // set up the step for the current stage
244            final double step     = (getMax() - getMin()) / n;
245            final double halfStep = step / 2.0;
246    
247            // integrate over all elementary steps
248            double midPoint = getMin() + halfStep;
249            double sum = 0.0;
250            for (int i = 0; i < n; ++i) {
251                for (int j = 0; j < abscissas.length; ++j) {
252                    sum += weights[j] * computeObjectiveValue(midPoint + halfStep * abscissas[j]);
253                }
254                midPoint += step;
255            }
256    
257            return halfStep * sum;
258    
259        }
260    
261    }