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
9
* http://www.apache.org/licenses/LICENSE-2.0
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.
17
package org.apache.commons.math.analysis;
19
import org.apache.commons.math.FunctionEvaluationException;
20
import org.apache.commons.math.MaxIterationsExceededException;
23
* Implements the <a href="http://mathworld.wolfram.com/RombergIntegration.html">
24
* Romberg Algorithm</a> for integration of real univariate functions. For
25
* reference, see <b>Introduction to Numerical Analysis</b>, ISBN 038795452X,
28
* Romberg integration employs k successvie refinements of the trapezoid
29
* rule to remove error terms less than order O(N^(-2k)). Simpson's rule
30
* is a special case of k = 2.</p>
32
* @version $Revision: 620312 $ $Date: 2008-02-10 12:28:59 -0700 (Sun, 10 Feb 2008) $
35
public class RombergIntegrator extends UnivariateRealIntegratorImpl {
37
/** serializable version identifier */
38
private static final long serialVersionUID = -1058849527738180243L;
41
* Construct an integrator for the given function.
43
* @param f function to integrate
45
public RombergIntegrator(UnivariateRealFunction f) {
50
* Integrate the function in the given interval.
52
* @param min the lower bound for the interval
53
* @param max the upper bound for the interval
54
* @return the value of integral
55
* @throws MaxIterationsExceededException if the maximum iteration count is exceeded
56
* or the integrator detects convergence problems otherwise
57
* @throws FunctionEvaluationException if an error occurs evaluating the
59
* @throws IllegalArgumentException if any parameters are invalid
61
public double integrate(double min, double max) throws MaxIterationsExceededException,
62
FunctionEvaluationException, IllegalArgumentException {
64
int i = 1, j, m = maximalIterationCount + 1;
65
// Array strcture here can be improved for better space
66
// efficiency because only the lower triangle is used.
67
double r, t[][] = new double[m][m], s, olds;
70
verifyInterval(min, max);
71
verifyIterationCount();
73
TrapezoidIntegrator qtrap = new TrapezoidIntegrator(this.f);
74
t[0][0] = qtrap.stage(min, max, 0);
76
while (i <= maximalIterationCount) {
77
t[i][0] = qtrap.stage(min, max, i);
78
for (j = 1; j <= i; j++) {
79
// Richardson extrapolation coefficient
80
r = (1L << (2 * j)) -1;
81
t[i][j] = t[i][j-1] + (t[i][j-1] - t[i-1][j-1]) / r;
84
if (i >= minimalIterationCount) {
85
if (Math.abs(s - olds) <= Math.abs(relativeAccuracy * olds)) {
93
throw new MaxIterationsExceededException(maximalIterationCount);
97
* Verifies that the iteration limits are valid and within the range.
99
* @throws IllegalArgumentException if not
101
protected void verifyIterationCount() throws IllegalArgumentException {
102
super.verifyIterationCount();
103
// at most 32 bisection refinements due to higher order divider
104
if (maximalIterationCount > 32) {
105
throw new IllegalArgumentException
106
("Iteration upper limit out of [0, 32] range: " +
107
maximalIterationCount);