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   * @version $Id: IterativeLegendreGaussIntegrator.java 1499765 2013-07-04 14:24:11Z erans $
45   * @since 3.1
46   */
47  
48  public class IterativeLegendreGaussIntegrator
49      extends BaseAbstractUnivariateIntegrator {
50      /** Factory that computes the points and weights. */
51      private static final GaussIntegratorFactory FACTORY
52          = new GaussIntegratorFactory();
53      /** Number of integration points (per interval). */
54      private final int numberOfPoints;
55  
56      /**
57       * Builds an integrator with given accuracies and iterations counts.
58       *
59       * @param n Number of integration points.
60       * @param relativeAccuracy Relative accuracy of the result.
61       * @param absoluteAccuracy Absolute accuracy of the result.
62       * @param minimalIterationCount Minimum number of iterations.
63       * @param maximalIterationCount Maximum number of iterations.
64       * @throws NotStrictlyPositiveException if minimal number of iterations
65       * or number of points are not strictly positive.
66       * @throws NumberIsTooSmallException if maximal number of iterations
67       * is smaller than or equal to the minimal number of iterations.
68       */
69      public IterativeLegendreGaussIntegrator(final int n,
70                                              final double relativeAccuracy,
71                                              final double absoluteAccuracy,
72                                              final int minimalIterationCount,
73                                              final int maximalIterationCount)
74          throws NotStrictlyPositiveException, NumberIsTooSmallException {
75          super(relativeAccuracy, absoluteAccuracy, minimalIterationCount, maximalIterationCount);
76          if (n <= 0) {
77              throw new NotStrictlyPositiveException(LocalizedFormats.NUMBER_OF_POINTS, n);
78          }
79         numberOfPoints = n;
80      }
81  
82      /**
83       * Builds an integrator with given accuracies.
84       *
85       * @param n Number of integration points.
86       * @param relativeAccuracy Relative accuracy of the result.
87       * @param absoluteAccuracy Absolute accuracy of the result.
88       * @throws NotStrictlyPositiveException if {@code n < 1}.
89       */
90      public IterativeLegendreGaussIntegrator(final int n,
91                                              final double relativeAccuracy,
92                                              final double absoluteAccuracy)
93          throws NotStrictlyPositiveException {
94          this(n, relativeAccuracy, absoluteAccuracy,
95               DEFAULT_MIN_ITERATIONS_COUNT, DEFAULT_MAX_ITERATIONS_COUNT);
96      }
97  
98      /**
99       * Builds an integrator with given iteration counts.
100      *
101      * @param n Number of integration points.
102      * @param minimalIterationCount Minimum number of iterations.
103      * @param maximalIterationCount Maximum number of iterations.
104      * @throws NotStrictlyPositiveException if minimal number of iterations
105      * is not strictly positive.
106      * @throws NumberIsTooSmallException if maximal number of iterations
107      * is smaller than or equal to the minimal number of iterations.
108      * @throws NotStrictlyPositiveException if {@code n < 1}.
109      */
110     public IterativeLegendreGaussIntegrator(final int n,
111                                             final int minimalIterationCount,
112                                             final int maximalIterationCount)
113         throws NotStrictlyPositiveException, NumberIsTooSmallException {
114         this(n, DEFAULT_RELATIVE_ACCURACY, DEFAULT_ABSOLUTE_ACCURACY,
115              minimalIterationCount, maximalIterationCount);
116     }
117 
118     /** {@inheritDoc} */
119     @Override
120     protected double doIntegrate()
121         throws MathIllegalArgumentException, TooManyEvaluationsException, MaxCountExceededException {
122         // Compute first estimate with a single step.
123         double oldt = stage(1);
124 
125         int n = 2;
126         while (true) {
127             // Improve integral with a larger number of steps.
128             final double t = stage(n);
129 
130             // Estimate the error.
131             final double delta = FastMath.abs(t - oldt);
132             final double limit =
133                 FastMath.max(getAbsoluteAccuracy(),
134                              getRelativeAccuracy() * (FastMath.abs(oldt) + FastMath.abs(t)) * 0.5);
135 
136             // check convergence
137             if (iterations.getCount() + 1 >= getMinimalIterationCount() &&
138                 delta <= limit) {
139                 return t;
140             }
141 
142             // Prepare next iteration.
143             final double ratio = FastMath.min(4, FastMath.pow(delta / limit, 0.5 / numberOfPoints));
144             n = FastMath.max((int) (ratio * n), n + 1);
145             oldt = t;
146             iterations.incrementCount();
147         }
148     }
149 
150     /**
151      * Compute the n-th stage integral.
152      *
153      * @param n Number of steps.
154      * @return the value of n-th stage integral.
155      * @throws TooManyEvaluationsException if the maximum number of evaluations
156      * is exceeded.
157      */
158     private double stage(final int n)
159         throws TooManyEvaluationsException {
160         // Function to be integrated is stored in the base class.
161         final UnivariateFunction f = new UnivariateFunction() {
162                 public double value(double x)
163                     throws MathIllegalArgumentException, TooManyEvaluationsException {
164                     return computeObjectiveValue(x);
165                 }
166             };
167 
168         final double min = getMin();
169         final double max = getMax();
170         final double step = (max - min) / n;
171 
172         double sum = 0;
173         for (int i = 0; i < n; i++) {
174             // Integrate over each sub-interval [a, b].
175             final double a = min + i * step;
176             final double b = a + step;
177             final GaussIntegrator g = FACTORY.legendreHighPrecision(numberOfPoints, a, b);
178             sum += g.integrate(f);
179         }
180 
181         return sum;
182     }
183 }