Numerically Evaluate Double and Triple Integrals
Numerically evaluate a double integral, resp. a triple integral by reducing it to a double integral.
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, ...)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. |
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 constantsc <= 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 a list with Q the integral and error the error term.
To avoid recursion, a possibly large matrix will be used and passed between subprograms. A more efficient implementation may be possible.
Copyright (c) 2008 Lawrence F. Shampine for Matlab code and description of the program; adapted and converted to R by Hans W Borchers.
Shampine, L. F. (2008). MATLAB Program for Quadrature in 2D. Proceedings of Applied Mathematics and Computation, 2008, pp. 266–274.
integral, cubature:adaptIntegrate
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
Please choose more modern alternatives, such as Google Chrome or Mozilla Firefox.