View Javadoc
1   /*
2    * Licensed to the Apache Software Foundation (ASF) under one or more
3    * contributor license agreements.  See the NOTICE file distributed with
4    * this work for additional information regarding copyright ownership.
5    * The ASF licenses this file to You under the Apache License, Version 2.0
6    * (the "License"); you may not use this file except in compliance with
7    * the License.  You may obtain a copy of the License at
8    *
9    *      http://www.apache.org/licenses/LICENSE-2.0
10   *
11   * Unless required by applicable law or agreed to in writing, software
12   * distributed under the License is distributed on an "AS IS" BASIS,
13   * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
14   * See the License for the specific language governing permissions and
15   * limitations under the License.
16   */
17  package org.apache.commons.math3.analysis.integration;
18  
19  import org.apache.commons.math3.exception.MathIllegalArgumentException;
20  import org.apache.commons.math3.exception.MaxCountExceededException;
21  import org.apache.commons.math3.exception.NotStrictlyPositiveException;
22  import org.apache.commons.math3.exception.NumberIsTooSmallException;
23  import org.apache.commons.math3.exception.TooManyEvaluationsException;
24  import org.apache.commons.math3.exception.util.LocalizedFormats;
25  import org.apache.commons.math3.util.FastMath;
26  
27  /**
28   * Implements the <a href="http://mathworld.wolfram.com/Legendre-GaussQuadrature.html">
29   * Legendre-Gauss</a> quadrature formula.
30   * <p>
31   * Legendre-Gauss integrators are efficient integrators that can
32   * accurately integrate functions with few function evaluations. A
33   * Legendre-Gauss integrator using an n-points quadrature formula can
34   * integrate 2n-1 degree polynomials exactly.
35   * </p>
36   * <p>
37   * These integrators evaluate the function on n carefully chosen
38   * abscissas in each step interval (mapped to the canonical [-1,1] interval).
39   * The evaluation abscissas are not evenly spaced and none of them are
40   * at the interval endpoints. This implies the function integrated can be
41   * undefined at integration interval endpoints.
42   * </p>
43   * <p>
44   * The evaluation abscissas x<sub>i</sub> are the roots of the degree n
45   * Legendre polynomial. The weights a<sub>i</sub> of the quadrature formula
46   * integrals from -1 to +1 &int; Li<sup>2</sup> where Li (x) =
47   * &prod; (x-x<sub>k</sub>)/(x<sub>i</sub>-x<sub>k</sub>) for k != i.
48   * </p>
49   * <p>
50   * @version $Id: LegendreGaussIntegrator.java 1455194 2013-03-11 15:45:54Z luc $
51   * @since 1.2
52   * @deprecated As of 3.1 (to be removed in 4.0). Please use
53   * {@link IterativeLegendreGaussIntegrator} instead.
54   */
55  @Deprecated
56  public class LegendreGaussIntegrator extends BaseAbstractUnivariateIntegrator {
57  
58      /** Abscissas for the 2 points method. */
59      private static final double[] ABSCISSAS_2 = {
60          -1.0 / FastMath.sqrt(3.0),
61           1.0 / FastMath.sqrt(3.0)
62      };
63  
64      /** Weights for the 2 points method. */
65      private static final double[] WEIGHTS_2 = {
66          1.0,
67          1.0
68      };
69  
70      /** Abscissas for the 3 points method. */
71      private static final double[] ABSCISSAS_3 = {
72          -FastMath.sqrt(0.6),
73           0.0,
74           FastMath.sqrt(0.6)
75      };
76  
77      /** Weights for the 3 points method. */
78      private static final double[] WEIGHTS_3 = {
79          5.0 / 9.0,
80          8.0 / 9.0,
81          5.0 / 9.0
82      };
83  
84      /** Abscissas for the 4 points method. */
85      private static final double[] ABSCISSAS_4 = {
86          -FastMath.sqrt((15.0 + 2.0 * FastMath.sqrt(30.0)) / 35.0),
87          -FastMath.sqrt((15.0 - 2.0 * FastMath.sqrt(30.0)) / 35.0),
88           FastMath.sqrt((15.0 - 2.0 * FastMath.sqrt(30.0)) / 35.0),
89           FastMath.sqrt((15.0 + 2.0 * FastMath.sqrt(30.0)) / 35.0)
90      };
91  
92      /** Weights for the 4 points method. */
93      private static final double[] WEIGHTS_4 = {
94          (90.0 - 5.0 * FastMath.sqrt(30.0)) / 180.0,
95          (90.0 + 5.0 * FastMath.sqrt(30.0)) / 180.0,
96          (90.0 + 5.0 * FastMath.sqrt(30.0)) / 180.0,
97          (90.0 - 5.0 * FastMath.sqrt(30.0)) / 180.0
98      };
99  
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 MathIllegalArgumentException if number of points is out of [2; 5]
132      * @exception NotStrictlyPositiveException if minimal number of iterations
133      * is not strictly positive
134      * @exception NumberIsTooSmallException if maximal number of iterations
135      * is lesser than or equal to the minimal number of iterations
136      */
137     public LegendreGaussIntegrator(final int n,
138                                    final double relativeAccuracy,
139                                    final double absoluteAccuracy,
140                                    final int minimalIterationCount,
141                                    final int maximalIterationCount)
142         throws MathIllegalArgumentException, NotStrictlyPositiveException, NumberIsTooSmallException {
143         super(relativeAccuracy, absoluteAccuracy, minimalIterationCount, maximalIterationCount);
144         switch(n) {
145         case 2 :
146             abscissas = ABSCISSAS_2;
147             weights   = WEIGHTS_2;
148             break;
149         case 3 :
150             abscissas = ABSCISSAS_3;
151             weights   = WEIGHTS_3;
152             break;
153         case 4 :
154             abscissas = ABSCISSAS_4;
155             weights   = WEIGHTS_4;
156             break;
157         case 5 :
158             abscissas = ABSCISSAS_5;
159             weights   = WEIGHTS_5;
160             break;
161         default :
162             throw new MathIllegalArgumentException(
163                     LocalizedFormats.N_POINTS_GAUSS_LEGENDRE_INTEGRATOR_NOT_SUPPORTED,
164                     n, 2, 5);
165         }
166 
167     }
168 
169     /**
170      * Build a Legendre-Gauss integrator with given accuracies.
171      * @param n number of points desired (must be between 2 and 5 inclusive)
172      * @param relativeAccuracy relative accuracy of the result
173      * @param absoluteAccuracy absolute accuracy of the result
174      * @exception MathIllegalArgumentException if number of points is out of [2; 5]
175      */
176     public LegendreGaussIntegrator(final int n,
177                                    final double relativeAccuracy,
178                                    final double absoluteAccuracy)
179         throws MathIllegalArgumentException {
180         this(n, relativeAccuracy, absoluteAccuracy,
181              DEFAULT_MIN_ITERATIONS_COUNT, DEFAULT_MAX_ITERATIONS_COUNT);
182     }
183 
184     /**
185      * Build a Legendre-Gauss integrator with given iteration counts.
186      * @param n number of points desired (must be between 2 and 5 inclusive)
187      * @param minimalIterationCount minimum number of iterations
188      * @param maximalIterationCount maximum number of iterations
189      * @exception MathIllegalArgumentException if number of points is out of [2; 5]
190      * @exception NotStrictlyPositiveException if minimal number of iterations
191      * is not strictly positive
192      * @exception NumberIsTooSmallException if maximal number of iterations
193      * is lesser than or equal to the minimal number of iterations
194      */
195     public LegendreGaussIntegrator(final int n,
196                                    final int minimalIterationCount,
197                                    final int maximalIterationCount)
198         throws MathIllegalArgumentException {
199         this(n, DEFAULT_RELATIVE_ACCURACY, DEFAULT_ABSOLUTE_ACCURACY,
200              minimalIterationCount, maximalIterationCount);
201     }
202 
203     /** {@inheritDoc} */
204     @Override
205     protected double doIntegrate()
206         throws MathIllegalArgumentException, TooManyEvaluationsException, MaxCountExceededException {
207 
208         // compute first estimate with a single step
209         double oldt = stage(1);
210 
211         int n = 2;
212         while (true) {
213 
214             // improve integral with a larger number of steps
215             final double t = stage(n);
216 
217             // estimate error
218             final double delta = FastMath.abs(t - oldt);
219             final double limit =
220                 FastMath.max(getAbsoluteAccuracy(),
221                              getRelativeAccuracy() * (FastMath.abs(oldt) + FastMath.abs(t)) * 0.5);
222 
223             // check convergence
224             if ((iterations.getCount() + 1 >= getMinimalIterationCount()) && (delta <= limit)) {
225                 return t;
226             }
227 
228             // prepare next iteration
229             double ratio = FastMath.min(4, FastMath.pow(delta / limit, 0.5 / abscissas.length));
230             n = FastMath.max((int) (ratio * n), n + 1);
231             oldt = t;
232             iterations.incrementCount();
233 
234         }
235 
236     }
237 
238     /**
239      * Compute the n-th stage integral.
240      * @param n number of steps
241      * @return the value of n-th stage integral
242      * @throws TooManyEvaluationsException if the maximum number of evaluations
243      * is exceeded.
244      */
245     private double stage(final int n)
246         throws TooManyEvaluationsException {
247 
248         // set up the step for the current stage
249         final double step     = (getMax() - getMin()) / n;
250         final double halfStep = step / 2.0;
251 
252         // integrate over all elementary steps
253         double midPoint = getMin() + halfStep;
254         double sum = 0.0;
255         for (int i = 0; i < n; ++i) {
256             for (int j = 0; j < abscissas.length; ++j) {
257                 sum += weights[j] * computeObjectiveValue(midPoint + halfStep * abscissas[j]);
258             }
259             midPoint += step;
260         }
261 
262         return halfStep * sum;
263 
264     }
265 
266 }