/* XXL: The eXtensible and fleXible Library for data processing
Copyright (C) 2000-2011 Prof. Dr. Bernhard Seeger
Head of the Database Research Group
Department of Mathematics and Computer Science
University of Marburg
Germany
This library is free software; you can redistribute it and/or
modify it under the terms of the GNU Lesser General Public
License as published by the Free Software Foundation; either
version 3 of the License, or (at your option) any later version.
This library is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
Lesser General Public License for more details.
You should have received a copy of the GNU Lesser General Public
License along with this library; If not, see <http://www.gnu.org/licenses/>.
http://code.google.com/p/xxl/
*/
package xxl.core.math.numerics.integration;
import xxl.core.functions.Function;
import xxl.core.math.Maths;
import xxl.core.math.functions.RealFunction;
import xxl.core.math.functions.RealFunctionArea;
import xxl.core.math.functions.RealFunctionFunction;
/** This class provides the numerical integration of a
* {@link xxl.core.math.functions.RealFunction real-valued function}
* that is defined over an interval <tt>[a,b]</tt>. For the integration
* the <tt>Simpson's Rule</tt> is applied. The <tt>Simpson's Rule</tt> belongs to
* the <tt>Newton-Cotes-Formulas</tt> and bases on three weights, providing
* an approximation of order 4. By partitioning the interval, the <tt>Simpson's Rule</tt>
* can be applied iteratively. Thereby, <tt>n</tt>, the number of intervals,
* has to be even.
*
* @see xxl.core.math.functions.RealFunction
* @see xxl.core.math.functions.RealFunctionArea
* @see xxl.core.math.numerics.integration.TrapezoidalRuleRealFunctionArea
*/
public class SimpsonsRuleRealFunctionArea extends RealFunctionArea {
/**
* error bound
*/
protected double epsilon;
/** Constructs a new Object of this type.
*
* @param realFunction {@link xxl.core.math.functions.RealFunction real-valued function} to integrate
* @param epsilon error bound
*/
public SimpsonsRuleRealFunctionArea(RealFunction realFunction, double epsilon) {
super(realFunction);
this.epsilon = epsilon;
}
/** Computes the area "under" the given real-valued function
* for the interval <tt>[a,b]</tt>.
*
* @param a left interval border of the area to compute
* @param b right interval border of the area to compute
* @return numerical approximation of the integral
* @throws IllegalArgumentException invalid arguments
*/
public double eval(double a, double b) throws IllegalArgumentException {
return SimpsonsRuleRealFunctionArea.simpsonx(a, b, new RealFunctionFunction(realFunction), epsilon);
}
/** Computes the area under a real-valued function for a given interval <tt>[a,b]</tt> using
* the Simpson's rule for numerical integration iteratively. Thereby, the number of steps
* has to be even.
*
* @param a left border of the integration interval
* @param b right border of the integration interval
* @param n number of steps used for computation
* @param real1DFunction one-dimensional real-valued function to integrate numerically
* @throws IllegalArgumentException if n is not even
* @return numerical approximation of the integral
*/
public static double simpson(double a, double b, Function real1DFunction, int n) throws IllegalArgumentException {
if (!Maths.isEven(n))
throw new IllegalArgumentException("number of computations (n) needs to be even!");
double h = (b - a) / n; // step-width
double sg = 0.0; // even index of the function values
double su = 0.0; // uneven index of the function values
double x = 0.0;
for (int i = 1; i <= n - 1; i++) {
x = a + i * h;
if (Maths.isEven(i)) {
sg += ((Number) real1DFunction.invoke(new Double(x))).doubleValue();
} else {
su += ((Number) real1DFunction.invoke(new Double(x))).doubleValue();
}
}
return (
((Number) real1DFunction.invoke(new Double(a))).doubleValue()
+ 4 * su
+ 2 * sg
+ ((Number) real1DFunction.invoke(new Double(b))).doubleValue())
* h
/ 3.0;
}
/** Computes the area under a real-valued function for a given interval <tt>[a,b]</tt> using
* the Simpson's rule for numerical integration iteratively until an error threshold is reached.
* For estimating the error, <tt>|S(n) - S (n/2)| <= epsilon * |S(n)|</tt> is used.
*
* @param a left border of the integration interval
* @param b right border of the integration interval
* @param epsilon error threshold used as |S(n) - S (n/2)| <= epsilon * |S(n)|
* @param real1DFunction real-valued (1D) function to integrate numerical
* @return numerical approximation of the integral for the given error bound
*/
public static double simpsonx(double a, double b, Function real1DFunction, double epsilon) {
int n = 2;
double sab =
((Number) real1DFunction.invoke(new Double(a))).doubleValue()
+ ((Number) real1DFunction.invoke(new Double(b))).doubleValue();
double sm = ((Number) real1DFunction.invoke(new Double(a + (b - a) * .5))).doubleValue();
double s = (sab + 4 * sm) * (b - a) / 6.0;
double s_old = 0.0;
do {
n = 2 * n;
s_old = s;
s = simpson(a, b, real1DFunction, n);
} while (Math.abs(s - s_old) > epsilon * Math.abs(s));
return s;
}
}