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   * @since 1.2
51   * @deprecated As of 3.1 (to be removed in 4.0). Please use
52   * {@link IterativeLegendreGaussIntegrator} instead.
53   */
54  @Deprecated
55  public class LegendreGaussIntegrator extends BaseAbstractUnivariateIntegrator {
56  
57      /** Abscissas for the 2 points method. */
58      private static final double[] ABSCISSAS_2 = {
59          -1.0 / FastMath.sqrt(3.0),
60           1.0 / FastMath.sqrt(3.0)
61      };
62  
63      /** Weights for the 2 points method. */
64      private static final double[] WEIGHTS_2 = {
65          1.0,
66          1.0
67      };
68  
69      /** Abscissas for the 3 points method. */
70      private static final double[] ABSCISSAS_3 = {
71          -FastMath.sqrt(0.6),
72           0.0,
73           FastMath.sqrt(0.6)
74      };
75  
76      /** Weights for the 3 points method. */
77      private static final double[] WEIGHTS_3 = {
78          5.0 / 9.0,
79          8.0 / 9.0,
80          5.0 / 9.0
81      };
82  
83      /** Abscissas for the 4 points method. */
84      private static final double[] ABSCISSAS_4 = {
85          -FastMath.sqrt((15.0 + 2.0 * FastMath.sqrt(30.0)) / 35.0),
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      };
90  
91      /** Weights for the 4 points method. */
92      private static final double[] WEIGHTS_4 = {
93          (90.0 - 5.0 * FastMath.sqrt(30.0)) / 180.0,
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      };
98  
99      /** 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 }