integral2: Numerically Evaluate Double and Triple Integrals

Description

Numerically evaluate a double integral, resp. a triple integral by reducing it to a double integral.

Usage

integral2(fun, xmin, xmax, ymin, ymax, sector = FALSE,
            reltol = 1e-6, abstol = 0, maxlist = 5000,
            singular = FALSE, vectorized = TRUE, ...)
integral3(fun, xmin, xmax, ymin, ymax, zmin, zmax,
            reltol = 1e-6, ...)

Arguments

fun

function

xmin, xmax

lower and upper limits of x.

ymin, ymax

lower and upper limits of y.

zmin, zmax

lower and upper limits of z.

sector

logical.

reltol

relative tolerance.

abstol

absolute tolerance.

maxlist

maximum length of the list of rectangles.

singular

logical; are there singularities at vertices.

vectorized

logical; is the function fully vectorized.

...

additional parameters to be passed to the function.

Value

Returns a list with Q the integral and error the error term.

Details

integral2 implements the `TwoD' algorithm, that is Gauss-Kronrod with (3, 7)-nodes on 2D rectangles.

The borders of the domain of integration must be finite. The limits of y, that is ymin and ymax, can be constants or scalar functions of x that describe the lower and upper boundaries. These functions must be vectorized.

integral2 attempts to satisfy ERRBND <= max(abstol,reltol*|q|)<="" code="">. This is absolute error control when |Q| is sufficiently small and relative error control when |Q| is larger.

The function fun itself must be fully vectorized: It must accept arrays X and Y and return an array Z = f(X,Y) of corresponding values. If option vectorized is set to FALSE the procedure will enforce this vectorized behavior.

With sector=TRUE the region is a generalized sector that is described in polar coordinates (r,theta) by

0 <= a="" <="theta"> -- a and b must be constants c <= r="" <="d -- c and d can be constants or ...

... functions of theta that describe the lower and upper boundaries. Functions must be vectorized. NOTE Polar coordinates are used only to describe the region -- the integrand is f(x,y) for both kinds of regions.

integral2 can be applied to functions that are singular on a boundary. With value singular=TRUE, this option causes integral2 to use transformations to weaken singularities for better performance.

integral3 reduces the triple integral by first integrating over the second and third variable with integral2, and then integrating over a single variable with integral.

References

Shampine, L. F. (2008). MATLAB Program for Quadrature in 2D. Proceedings of Applied Mathematics and Computation, 2008, pp. 266--274.

Examples

Run this code

fun <- function(x, y) cos(x) * cos(y)
integral2(fun, 0, 1, 0, 1, reltol = 1e-10)
# $Q:     0.708073418273571  # 0.70807341827357119350 = sin(1)^2
# $error: 8.618277e-19       # 1.110223e-16

##  Compute the volume of a sphere
f <- function(x, y) sqrt(1 -x^2 - y^2)
xmin <- 0; xmax <- 1
ymin <- 0; ymax <- function(x) sqrt(1 - x^2)
I <- integral2(f, xmin, xmax, ymin, ymax)
I$Q                             # 0.5236076 - pi/6 => 8.800354e-06

##  Compute the volume over a sector
I <- integral2(f, 0,pi/2, 0,1, sector = TRUE)
I$Q                             # 0.5236308 - pi/6 => 3.203768e-05

##  Integrate 1/( sqrt(x + y)*(1 + x + y)^2 ) over the triangle
##   0 <= x <= 1, 0 <= y <= 1 - x.  The integrand is infinite at (0,0).
f <- function(x,y) 1/( sqrt(x + y) * (1 + x + y)^2 )
ymax <- function(x) 1 - x
I <- integral2(f, 0,1, 0,ymax)
I$Q + 1/2 - pi/4                # -3.247091e-08

##  Compute this integral as a sector
rmax <- function(theta) 1/(sin(theta) + cos(theta))
I <- integral2(f, 0,pi/2, 0,rmax, sector = TRUE, singular = TRUE)
I$Q + 1/2 - pi/4                # -4.998646e-11

##  Compute a triple integral
fun <- function(x, y, z) y*sin(x) + z*cos(x)
integral3(fun, 0,pi, 0,1, -1,1) # - 2.0 => 0.0

Run the code above in your browser using DataLab

Get 50% off unlimited learning