# Documentation

Numerically evaluate double integral, tiled method

## Syntax

`q = quad2d(fun,a,b,c,d)[q,errbnd] = quad2d(...)q = quad2d(fun,a,b,c,d,param1,val1,param2,val2,...)`

## Description

`q = quad2d(fun,a,b,c,d)` approximates the integral of `fun(x,y)` over the planar region $a\le x\le b$ and $c\left(x\right)\le y\le d\left(x\right)$. `fun` is a function handle, `c` and `d` may each be a scalar or a function handle.

All input functions must be vectorized. The function `Z=fun(X,Y)` must accept 2-D matrices `X` and `Y` of the same size and return a matrix `Z` of corresponding values. The functions `ymin=c(X)` and `ymax=d(X)` must accept matrices and return matrices of the same size with corresponding values.

`[q,errbnd] = quad2d(...)`. `errbnd` is an approximate upper bound on the absolute error, `|Q - I|`, where `I` denotes the exact value of the integral.

`q = quad2d(fun,a,b,c,d,param1,val1,param2,val2,...)` performs the integration as above with specified values of optional parameters:

 `AbsTol` absolute error tolerance `RelTol` relative error tolerance

`quad2d` 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. A default tolerance value is used when a tolerance is not specified. The default value of `AbsTol` is 1e-5. The default value of `RelTol` is `100*eps(class(Q))`. This is also the minimum value of `RelTol`. Smaller `RelTol` values are automatically increased to the default value.

 `MaxFunEvals` Maximum allowed number of evaluations of `fun` reached.

The `MaxFunEvals` parameter limits the number of vectorized calls to `fun`. The default is 2000.

 `FailurePlot` Generate a plot if `MaxFunEvals` is reached.

Setting `FailurePlot` to `true` generates a graphical representation of the regions needing further refinement when `MaxFunEvals` is reached. No plot is generated if the integration succeeds before reaching `MaxFunEvals`. These (generally) 4-sided regions are mapped to rectangles internally. Clusters of small regions indicate the areas of difficulty. The default is `false`.

 `Singular` Problem may have boundary singularities

With `Singular` set to `true`, `quad2d` will employ transformations to weaken boundary singularities for better performance. The default is `true`. Setting `Singular` to `false` will turn these transformations off, which may provide a performance benefit on some smooth problems.

## Examples

### Example 1

Integrate $y\mathrm{sin}\left(x\right)+x\mathrm{cos}\left(y\right)$ over $\pi \le x\le 2\pi$, $0\le y\le \pi$. The true value of the integral is $-{\pi }^{2}$.

`Q = quad2d(@(x,y) y.*sin(x)+x.*cos(y),pi,2*pi,0,pi)`

### Example 2

Integrate ${\left[{\left(x+y\right)}^{1/2}{\left(1+x+y\right)}^{2}\right]}^{-1}$ over the triangle $0\le x\le 1$ and $0\le y\le 1-x$. The integrand is infinite at (0,0). The true value of the integral is $\pi /4-1/2$.

`fun = @(x,y) 1./(sqrt(x + y) .* (1 + x + y).^2 )`

In Cartesian coordinates:

```ymax = @(x) 1 - x; Q = quad2d(fun,0,1,0,ymax)```
In polar coordinates:
```polarfun = @(theta,r) fun(r.*cos(theta),r.*sin(theta)).*r; rmax = @(theta) 1./(sin(theta) + cos(theta)); Q = quad2d(polarfun,0,pi/2,0,rmax)```

## Limitations

`quad2d` begins by mapping the region of integration to a rectangle. Consequently, it may have trouble integrating over a region that does not have four sides or has a side that cannot be mapped smoothly to a straight line. If the integration is unsuccessful, some helpful tactics are leaving `Singular` set to its default value of `true`, changing between Cartesian and polar coordinates, or breaking the region of integration into pieces and adding the results of integration over the pieces.

For instance:

```fun = @(x,y)abs(x.^2 + y.^2 - 0.25); c = @(x)-sqrt(1 - x.^2); d = @(x)sqrt(1 - x.^2); quad2d(fun,-1,1,c,d,'AbsTol',1e-8,... 'FailurePlot',true,'Singular',false); ```
```Warning: Reached the maximum number of function evaluations (2000). The result fails the global error test. ```

The failure plot shows two areas of difficulty, near the points `(-1,0)` and `(1,0)` and near the circle .

Changing the value of `Singular` to `true` will cope with the geometric singularities at `(-1,0)` and `(1,0)`. The larger shaded areas may need refinement but are probably not areas of difficulty.

```Q = quad2d(fun,-1,1,c,d,'AbsTol',1e-8, ... 'FailurePlot',true,'Singular',true); ```
```Warning: Reached the maximum number of function evaluations (2000). The result passes the global error test. ```

From here you can take advantage of symmetry:

```Q = 4*quad2d(fun,0,1,0,d,'Abstol',1e-8,... 'Singular',true,'FailurePlot',true) ```
```Q = 0.9817 ```

However, the code is still working very hard near the singularity. It may not be able to provide higher accuracy:

```Q = 4*quad2d(fun,0,1,0,d,'Abstol',1e-10,... 'Singular',true,'FailurePlot',true); ```
```Warning: Reached the maximum number of function evaluations (2000). The result passes the global error test. ```

At higher accuracy, a change in coordinates may work better.

```polarfun = @(theta,r) fun(r.*cos(theta),r.*sin(theta)).*r; Q = 4*quad2d(polarfun,0,pi/2,0,1,'AbsTol',1e-10); ```

It is best to put the singularity on the boundary by splitting the region of integration into two parts:

```Q1 = 4*quad2d(polarfun,0,pi/2,0,0.5,'AbsTol',5e-11); Q2 = 4*quad2d(polarfun,0,pi/2,0.5,1,'AbsTol',5e-11); Q = Q1 + Q2; ```

## References

[1] L.F. Shampine, "Matlab Program for Quadrature in 2D." Applied Mathematics and Computation. Vol. 202, Issue 1, 2008, pp. 266–274.