# 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

 ```1 2 3 4 5 6``` ```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`.

### Value

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

### Note

To avoid recursion, a possibly large matrix will be used and passed between subprograms. A more efficient implementation may be possible.

### 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.

### References

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

`integral`, `cubature:adaptIntegrate`

### Examples

 ``` 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53``` ```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 ```

Search within the pracma package
Search all R packages, documentation and source code

Questions? Problems? Suggestions? or email at ian@mutexlabs.com.

Please suggest features or report bugs with the GitHub issue tracker.

All documentation is copyright its authors; we didn't write any of that.