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