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