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.math.analysis.integration;
18  
19  import org.apache.commons.math.ConvergenceException;
20  import org.apache.commons.math.FunctionEvaluationException;
21  import org.apache.commons.math.MathRuntimeException;
22  import org.apache.commons.math.MaxIterationsExceededException;
23  import org.apache.commons.math.analysis.UnivariateRealFunction;
24  
25  /**
26   * Implements the <a href="http://mathworld.wolfram.com/Legendre-GaussQuadrature.html">
27   * Legendre-Gauss</a> quadrature formula.
28   * <p>
29   * Legendre-Gauss integrators are efficient integrators that can
30   * accurately integrate functions with few functions evaluations. A
31   * Legendre-Gauss integrator using an n-points quadrature formula can
32   * integrate exactly 2n-1 degree polynomialss.
33   * </p>
34   * <p>
35   * These integrators evaluate the function on n carefully chosen
36   * abscissas in each step interval (mapped to the canonical [-1  1] interval).
37   * The evaluation abscissas are not evenly spaced and none of them are
38   * at the interval endpoints. This implies the function integrated can be
39   * undefined at integration interval endpoints.
40   * </p>
41   * <p>
42   * The evaluation abscissas x<sub>i</sub> are the roots of the degree n
43   * Legendre polynomial. The weights a<sub>i</sub> of the quadrature formula
44   * integrals from -1 to +1 &int; Li<sup>2</sup> where Li (x) =
45   * &prod; (x-x<sub>k</sub>)/(x<sub>i</sub>-x<sub>k</sub>) for k != i.
46   * </p>
47   * <p>
48   * @version $Revision: 799857 $ $Date: 2009-08-01 09:07:12 -0400 (Sat, 01 Aug 2009) $
49   * @since 1.2
50   */
51  
52  public class LegendreGaussIntegrator extends UnivariateRealIntegratorImpl {
53  
54      /** Abscissas for the 2 points method. */
55      private static final double[] ABSCISSAS_2 = {
56          -1.0 / Math.sqrt(3.0),
57           1.0 / Math.sqrt(3.0)
58      };
59  
60      /** Weights for the 2 points method. */
61      private static final double[] WEIGHTS_2 = {
62          1.0,
63          1.0
64      };
65  
66      /** Abscissas for the 3 points method. */
67      private static final double[] ABSCISSAS_3 = {
68          -Math.sqrt(0.6),
69           0.0,
70           Math.sqrt(0.6)
71      };
72  
73      /** Weights for the 3 points method. */
74      private static final double[] WEIGHTS_3 = {
75          5.0 / 9.0,
76          8.0 / 9.0,
77          5.0 / 9.0
78      };
79  
80      /** Abscissas for the 4 points method. */
81      private static final double[] ABSCISSAS_4 = {
82          -Math.sqrt((15.0 + 2.0 * Math.sqrt(30.0)) / 35.0),
83          -Math.sqrt((15.0 - 2.0 * Math.sqrt(30.0)) / 35.0),
84           Math.sqrt((15.0 - 2.0 * Math.sqrt(30.0)) / 35.0),
85           Math.sqrt((15.0 + 2.0 * Math.sqrt(30.0)) / 35.0)
86      };
87  
88      /** Weights for the 4 points method. */
89      private static final double[] WEIGHTS_4 = {
90          (90.0 - 5.0 * Math.sqrt(30.0)) / 180.0,
91          (90.0 + 5.0 * Math.sqrt(30.0)) / 180.0,
92          (90.0 + 5.0 * Math.sqrt(30.0)) / 180.0,
93          (90.0 - 5.0 * Math.sqrt(30.0)) / 180.0
94      };
95  
96      /** Abscissas for the 5 points method. */
97      private static final double[] ABSCISSAS_5 = {
98          -Math.sqrt((35.0 + 2.0 * Math.sqrt(70.0)) / 63.0),
99          -Math.sqrt((35.0 - 2.0 * Math.sqrt(70.0)) / 63.0),
100          0.0,
101          Math.sqrt((35.0 - 2.0 * Math.sqrt(70.0)) / 63.0),
102          Math.sqrt((35.0 + 2.0 * Math.sqrt(70.0)) / 63.0)
103     };
104 
105     /** Weights for the 5 points method. */
106     private static final double[] WEIGHTS_5 = {
107         (322.0 - 13.0 * Math.sqrt(70.0)) / 900.0,
108         (322.0 + 13.0 * Math.sqrt(70.0)) / 900.0,
109         128.0 / 225.0,
110         (322.0 + 13.0 * Math.sqrt(70.0)) / 900.0,
111         (322.0 - 13.0 * Math.sqrt(70.0)) / 900.0
112     };
113 
114     /** Abscissas for the current method. */
115     private final double[] abscissas;
116 
117     /** Weights for the current method. */
118     private final double[] weights;
119 
120     /** Build a Legendre-Gauss integrator.
121      * @param n number of points desired (must be between 2 and 5 inclusive)
122      * @param defaultMaximalIterationCount maximum number of iterations
123      * @exception IllegalArgumentException if the number of points is not
124      * in the supported range
125      */
126     public LegendreGaussIntegrator(final int n, final int defaultMaximalIterationCount)
127         throws IllegalArgumentException {
128         super(defaultMaximalIterationCount);
129         switch(n) {
130         case 2 :
131             abscissas = ABSCISSAS_2;
132             weights   = WEIGHTS_2;
133             break;
134         case 3 :
135             abscissas = ABSCISSAS_3;
136             weights   = WEIGHTS_3;
137             break;
138         case 4 :
139             abscissas = ABSCISSAS_4;
140             weights   = WEIGHTS_4;
141             break;
142         case 5 :
143             abscissas = ABSCISSAS_5;
144             weights   = WEIGHTS_5;
145             break;
146         default :
147             throw MathRuntimeException.createIllegalArgumentException(
148                     "{0} points Legendre-Gauss integrator not supported, " +
149                     "number of points must be in the {1}-{2} range",
150                     n, 2, 5);
151         }
152 
153     }
154 
155     /** {@inheritDoc} */
156     @Deprecated
157     public double integrate(final double min, final double max)
158         throws ConvergenceException,  FunctionEvaluationException, IllegalArgumentException {
159         return integrate(f, min, max);
160     }
161 
162     /** {@inheritDoc} */
163     public double integrate(final UnivariateRealFunction f,
164             final double min, final double max)
165         throws ConvergenceException,  FunctionEvaluationException, IllegalArgumentException {
166         
167         clearResult();
168         verifyInterval(min, max);
169         verifyIterationCount();
170 
171         // compute first estimate with a single step
172         double oldt = stage(f, min, max, 1);
173 
174         int n = 2;
175         for (int i = 0; i < maximalIterationCount; ++i) {
176 
177             // improve integral with a larger number of steps
178             final double t = stage(f, min, max, n);
179 
180             // estimate error
181             final double delta = Math.abs(t - oldt);
182             final double limit =
183                 Math.max(absoluteAccuracy,
184                          relativeAccuracy * (Math.abs(oldt) + Math.abs(t)) * 0.5);
185 
186             // check convergence
187             if ((i + 1 >= minimalIterationCount) && (delta <= limit)) {
188                 setResult(t, i);
189                 return result;
190             }
191 
192             // prepare next iteration
193             double ratio = Math.min(4, Math.pow(delta / limit, 0.5 / abscissas.length));
194             n = Math.max((int) (ratio * n), n + 1);
195             oldt = t;
196 
197         }
198 
199         throw new MaxIterationsExceededException(maximalIterationCount);
200 
201     }
202 
203     /**
204      * Compute the n-th stage integral.
205      * @param f the integrand function
206      * @param min the lower bound for the interval
207      * @param max the upper bound for the interval
208      * @param n number of steps
209      * @return the value of n-th stage integral
210      * @throws FunctionEvaluationException if an error occurs evaluating the
211      * function
212      */
213     private double stage(final UnivariateRealFunction f,
214                          final double min, final double max, final int n)
215         throws FunctionEvaluationException {
216 
217         // set up the step for the current stage
218         final double step     = (max - min) / n;
219         final double halfStep = step / 2.0;
220 
221         // integrate over all elementary steps
222         double midPoint = min + halfStep;
223         double sum = 0.0;
224         for (int i = 0; i < n; ++i) {
225             for (int j = 0; j < abscissas.length; ++j) {
226                 sum += weights[j] * f.value(midPoint + halfStep * abscissas[j]);
227             }
228             midPoint += step;
229         }
230 
231         return halfStep * sum;
232 
233     }
234 
235 }