public static double Integrate(
	Function<double, double> integrand,
	Interval range
)

Public Shared Function Integrate ( _
	integrand As Function(Of Double, Double), _
	range As Interval _
) As Double

public:
static double Integrate(
	Function<double, double>^ integrand, 
	Interval range
)

Return Value

Integral values are accurate to within about a digit of full double precision.

To do integrals over infinite regions, simply set the lower bound of the range to NegativeInfinity or the upper bound to PositiveInfinity.

Our numerical integrator uses a Gauss-Kronrod rule that can integrate efficiently, combined with an adaptive strategy that limits function evaluations to those regions required to achieve the desired accuracy.

Our integrator handles smooth functions extremely efficiently. It handles integrands with discontinuities, or discontinuities of derivatives, at the price of slightly more evaluations of the integrand. It handles oscilatory functions, so long as not too many periods contribute significantly to the integral. It can integrate logarithmic and mild power-law singularities.

Strong power-law singularities will cause the alrorighm to fail with a NonconvergenceException. This is unavoidable for essentially any double-precision numerical integrator. Consider, for example, the integrable singularity 1/√x. Since ε = ∫₀^δ x^-1/2 dx = 2 δ^1/2, points within δ ∼ 10^-16 of the end-points, which as a close as you can get to a point in double precision without being on top of it, contribute at the ε ∼ 10^-8 level to our integral, well beyond limit that nearly-full double precision requires. Said differently, to know the value of the integral to ε ∼ 10^-16 prescision, we would need to evaluate the contributions of points within δ ∼ 10^-32 of the endpoints, far closer than we can get.

If you need to evaluate an integral with such a strong singularity, make an analytic change of variable to absorb the singularity before attempting numerical integration. For example, to evaluate I = ∫₀^b f(x) x^-1/2 dx, substitute y = x^1/2 to obtain I = 2 ∫₀^√b f(y²) dy.