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    *      https://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  
18  /*
19   * This is not the original file distributed by the Apache Software Foundation
20   * It has been modified by the Hipparchus project
21   */
22  package org.hipparchus.analysis.integration;
23  
24  import org.hipparchus.analysis.UnivariateFunction;
25  import org.hipparchus.analysis.integration.gauss.GaussIntegrator;
26  import org.hipparchus.analysis.integration.gauss.GaussIntegratorFactory;
27  import org.hipparchus.exception.LocalizedCoreFormats;
28  import org.hipparchus.exception.MathIllegalArgumentException;
29  import org.hipparchus.exception.MathIllegalStateException;
30  import org.hipparchus.util.FastMath;
31  
32  /**
33   * This algorithm divides the integration interval into equally-sized
34   * sub-interval and on each of them performs a
35   * <a href="http://mathworld.wolfram.com/Legendre-GaussQuadrature.html">
36   * Legendre-Gauss</a> quadrature.
37   * Because of its <em>non-adaptive</em> nature, this algorithm can
38   * converge to a wrong value for the integral (for example, if the
39   * function is significantly different from zero toward the ends of the
40   * integration interval).
41   * In particular, a change of variables aimed at estimating integrals
42   * over infinite intervals as proposed
43   * <a href="http://en.wikipedia.org/w/index.php?title=Numerical_integration#Integrals_over_infinite_intervals">
44   *  here</a> should be avoided when using this class.
45   *
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 MathIllegalArgumentException if minimal number of iterations
65       * or number of points are not strictly positive.
66       * @throws MathIllegalArgumentException 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 MathIllegalArgumentException {
75          super(relativeAccuracy, absoluteAccuracy, minimalIterationCount, maximalIterationCount);
76          if (n <= 0) {
77              throw new MathIllegalArgumentException(LocalizedCoreFormats.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 MathIllegalArgumentException if {@code n < 1}.
89       */
90      public IterativeLegendreGaussIntegrator(final int n,
91                                              final double relativeAccuracy,
92                                              final double absoluteAccuracy)
93          throws MathIllegalArgumentException {
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 MathIllegalArgumentException if minimal number of iterations
105      * is not strictly positive.
106      * @throws MathIllegalArgumentException if maximal number of iterations
107      * is smaller than or equal to the minimal number of iterations.
108      * @throws MathIllegalArgumentException if {@code n < 1}.
109      */
110     public IterativeLegendreGaussIntegrator(final int n,
111                                             final int minimalIterationCount,
112                                             final int maximalIterationCount)
113         throws MathIllegalArgumentException {
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, MathIllegalStateException {
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.increment();
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 MathIllegalStateException if the maximum number of evaluations
156      * is exceeded.
157      */
158     private double stage(final int n)
159         throws MathIllegalStateException {
160         // Function to be integrated is stored in the base class.
161         final UnivariateFunction f = new UnivariateFunction() {
162                 /** {@inheritDoc} */
163                 @Override
164                 public double value(double x)
165                     throws MathIllegalArgumentException, MathIllegalStateException {
166                     return computeObjectiveValue(x);
167                 }
168             };
169 
170         final double min = getMin();
171         final double max = getMax();
172         final double step = (max - min) / n;
173 
174         double sum = 0;
175         for (int i = 0; i < n; i++) {
176             // Integrate over each sub-interval [a, b].
177             final double a = min + i * step;
178             final double b = a + step;
179             final GaussIntegrator g = FACTORY.legendreHighPrecision(numberOfPoints, a, b);
180             sum += g.integrate(f);
181         }
182 
183         return sum;
184     }
185 }