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.analysis.UnivariateFunction;
20  import org.apache.commons.math3.analysis.integration.gauss.GaussIntegratorFactory;
21  import org.apache.commons.math3.analysis.integration.gauss.GaussIntegrator;
22  import org.apache.commons.math3.exception.MathIllegalArgumentException;
23  import org.apache.commons.math3.exception.MaxCountExceededException;
24  import org.apache.commons.math3.exception.NotStrictlyPositiveException;
25  import org.apache.commons.math3.exception.NumberIsTooSmallException;
26  import org.apache.commons.math3.exception.TooManyEvaluationsException;
27  import org.apache.commons.math3.exception.util.LocalizedFormats;
28  import org.apache.commons.math3.util.FastMath;
29  
30  /**
31   * This algorithm divides the integration interval into equally-sized
32   * sub-interval and on each of them performs a
33   * <a href="http://mathworld.wolfram.com/Legendre-GaussQuadrature.html">
34   * Legendre-Gauss</a> quadrature.
35   * Because of its <em>non-adaptive</em> nature, this algorithm can
36   * converge to a wrong value for the integral (for example, if the
37   * function is significantly different from zero toward the ends of the
38   * integration interval).
39   * In particular, a change of variables aimed at estimating integrals
40   * over infinite intervals as proposed
41   * <a href="http://en.wikipedia.org/w/index.php?title=Numerical_integration#Integrals_over_infinite_intervals">
42   *  here</a> should be avoided when using this class.
43   *
44   * @since 3.1
45   */
46  
47  public class IterativeLegendreGaussIntegrator
48      extends BaseAbstractUnivariateIntegrator {
49      /** Factory that computes the points and weights. */
50      private static final GaussIntegratorFactory FACTORY
51          = new GaussIntegratorFactory();
52      /** Number of integration points (per interval). */
53      private final int numberOfPoints;
54  
55      /**
56       * Builds an integrator with given accuracies and iterations counts.
57       *
58       * @param n Number of integration points.
59       * @param relativeAccuracy Relative accuracy of the result.
60       * @param absoluteAccuracy Absolute accuracy of the result.
61       * @param minimalIterationCount Minimum number of iterations.
62       * @param maximalIterationCount Maximum number of iterations.
63       * @throws NotStrictlyPositiveException if minimal number of iterations
64       * or number of points are not strictly positive.
65       * @throws NumberIsTooSmallException if maximal number of iterations
66       * is smaller than or equal to the minimal number of iterations.
67       */
68      public IterativeLegendreGaussIntegrator(final int n,
69                                              final double relativeAccuracy,
70                                              final double absoluteAccuracy,
71                                              final int minimalIterationCount,
72                                              final int maximalIterationCount)
73          throws NotStrictlyPositiveException, NumberIsTooSmallException {
74          super(relativeAccuracy, absoluteAccuracy, minimalIterationCount, maximalIterationCount);
75          if (n <= 0) {
76              throw new NotStrictlyPositiveException(LocalizedFormats.NUMBER_OF_POINTS, n);
77          }
78         numberOfPoints = n;
79      }
80  
81      /**
82       * Builds an integrator with given accuracies.
83       *
84       * @param n Number of integration points.
85       * @param relativeAccuracy Relative accuracy of the result.
86       * @param absoluteAccuracy Absolute accuracy of the result.
87       * @throws NotStrictlyPositiveException if {@code n < 1}.
88       */
89      public IterativeLegendreGaussIntegrator(final int n,
90                                              final double relativeAccuracy,
91                                              final double absoluteAccuracy)
92          throws NotStrictlyPositiveException {
93          this(n, relativeAccuracy, absoluteAccuracy,
94               DEFAULT_MIN_ITERATIONS_COUNT, DEFAULT_MAX_ITERATIONS_COUNT);
95      }
96  
97      /**
98       * Builds an integrator with given iteration counts.
99       *
100      * @param n Number of integration points.
101      * @param minimalIterationCount Minimum number of iterations.
102      * @param maximalIterationCount Maximum number of iterations.
103      * @throws NotStrictlyPositiveException if minimal number of iterations
104      * is not strictly positive.
105      * @throws NumberIsTooSmallException if maximal number of iterations
106      * is smaller than or equal to the minimal number of iterations.
107      * @throws NotStrictlyPositiveException if {@code n < 1}.
108      */
109     public IterativeLegendreGaussIntegrator(final int n,
110                                             final int minimalIterationCount,
111                                             final int maximalIterationCount)
112         throws NotStrictlyPositiveException, NumberIsTooSmallException {
113         this(n, DEFAULT_RELATIVE_ACCURACY, DEFAULT_ABSOLUTE_ACCURACY,
114              minimalIterationCount, maximalIterationCount);
115     }
116 
117     /** {@inheritDoc} */
118     @Override
119     protected double doIntegrate()
120         throws MathIllegalArgumentException, TooManyEvaluationsException, MaxCountExceededException {
121         // Compute first estimate with a single step.
122         double oldt = stage(1);
123 
124         int n = 2;
125         while (true) {
126             // Improve integral with a larger number of steps.
127             final double t = stage(n);
128 
129             // Estimate the error.
130             final double delta = FastMath.abs(t - oldt);
131             final double limit =
132                 FastMath.max(getAbsoluteAccuracy(),
133                              getRelativeAccuracy() * (FastMath.abs(oldt) + FastMath.abs(t)) * 0.5);
134 
135             // check convergence
136             if (iterations.getCount() + 1 >= getMinimalIterationCount() &&
137                 delta <= limit) {
138                 return t;
139             }
140 
141             // Prepare next iteration.
142             final double ratio = FastMath.min(4, FastMath.pow(delta / limit, 0.5 / numberOfPoints));
143             n = FastMath.max((int) (ratio * n), n + 1);
144             oldt = t;
145             iterations.incrementCount();
146         }
147     }
148 
149     /**
150      * Compute the n-th stage integral.
151      *
152      * @param n Number of steps.
153      * @return the value of n-th stage integral.
154      * @throws TooManyEvaluationsException if the maximum number of evaluations
155      * is exceeded.
156      */
157     private double stage(final int n)
158         throws TooManyEvaluationsException {
159         // Function to be integrated is stored in the base class.
160         final UnivariateFunction f = new UnivariateFunction() {
161                 public double value(double x)
162                     throws MathIllegalArgumentException, TooManyEvaluationsException {
163                     return computeObjectiveValue(x);
164                 }
165             };
166 
167         final double min = getMin();
168         final double max = getMax();
169         final double step = (max - min) / n;
170 
171         double sum = 0;
172         for (int i = 0; i < n; i++) {
173             // Integrate over each sub-interval [a, b].
174             final double a = min + i * step;
175             final double b = a + step;
176             final GaussIntegrator g = FACTORY.legendreHighPrecision(numberOfPoints, a, b);
177             sum += g.integrate(f);
178         }
179 
180         return sum;
181     }
182 }