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