Coverage Report - org.apache.commons.math.analysis.RombergIntegrator

Classes in this File	Line Coverage	Branch Coverage	Complexity
RombergIntegrator	96%	100%	3.6666666666666665;3.667

1		/*
2		* Copyright 2005 The Apache Software Foundation.
3		*
4		* Licensed under the Apache License, Version 2.0 (the "License");
5		* you may not use this file except in compliance with the License.
6		* You may obtain a copy of the License at
7		*
8		* http://www.apache.org/licenses/LICENSE-2.0
9		*
10		* Unless required by applicable law or agreed to in writing, software
11		* distributed under the License is distributed on an "AS IS" BASIS,
12		* WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
13		* See the License for the specific language governing permissions and
14		* limitations under the License.
15		*/
16		package org.apache.commons.math.analysis;
17
18		import org.apache.commons.math.ConvergenceException;
19		import org.apache.commons.math.FunctionEvaluationException;
20
21		/**
22		* Implements the <a href="http://mathworld.wolfram.com/RombergIntegration.html">
23		* Romberg Algorithm</a> for integration of real univariate functions. For
24		* reference, see <b>Introduction to Numerical Analysis</b>, ISBN 038795452X,
25		* chapter 3.
26		* <p>
27		* Romberg integration employs k successvie refinements of the trapezoid
28		* rule to remove error terms less than order O(N^(-2k)). Simpson's rule
29		* is a special case of k = 2.
30		*
31		* @version $Revision$ $Date: 2005-08-24 15:27:44 -0700 (Wed, 24 Aug 2005) $
32		*/
33		public class RombergIntegrator extends UnivariateRealIntegratorImpl {
34
35		/** serializable version identifier */
36		static final long serialVersionUID = -1058849527738180243L;
37
38		/**
39		* Construct an integrator for the given function.
40		*
41		* @param f function to integrate
42		*/
43		public RombergIntegrator(UnivariateRealFunction f) {
44	6	super(f, 32);
45	6	}
46
47		/**
48		* Integrate the function in the given interval.
49		*
50		* @param min the lower bound for the interval
51		* @param max the upper bound for the interval
52		* @return the value of integral
53		* @throws ConvergenceException if the maximum iteration count is exceeded
54		* or the integrator detects convergence problems otherwise
55		* @throws FunctionEvaluationException if an error occurs evaluating the
56		* function
57		* @throws IllegalArgumentException if any parameters are invalid
58		*/
59		public double integrate(double min, double max) throws ConvergenceException,
60		FunctionEvaluationException, IllegalArgumentException {
61
62	16	int i = 1, j, m = maximalIterationCount + 1;
63		// Array strcture here can be improved for better space
64		// efficiency because only the lower triangle is used.
65	16	double r, t[][] = new double[m][m], s, olds;
66
67	16	clearResult();
68	16	verifyInterval(min, max);
69	14	verifyIterationCount();
70
71	10	TrapezoidIntegrator qtrap = new TrapezoidIntegrator(this.f);
72	10	t[0][0] = qtrap.stage(min, max, 0);
73	10	olds = t[0][0];
74	34	while (i <= maximalIterationCount) {
75	34	t[i][0] = qtrap.stage(min, max, i);
76	112	for (j = 1; j <= i; j++) {
77		// Richardson extrapolation coefficient
78	78	r = (1L << (2 * j)) -1;
79	78	t[i][j] = t[i][j-1] + (t[i][j-1] - t[i-1][j-1]) / r;
80		}
81	34	s = t[i][i];
82	34	if (i >= minimalIterationCount) {
83	14	if (Math.abs(s - olds) <= Math.abs(relativeAccuracy * olds)) {
84	10	setResult(s, i);
85	10	return result;
86		}
87		}
88	24	olds = s;
89	24	i++;
90		}
91	0	throw new ConvergenceException("Maximum number of iterations exceeded.");
92		}
93
94		/**
95		* Verifies that the iteration limits are valid and within the range.
96		*
97		* @throws IllegalArgumentException if not
98		*/
99		protected void verifyIterationCount() throws IllegalArgumentException {
100	20	super.verifyIterationCount();
101		// at most 32 bisection refinements due to higher order divider
102	18	if (maximalIterationCount > 32) {
103	2	throw new IllegalArgumentException
104		("Iteration upper limit out of [0, 32] range: " +
105		maximalIterationCount);
106		}
107	16	}
108		}