integral2: Numerically Evaluate Double and Triple Integrals In pracma: Practical Numerical Math Functions

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.

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) ##  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) ##  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) ##  8.37758 # 8.377579076269617

Example output

\$Q
 0.7080734

\$error
 4.309141e-19

 0.5236076
 0.5236226
 -3.247091e-08
 -4.998646e-11
 2
 5.206447
 47416.76
 8.37758

pracma documentation built on Dec. 11, 2021, 9:57 a.m.