# Numerically Evaluate Double and Triple Integrals Numerically evaluate a double integral, resp. a triple integral by reducing it to a double integral. ```r 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. ## 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|)`. 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 <= b` -- 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` also accepts functions for the inner interval limits. `ymin, ymax` must be constants or functions of one variable (`x`), `zmin, zmax` constants or functions of two variables (`x, y`), all functions vectorized. The triple integral will be first integrated over the second and third variable with `integral2`, and then integrated over a single variable with `integral`. ## Returns Returns a list with `Q` the integral and `error` the error term. ## References Shampine, L. F. (2008). MATLAB Program for Quadrature in 2D. Proceedings of Applied Mathematics and Computation, 2008, pp. 266--274. ## Author(s) Copyright (c) 2008 Lawrence F. Shampine for Matlab code and description of the program; adapted and converted to R by Hans W Borchers. ## Note To avoid recursion, a possibly large matrix will be used and passed between subprograms. A more efficient implementation may be possible. ## See Also `integral`, `cubature:adaptIntegrate` ## Examples ```r 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 ## Examples of computing triple integrals f0 <- function(x, y, z) y*sin(x) + z*cos(x) integral3(f0, 0, pi, 0,1, -1,1) # - 2.0 => 0.0 f1 <- function(x, y, z) exp(x+y+z) integral3(f1, 0, 1, 1, 2, 0, 0.5) ## [1] 5.206447 # 5.20644655 f2 <- function(x, y, z) x^2 + y^2 + z a <- 2; b <- 4 ymin <- function(x) x - 1 ymax <- function(x) x + 6 zmin <- -2 zmax <- function(x, y) 4 + y^2 integral3(f2, a, b, ymin, ymax, zmin, zmax) ## [1] 47416.75556 # 47416.7555556 f3 <- function(x, y, z) sqrt(x^2 + y^2) a <- -2; b <- 2 ymin <- function(x) -sqrt(4-x^2) ymax <- function(x) sqrt(4-x^2) zmin <- function(x, y) sqrt(x^2 + y^2) zmax <- 2 integral3(f3, a, b, ymin, ymax, zmin, zmax) ## [1] 8.37758 # 8.377579076269617 ```

integral2 function