001 /*
002 * Licensed to the Apache Software Foundation (ASF) under one or more
003 * contributor license agreements. See the NOTICE file distributed with
004 * this work for additional information regarding copyright ownership.
005 * The ASF licenses this file to You under the Apache License, Version 2.0
006 * (the "License"); you may not use this file except in compliance with
007 * the License. You may obtain a copy of the License at
008 *
009 * http://www.apache.org/licenses/LICENSE-2.0
010 *
011 * Unless required by applicable law or agreed to in writing, software
012 * distributed under the License is distributed on an "AS IS" BASIS,
013 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
014 * See the License for the specific language governing permissions and
015 * limitations under the License.
016 */
017 package org.apache.commons.math3.analysis.integration;
018
019 import org.apache.commons.math3.analysis.UnivariateFunction;
020 import org.apache.commons.math3.analysis.integration.gauss.GaussIntegratorFactory;
021 import org.apache.commons.math3.analysis.integration.gauss.GaussIntegrator;
022 import org.apache.commons.math3.exception.MaxCountExceededException;
023 import org.apache.commons.math3.exception.NotStrictlyPositiveException;
024 import org.apache.commons.math3.exception.NumberIsTooSmallException;
025 import org.apache.commons.math3.exception.TooManyEvaluationsException;
026 import org.apache.commons.math3.util.FastMath;
027
028 /**
029 * This algorithm divides the integration interval into equally-sized
030 * sub-interval and on each of them performs a
031 * <a href="http://mathworld.wolfram.com/Legendre-GaussQuadrature.html">
032 * Legendre-Gauss</a> quadrature.
033 *
034 * @version $Id: IterativeLegendreGaussIntegrator.java 1416643 2012-12-03 19:37:14Z tn $
035 * @since 3.1
036 */
037
038 public class IterativeLegendreGaussIntegrator
039 extends BaseAbstractUnivariateIntegrator {
040 /** Factory that computes the points and weights. */
041 private static final GaussIntegratorFactory FACTORY
042 = new GaussIntegratorFactory();
043 /** Number of integration points (per interval). */
044 private final int numberOfPoints;
045
046 /**
047 * Builds an integrator with given accuracies and iterations counts.
048 *
049 * @param n Number of integration points.
050 * @param relativeAccuracy Relative accuracy of the result.
051 * @param absoluteAccuracy Absolute accuracy of the result.
052 * @param minimalIterationCount Minimum number of iterations.
053 * @param maximalIterationCount Maximum number of iterations.
054 * @throws NotStrictlyPositiveException if minimal number of iterations
055 * is not strictly positive.
056 * @throws NumberIsTooSmallException if maximal number of iterations
057 * is smaller than or equal to the minimal number of iterations.
058 */
059 public IterativeLegendreGaussIntegrator(final int n,
060 final double relativeAccuracy,
061 final double absoluteAccuracy,
062 final int minimalIterationCount,
063 final int maximalIterationCount)
064 throws NotStrictlyPositiveException, NumberIsTooSmallException {
065 super(relativeAccuracy, absoluteAccuracy, minimalIterationCount, maximalIterationCount);
066 numberOfPoints = n;
067 }
068
069 /**
070 * Builds an integrator with given accuracies.
071 *
072 * @param n Number of integration points.
073 * @param relativeAccuracy Relative accuracy of the result.
074 * @param absoluteAccuracy Absolute accuracy of the result.
075 */
076 public IterativeLegendreGaussIntegrator(final int n,
077 final double relativeAccuracy,
078 final double absoluteAccuracy) {
079 this(n, relativeAccuracy, absoluteAccuracy,
080 DEFAULT_MIN_ITERATIONS_COUNT, DEFAULT_MAX_ITERATIONS_COUNT);
081 }
082
083 /**
084 * Builds an integrator with given iteration counts.
085 *
086 * @param n Number of integration points.
087 * @param minimalIterationCount Minimum number of iterations.
088 * @param maximalIterationCount Maximum number of iterations.
089 * @throws NotStrictlyPositiveException if minimal number of iterations
090 * is not strictly positive.
091 * @throws NumberIsTooSmallException if maximal number of iterations
092 * is smaller than or equal to the minimal number of iterations.
093 */
094 public IterativeLegendreGaussIntegrator(final int n,
095 final int minimalIterationCount,
096 final int maximalIterationCount) {
097 this(n, DEFAULT_RELATIVE_ACCURACY, DEFAULT_ABSOLUTE_ACCURACY,
098 minimalIterationCount, maximalIterationCount);
099 }
100
101 /** {@inheritDoc} */
102 @Override
103 protected double doIntegrate()
104 throws TooManyEvaluationsException, MaxCountExceededException {
105 // Compute first estimate with a single step.
106 double oldt = stage(1);
107
108 int n = 2;
109 while (true) {
110 // Improve integral with a larger number of steps.
111 final double t = stage(n);
112
113 // Estimate the error.
114 final double delta = FastMath.abs(t - oldt);
115 final double limit =
116 FastMath.max(getAbsoluteAccuracy(),
117 getRelativeAccuracy() * (FastMath.abs(oldt) + FastMath.abs(t)) * 0.5);
118
119 // check convergence
120 if (iterations.getCount() + 1 >= getMinimalIterationCount() &&
121 delta <= limit) {
122 return t;
123 }
124
125 // Prepare next iteration.
126 final double ratio = FastMath.min(4, FastMath.pow(delta / limit, 0.5 / numberOfPoints));
127 n = FastMath.max((int) (ratio * n), n + 1);
128 oldt = t;
129 iterations.incrementCount();
130 }
131 }
132
133 /**
134 * Compute the n-th stage integral.
135 *
136 * @param n Number of steps.
137 * @return the value of n-th stage integral.
138 * @throws TooManyEvaluationsException if the maximum number of evaluations
139 * is exceeded.
140 */
141 private double stage(final int n)
142 throws TooManyEvaluationsException {
143 // Function to be integrated is stored in the base class.
144 final UnivariateFunction f = new UnivariateFunction() {
145 public double value(double x) {
146 return computeObjectiveValue(x);
147 }
148 };
149
150 final double min = getMin();
151 final double max = getMax();
152 final double step = (max - min) / n;
153
154 double sum = 0;
155 for (int i = 0; i < n; i++) {
156 // Integrate over each sub-interval [a, b].
157 final double a = min + i * step;
158 final double b = a + step;
159 final GaussIntegrator g = FACTORY.legendreHighPrecision(numberOfPoints, a, b);
160 sum += g.integrate(f);
161 }
162
163 return sum;
164 }
165 }