Documentation

This is machine translation

Translated by
Mouse over text to see original. Click the button below to return to the English verison of the page.

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; ```