Numerically evaluate double integral, tiled method

`q = quad2d(fun,a,b,c,d)`

[q,errbnd] = quad2d(...)

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

`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(x)\le y\le d(x)$$. `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.

Integrate $$y\mathrm{sin}(x)+x\mathrm{cos}(y)$$ 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)

Integrate $${[{(x+y)}^{1/2}{(1+x+y)}^{2}]}^{-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)

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)

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

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

Was this topic helpful?