Documentation

quad2d

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 offun(x,y)over the planar region a x b and c ( x ) y d ( x ) .funis a function handle,canddmay each be a scalar or a function handle.

All input functions must be vectorized. The functionZ=fun(X,Y)must accept 2-D matricesXandYof the same size and return a matrixZof corresponding values. The functionsymin=c(X)andymax=d(X)must accept matrices and return matrices of the same size with corresponding values.

[q,errbnd] = quad2d(...).errbndis an approximate upper bound on the absolute error,|Q - I|, whereIdenotes 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

quad2dattempts to satisfyERRBND <= 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 ofAbsTolis 1e-5. The default value ofRelTolis100*eps(class(Q)). This is also the minimum value ofRelTol. SmallerRelTolvalues are automatically increased to the default value.

MaxFunEvals Maximum allowed number of evaluations offunreached.

TheMaxFunEvalsparameter limits the number of vectorized calls tofun. The default is 2000.

FailurePlot Generate a plot ifMaxFunEvalsis reached.

SettingFailurePlottotruegenerates a graphical representation of the regions needing further refinement whenMaxFunEvalsis reached. No plot is generated if the integration succeeds before reachingMaxFunEvals. These (generally) 4-sided regions are mapped to rectangles internally. Clusters of small regions indicate the areas of difficulty. The default isfalse.

Singular Problem may have boundary singularities

WithSingularset totrue,quad2dwill employ transformations to weaken boundary singularities for better performance. The default istrue. SettingSingulartofalsewill turn these transformations off, which may provide a performance benefit on some smooth problems.

Examples

collapse all

Integrate

overand.

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

Compare the result to the true value of the integral,.

-pi^2
ans = -9.8696

Integrate the function

over the regionand. This integrand is infinite at the origin (0,0), which lies on the boundary of the integration region.

fun = @(x,y) 1./(sqrt(x + y) .* (1 + x + y).^2 ); ymax = @(x) 1 - x; Q = quad2d(fun,0,1,0,ymax)
Q = 0.2854

The true value of the integral is.

pi/4 - 0.5
ans = 0.2854

quad2dbegins 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 leavingSingularset to its default value oftrue, 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 ofSingulartotruewill 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',真正的);
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',真正的)
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',真正的);
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.

Extended Capabilities

Was this topic helpful?