计算量不小 区间不宜过大

zlabel('f5')

title('f5')

......................

...........................

matlab求二重积分

g = @(x,y)

1.*((x-2.5).^2+(y-1).^2<=1)+0.*((x-2.5).^2+(y-1).^2>1);

f = @(x,y)

1.*((x-2).^2+(y-1).^2<=1)+0.*((x-2).^2+(y-1).^2>1);

z = @(x,y) f(x,y)+g(x,y)

Q = dblquad(z,0,4,-1,3)

...........

% 以下代码在7.1版以上均可运行。

r=1;

a=2; % 输入a的值

Phi=3; % 输入Phi的值

f1 = @(Theta_2,Phi_2) sin(Theta_2).*sin(Theta_2).*cos(Phi_2);

f2 = @(Theta_2,Phi_2) sqrt(r^2+a^2-2*r*a.*sin(Theta_2).*cos(Phi-Phi_2));

f3 = @(Theta_2,Phi_2) f1(Theta_2,Phi_2)./f2(Theta_2,Phi_2);

f5 = dblquad(f3,0,pi,0,pi)

......................

...........................

..................

dblquad-

Numerically evaluate double

integral over rectangle

Syntax

q = dblquad(fun,xmin,xmax,ymin,ymax)

q = dblquad(fun,xmin,xmax,ymin,ymax,tol)

q = dblquad(fun,xmin,xmax,ymin,ymax,tol,method)

Description

q = dblquad(fun,xmin,xmax,ymin,ymax)

callsthe quad

function to evaluatethe double integral fun(x,y) over the

rectangle xmin<= x <= xmax,

ymin <=y <= ymax.

fun isa function handle. See Function

Handles in the MATLAB Programming documentationfor more

information. fun(x,y) must accept a vector x anda

scalar y and return a vector of values of

theintegrand.

Parameterizing

Functions, in the MATLAB Mathematicsdocumentation, explains how

to provide additional parameters to thefunction fun, if

necessary.

q = dblquad(fun,xmin,xmax,ymin,ymax,tol)

usesa tolerance tol instead of the default, which is

1.0e-6.

q =

dblquad(fun,xmin,xmax,ymin,ymax,tol,method) usesthe

quadrature function specified as method, insteadof the

default quad. Valid values for method are

@quadl orthe function handle of a user-defined quadrature

method that has thesame calling sequence as quad and

quadl.

Examples

Pass function handle @integrnd to dblquad:

Q = dblquad(@integrnd,pi,2*pi,0,pi);

where the function integrnd.m is:

function z = integrnd(x, y)

z = y*sin(x)+x*cos(y);

Pass anonymous function handle F to

dblquad:

F = @(x,y)y*sin(x)+x*cos(y);

Q = dblquad(F,pi,2*pi,0,pi);

The integrnd function integrates

y*sin(x)+x*cos(y) overthe square pi <=

x <= 2*pi, 0 <= y

<= pi. Note that the integrand can be evaluated

with a vector x anda scalar y.

Nonsquare regions can be handled by setting the integrand tozero

outside of the region. For example, the volume of a

hemisphereis:

dblquad(@(x,y)sqrt(max(1-(x.^2+y.^2),0)), -1, 1, -1, 1)

or

dblquad(@(x,y)sqrt(1-(x.^2+y.^2)).*(x.^2+y.^2<=1), -1, 1, -1, 1)

See Also

...............

quad2d-

Numerically evaluate double

integral over planar region

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

theintegral of fun(x,y) over the planar region

and . fun is afunction handle,

c and d mayeach be a scalar or a function

handle.

All input functions must be vectorized. The function

Z=fun(X,Y) mustaccept 2-D matrices X and

Y ofthe same size and return a matrix Z of

correspondingvalues. The functions ymin=c(X) and

ymax=d(X) mustaccept matrices and return matrices of the

same size with correspondingvalues.

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

isan 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,...)

performsthe 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 controlwhen

|Q| is sufficiently small and relative errorcontrol when

|Q| is larger. A default tolerancevalue is used when a

tolerance is not specified. The default valueof AbsTol is

1e-5. The default value of RelTol is

100*eps(class(Q)).This is also the minimum value of

RelTol. Smaller RelTol valuesare automatically

increased to the default value.

MaxFunEvals

Maximum allowed number of evaluations of

fun reached.

The MaxFunEvals parameter limits the numberof

vectorized calls to fun. The default is 2000.

FailurePlot

Generate a plot if MaxFunEvals is

reached.

Setting FailurePlot to true generatesa

graphical representation of the regions needing further

refinementwhen MaxFunEvals is reached. No plot is

generatedif 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 defaultis false.

Singular

Problem may have boundary singularities

With Singular set to true, quad2d

will employ transformations to weaken boundary singularities for

better performance. The defaultis true. Setting

Singular to false willturn these transformations

off, which may provide a performance benefiton some smooth

problems.

Examples

Example 1

Integrate over , . The true value of the integralis

.

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

Example 2

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)

Limitations

quad2d

begins by mappingthe region of integration to a rectangle.

Consequently, it may havetrouble integrating over a region that

does not have four sides orhas a side that cannot be mapped

smoothly to a straight line. Ifthe integration is unsuccessful,

some helpful tactics are leaving Singular set to its

default value of true, changing betweenCartesian and polar

coordinates, or breaking the region of integrationinto pieces and

adding the results of integration over the pieces.

For example:

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 twoareas of difficulty, near the points

(-1,0) and (1,0) andnear the circle

:

Changing the value of Singular to true

willcope with the geometric singularities at (-1,0) and

(1,0).The larger shaded areas may need refinement but are

probably not areasof 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 youcan take advantage of symmetry:

Q = 4*quad2d(fun,0,1,0,d,'Abstol',1e-8,...

'Singular',true, 'FailurePlot',true)

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

splittingthe 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 Quadraturein

2D."Applied Mathematics and Computation. Vol.

202,Issue 1, 2008, pp. 266–274.

See Also