## Documentation Center |

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
and
. `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 over , . The true value of the integral is .

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

Integrate over the triangle and . The integrand is infinite at (0,0). The true value of the integral is .

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)

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

`dblquad` | `function_handle` | `integral` | `integral2` | `integral3` | `quad` | `quadgk` | `quadl` | `quadv` | `triplequad`

Was this topic helpful?