第10章重积分

重积分针对多元函数，比较典型的是二重积分和三重积分。二重积分是二元函数在空间上的积分，同定积分类似，是某种特定形式的和的极限，本质是求曲顶柱体体积。重积分有着广泛的应用，可以用来计算曲面的面积、平面薄片重心等。三重积分的几何意义和物理意义是不均匀的空间物体的质量。知道重积分的这些意义，有助于理解重积分的概念和应用场景，也就能明确学习重积分的意义了。
10.1本章目标
本章将用MATLAB实现以下操作： 
(1) 二重积分的数值计算。
(2) 极坐标下二重积分的数值计算。
(3) 三重积分的数值计算。
(4) 柱面坐标系及球面坐标系下三重积分的数值计算。
10.2相关命令
本章涉及的MATLAB命令如下。
(1) integral2： 对二重积分进行数值计算。用法如下： 

 q=integral2(fun,xmin,xmax,ymin,ymax)： 在平面区域xmin≤x ≤xmax和ymin(x)≤y≤ymax(x)上逼近函数z=fun(x,y)的积分。

 q=integral2(fun,xmin,xmax,ymin,ymax,Name,Value)： 指定具有一个或多个Name,Value对组参数的其他选项。
(2) integral3： 对三重积分进行数值计算。用法如下： 

 q=integral3(fun,xmin,xmax,ymin,ymax,zmin,zmax)： 在区域xmin≤x≤xmax、ymin(x)≤y≤ymax(x)和zmin(x,y)≤z≤zmax(x,y)逼近函数z=fun(x,y,z)的积分。

 q=integral3(fun,xmin,xmax,ymin,ymax,zmin,zmax,Name,Value)： 指定具有一个或多个Name,Value对组参数的其他选项。




10.3二重积分的计算
设函数f(x,y)是有界闭区域D上的有界函数。将闭区域D任意分成n个小闭区域Δσ1,Δσ2,…,Δσn，其中Δσi表示第i个小闭区域。在每个Δσi上任取一点(ξi,ηi)，作乘积f(ξi,ηi)Δσi(i=1,2,…,n)，并作和∑ni=1f(ξi,ηi)Δσi，如果当各小闭区域的直径中的最大值λ趋于零时，这个和的极限总存在，那么称此极限为函数f(x,y)在闭区域D上的二重积分，记作Df(x,y)dσ，即

Df(x,y)dσ=limλ→0∑ni=1fξi,ηiΔσi

其中，f(x,y)叫作被积函数，f(x,y)dσ叫作被积表达式，dσ叫作面积元素，x与y叫作积分变量，D叫作积分区域，∑ni=1f(ξi,ηi)Δσi叫作积分和。
在二重积分的定义中对闭区域D的划分是任意的，如果在直角坐标系中用平行于坐标轴的直线网来划分D，那么除了包含边界点的一些小闭区域外，其余的小闭区域都是矩形闭区域。设矩形区域Δσi的边长为Δxj和Δyk，则Δσi=Δxj·Δyk。因此在直角坐标系中，有时也把面积元素dσ记作dxdy，而把二重积分记作

Df(x,y)dxdy

其中，dxdy叫作直角坐标系中的面积元素。
根据上面二重积分直角坐标系中的定义，可以编写近似的用定义方法求解的MATLAB函数，具体方法是先获取积分区域D在x和y方向上的最大值、最小值，从而将积分区域扩展成矩形区域［a,b］×［c,d］，判断区域D在矩形区域［a,b］×［c,d］中的部分，通过上述二重积分定义式求解。按照定义来计算二重积分对于少部分简单的积分区域和被积函数来说是可行的，因为使用情况较少，具体的计算方法函数代码就不在此呈现。在MATLAB中有专门的函数来更方便地实现二重积分。
10.3.1二重积分的数值计算
首先考虑特殊的矩形区域的二重积分。
计算下面的二重积分问题

∫ba∫dcf(x,y)dxdy

矩形区域的二重积分的数值解可以直接使用MATLAB提供的integral2()函数直接求出。
可以通过以下命令调出integral2()函数的代码了解具体的计算方法： 

edit(which('integral2.m'))

这里更关注如何使用该函数来实现二重积分的求解。
例101试求出二重积分

I=∫1-1∫2-2e-x22sin(x2+y)dxdy

解： 



f=@(x,y)exp(-x.^2/2.*sin(x.^2+y));

y=integral2(f,-2,2,-1,1)



运行程序，得到如下结果： 

y = 8.1235

同时可知，integral2（）也可以用于非矩形区域二重积分的计算。
例102计算函数f(x,y)=1x+y(1+x+y)2在0≤x≤1，0≤y≤1-x所围成的区域上的积分。
解： 
创建匿名函数： 



fun = @(x,y) 1./( sqrt(x + y) .* (1 + x + y).^2 );



对0≤x≤1，0≤y≤1-x限定的三角形区域计算积分： 



ymax = @(x) 1 - x;

q = integral2(fun,0,1,0,ymax)



运行程序，得到如下结果： 

q = 0.2854

10.3.2直角坐标计算
为了解决更一般的重积分计算问题，接下来使用int函数，将二重积分化为二次积分来计算。依据高等数学中的知识有等式

Df(x,y)dσ=∫ba∫φ2(x)φ1(x)f(x,y)dydx

上式右端的积分是先对y、后对x的二次积分，也常记作

∫badx∫φ2(x)φ1(x)f(x,y)dy

这就是把二重积分化为先对y、后对x的二次积分的公式。这样的积分区域是X型的，可用φ1(x)≤y≤φ2(x)，a≤x≤b表示。若积分区域是Y型的则用公式

Df(x,y)dσ=∫dcdy∫ψ2(y)ψ1(y)f(x,y)dx

来计算。如果积分区域既是X型的又是Y型的，则两个不同次序的二次积分相等。
将二重积分化为二次积分时，确定积分限是一个关键。积分限是根据积分区域D来确定的，先画出积分区域D的图形。假如积分区域是X型的，在区间［a,b］上任意取定一个x值，积分区域上以这个x值为横坐标的点在一段直线上，这段直线平行于y轴，该线段上点的纵坐标就是先把x看作常量而对y积分时的下限和上限。再把x看作变量而对x积分时，积分区间就是［a,b］。
例103计算Dxydσ，其中D是由直线y=1、x=2及y=x所围成的闭区域。
解： 
第一步，绘制积分区域： 



y1=1;

fplot(y1);

hold on;

y2=@(x)x;

fplot(y2);

hold on;

fimplicit(@(x,y) x-2);

fill（［1,2,2,1］,［1,2,1,1］,'r'）;

xlabel('x');

ylabel('y');



运行程序，得到如图101所示的积分区域。


图101由直线y=1、x=2及y=x所围成的积分区域


第二步，确定积分上、下限： 



fun=@(x)x-1;

A=fzero(fun,0);

B=2;



第三步，计算积分值： 



syms x y

I=int(int(x*y,y,y1,y2),x,A,B)



运行程序，得到如下结果： 


I = 98

在化二重积分为二次积分时，为了计算方便，需要选择恰当的二次积分的顺序。这时既要考虑积分区域的形状，又要考虑被积函数的特性。
例104计算Dxydσ，其中D是由抛物线y2=x及直线y=x-2所围成的闭区域。
解： 
第一步，绘制积分区域： 



syms x y;

fimplicit(y-x+2);

hold on;

fimplicit(y^2-x);

xlabel('x');

ylabel('y');



运行程序，得到如图102所示的积分区域。


图102由抛物线y2=x及直线y=x-2所围成的积分区域


D既是X型的又是Y型的。若利用X型的公式计算，则由于在区间［0,1］及［1,4］上表示φ1(x)的式子不同，所以要用经过交点(1,-1)且平行于y轴的直线x=1把区域D分成D1和D2两部分，其中

D1=(x,y)-x≤y≤x,0≤x≤1
D2=(x,y)x-2≤y≤x,1≤x≤4

而如果D是Y型的，不用将区域分成两部分计算，是较为简便的。
第二步，确定积分上、下限： 



A=double(solve((x-2)^2-x,0));



第三步，计算积分值： 



y2=@(x)x-2;

x1=@(y)y^2;

x2=@(y)y+2;

I=int(int(x*y,x,x1,x2),y,y2(A(1)),y2(A(2)))



运行程序，得到如下结果： 

I = 458

10.3.3极坐标计算
当二重积分的积分区域D的边界曲线用极坐标方程来表示比较方便，且某些被积函数用极坐标变量ρ、θ表达比较简单时，就可以考虑用极坐标来计算二重积分Df(x,y)dσ。

Df(x,y)dσ=Df(ρcosθ,ρsinθ)ρdρdθ

上式就是二重积分的变量从直角坐标变换为极坐标的变换公式，其中ρdρdθ就是极坐标系中的面积元素。
极坐标系中的二重积分，同样可以化为二次积分来计算。设积分区域D可以用不等式
φ1(x)≤ρ≤φ2(x)，α≤θ≤β

来表示，其中函数φ1(θ)、φ1(θ)≤ρ≤φ2(θ)在区间［α,β］上连续。极坐标系中的二重积分化为二次积分的公式为

Dfρcosθ,ρsinθρdρdθ=∫βαdθ∫φ2(θ)φ1(θ)
f(ρcosθ,ρsinθ)ρdρ

例105计算De-x2-y2dxdy，其中D是由圆心在原点、半径为a的圆周所围成的闭区域。
解： 
本题如果用直角坐标计算，因为积分∫e-x2dx不能用初等函数表示，所以不容易计算。如果利用极坐标计算，区域D可表示为0≤ρ≤a，0≤θ≤2π，且比较利于积分的计算，具体实现代码为： 



syms rho theta a;

x=rho*cos(theta);

y=rho*sin(theta);

f=exp(-x^2-y^2);

I=int(int(f*rho,rho,0,a),theta,0,2*pi)



运行程序，得到如下结果： 

I = -π(e-a2-1)

10.3.4二重积分换元法
将二重积分的变量从直角坐标变换为极坐标是二重积分换元法的一种特殊情形。把平面上的点同时用直角坐标(x,y)、极坐标(ρ,θ)表示，它们之间的关系是
x=ρcosθ
y=ρsinθ

也可以看成建立了ρOθ平面和xOy平面上的一一对应关系。
一般的二重积分换元法有如下定理： 设f(x,y)在xOy平面上的闭区域D上连续，若变换
T： x=xu,v，y=y(u,v)

将uOv平面上的闭区域D′变为xOy平面上的D，且满足
(1) x(u,v)，y(u,v)在D′上具有一阶连续偏导数。
(2) 在D′上雅可比式
J(u,v)=(x,y)(u,v)≠0
(3) 变换T： D′→D是一对一的，则有

Df(x,y)dxdy=Dfx(u,v),y(u,v)J(u,v)dudv

该公式称为二重积分的换元公式。
例106计算Dey-xy+xdxdy，其中D是由x轴、y轴和直线x+y=2所围成的闭区域。
解： 
令u=y-x，v=y+x，则x=v-u2，y=v+u2。
作变换x=v-u2，y=v+u2，绘制xOy平面上的闭区域D和uOv平面上的对应区域D′。
第一步，绘制积分区域： 



syms x y u v;

U=y-x;

V=y+x;

S=solve(u-U,v-V,x,y);

subplot 121;

ezplot(x+y-2,［-1,3］);

hold on;

plot(［0,0,3］,［3,0,0］);

fill(［0,0,2,0］,［2,0,0,2］,'r');

axis equal tight;

title('原积分区域图');

subplot 122;

ezplot(S.x,［-2,2］);

hold on;

ezplot(S.y,［-2,2］);

v1=solve(S.x+S.y-2,v);

ezplot(v1,［-2,2］);

fill(［0,2,-2,0］,［0,2,2,0］,'r');

axis equal tight;

title('变换后的积分区域图');



运行程序，得到如图103所示的积分区域。


图103例106的积分区域


第二步，根据变换后的积分区域计算二重积分： 



f=exp((y-x)/(y+x));

f=subs(f,{x,y},{S.x,S.y});

J=det(jacobian(［S.x;S.y］,［u,v］));

I=int(int(f*abs(J),u,-v,v),v,0,2)



运行程序，得到如下结果： 

I = e-e-1

10.4三重积分
设f(x,y,z)是空间有界闭区域Ω上的有界函数。将Ω任意分成n个小闭区域

Δv1,Δv2,…,Δvn

其中，Δvi表示第i个小闭区域，也表示它的体积。在每个Δvi上任取一点(ξi,ηi,ζi)，作乘积f(ξi,ηi,ζi)Δvi(i=1,2,…,n)，并作和∑ni=1f(ξi,ηi,ζi)Δvi，如果当各小闭区域直径中的最大值λ→0时，这个和的极限总存在，且与闭区域Ω的分法及点(ξi,ηi,ζi)的取法无关，那么称此极限为函数f(x,y,z)在闭区域Ω上的三重积分，记作Ωf(x,y,z)dv，即

Ωf(x,y,z)dv=limλ→0∑ni=1f(ξi,ηi,ζi)Δvi

其中，f(x,y,z)叫作被积函数，dv叫作体积元素，Ω叫作积分区域。
在直角坐标系中，有时也把体积元素dv记作dxdydz，而把三重积分记作

Ωf(x,y,z)dxdydz

其中，dxdydz叫作直角坐标系中的体积元素。
当函数f(x,y,z)在闭区域连续时，limλ→0∑ni=1f(ξi,ηi,ζi)Δvi必存在，也就是函数在闭区域上的三重积分必定存在。三重积分的性质与二重积分的性质相似。
10.4.1利用直角坐标计算三重积分
下面介绍利用直角坐标计算三重积分。
若平行于z轴且穿过闭区域内部的直线与闭区域Ω的边界曲面S相交不多于两点，可以把Ω投影到xOy面上。设积分区域可以表示为

Ω=(x,y,z)z1(x,y)≤z≤z2(x,y),(x,y)∈Dxy

先将x、y看作定值，将f(x,y,z)只看作z的函数，在区间［z1(x,y),z2(x,y)］上对z积分，积分的结果是x、y的函数，记为F(x,y),即

F(x,y)=∫z2(x,y)z1(x,y)f(x,y,z)dz

然后计算F(x,y)在闭区域Dxy上的二重积分

DxyF(x,y)dσ=Dxy∫z2(x,y)z1(x,y)f(x,y,z)dzdσ

设闭区域

Dxy=(x,y)y1(x)≤y≤y2(x),a≤x≤b

把这个二重积分化为二次积分，于是有三重积分的计算公式

Ωf(x,y,z)dv=∫badx∫y2(x)y1(x)dy∫z2(x,y)z1(x,y)f(x,y,z)dz

此公式把三重积分化为先对z、次对y、最后对x的三次积分。其他积分次序的公式也类似。
例107计算Ωxdxdydz，其中Ω为三个坐标面及平面x+2y+z=1所围成的闭区域。
解： 
第一步，绘制积分区域： 



［X,Y］=meshgrid(linspace(0,1,30));

mesh(X,Y,1-X-2*Y);

hold on;

mesh(X,Y,zeros(size(X)));

mesh(X,Y,zeros(size(Y)));

hidden off;

view(［60,10］);

xlabel('x');






ylabel('y');

zlabel('z');



运行程序，得到如图104所示的积分区域。


图104例107的三重积分区域


第二步，计算三重积分： 



syms x y z;

I=int(int(int(x,z,0,1-x-2*y),y,0,(1/2)*(1-x)),x,0,1)



运行程序，得到如下结果： 

I = 148

例108计算三重积分Ωy1-x2dxdydz，其中Ω由曲面y=-1-x2-z2、x2+z2=1和平面y=1围成。
解： 
第一步，绘制积分区域： 



［theta,rho］=meshgrid(linspace(0,2*pi),linspace(0,1-eps));

［X,Z］=pol2cart(theta,rho);

surf(X,-sqrt(1-X.^2-Z.^2),Z);

hold on;

［X1,Y1,Z1］=cylinder(ones(1,20),40);

surf(X1,Z1,Y1);

surf(X,ones(size(X)),Z);

shading flat;

view(［-100,30］);






axis equal;

xlabel('x');

ylabel('y');

zlabel('z');



运行程序，得到如图105所示的积分区域。


图105例108中的积分区域


容易看出,此处积分区域按y、z、x的次序三次积分。
第二步，计算三重积分： 



syms x y z;

I=int(int(int(y*sqrt(1-x^2),y,-sqrt(1-x^2-z^2),1),z,-sqrt(1-x^2),sqrt(1-x^2)),x,-1,1)



运行程序，得到如下结果： 

I = 2845

10.4.2利用柱面坐标计算三重积分
如图106所示，设M(x,y,z)为空间内一点，并设点M在xOy面上的投影P的极坐标为ρ、θ，则这样的三个数ρ、θ、z就叫作点M的柱面坐标，这里规定ρ、θ、z的变化范围为

0≤ρ＜+∞
0≤θ≤2π
-∞＜z＜+∞



图106点M及其投影P在空间中的位置关系示意图


三组坐标面分别为： 
ρ=常数，即以z轴为轴的圆柱面；
θ=常数，即过z轴的半平面；
z=常数，即与xOy面平行的平面。
转换关系为
x=ρcosθ
y=ρsinθ
z=z
能够得到柱面坐标系中的体积元素为
dv=ρdρdθdz
把三重积分的变量从直角坐标变换为柱面坐标的公式

Ωf(x,y,z)dxdydz=ΩF(ρ,θ,z)ρdρdθdz

其中，F(ρ,θ,z)=f(ρcosθ,ρsinθ,z)。
转换为柱面坐标系后的三重积分可以根据ρ、θ、z在积分区域中的变化范围确定化为三次积分的先后顺序。
例109利用柱面坐标计算三重积分Ωzdxdydz，其中Ω是由曲面z=x2+y2与平面z=4所围成的闭区域。
解： 
把闭区域投影到xOy面上，得半径为2的圆形闭区域
Dxy=(ρ,θ)0≤ρ≤2,0≤θ≤2π
在Dxy内任取一点(ρ,θ)，过此点作平行于z轴的直线，此直线通过曲面z=x2+y2穿入闭区域内，然后通过平面z=4穿出Ω外。因此闭区域Ω可用不等式
ρ2≤z≤4， 0≤ρ≤2， 0≤θ≤2π

来表示。
于是编写程序计算： 



syms rho theta z a

x=rho*cos(theta);

y=rho*sin(theta);

f=z*sqrt(x^2+y^2);

I=int(int(int(z,z,rho^2,4)*rho,rho,0,2),theta,0,2*pi)



运行程序，得到如下结果： 

I = 64π3

10.4.3利用球面坐标计算三重积分



图107点M在空间中的位置示意图


如图107所示，设M(x,y,z)为空间内的一点，则点M也可以用r、φ、θ来确定，其中r为原点O与点M之间的距离，φ为有向线段OM与z轴正向所夹的角，θ为从正z轴来看自x轴按逆时针方向转到有向线段OP的角，这里点P为点M在xOy面上的投影。
r，φ，θ叫作点M的球面坐标，变化范围是： 
0≤r<+∞
0≤φ≤π
0≤θ≤2π
三组坐标面分别为： 
r=常数，即以原点为心得到球面；
φ=常数，即以原点为顶点、z轴为轴的圆锥面；
θ=常数，即过z轴的半平面。
直角坐标与球面坐标的关系为： 

x=rsinφcosθ
y=rsinφsinθ
z=rcosφ
球面坐标系中的体积元素： 
dv=r2sinφdrdφdθ
有关系式

Ωf(x,y,z)dv=ΩF(r,φ,θ)r2sinφdrdφdθ

其中，F(r,φ,θ)=f(rsinφcosθ,rsinφsinθ,rcosφ)。这就是把三重积分的变量从直角坐标变换为球面坐标的公式。
例1010求半径为a的球面与半顶角为α的内接锥面所围成的几何体(如图108所示)的体积。


图108球面与内接锥面所围成的几何体示意图


解： 



syms rho theta phi a alpha;

x=rho*sin(phi)*cos(theta);

y=rho*sin(phi)*sin(theta);

z=rho*cos(phi);

I=int(int(int(rho^2*sin(phi),rho,0,2*a*cos(phi)),phi,0,alpha),theta,0,2*pi)



运行程序，得到如下结果： 

I = -2πa32cos4a3-23

10.5拓展内容
利用重积分计算几何图形面积或体积的例子比较多，以下将补充一些重积分应用的案例。
10.5.1重积分补充案例
例1011计算曲线f=x2与直线g=3x所围成区域的面积。
解： 



syms x y

f = x^2; g = 3*x;






figure

fplot(f, ［0 3］, 'LineWidth', 2); grid on; hold on

fplot(g, ［0 3］, 'LineWidth', 2)

xlabel('x');

ylabel('y');

axis(［0 3 0 9］); xticks(0:3); yticks(0:3:9)



运行程序，得到如图109所示的积分区域。



图109曲线f=x2与直线g=3x所围成的区域





I1 = int(int(x*y, y, x^2, 3*x), x, 0, 3)



I1 = 2438



I1_appx = double(I1)



I1_appx = 30.3750



I2 =int(int(x*y, x, y/3, sqrt(y)), y, 0, 9)



I2 = 2438



I2_appx = double(I2)



I2_appx = 30.3750
例1012计算四分之一圆盘的积分。
解： 



syms x y

figure

t = linspace(0,pi/2); xt = cos(t); yt = sin(t);

X = ［0 xt 0］; Y = ［0 yt 0］;

fill(X, Y, 'y'); grid on ; hold on

plot(X, Y, 'LineWidth', 2)

plot(［-0.5 1.5］, ［0 0］, 'k')

plot(［0 0］, ［-0.5 1.5］, 'k')

xlabel('x');

ylabel('y');

xticks(0:1); yticks(0:1)

axis equal;

axis(［-0.5 1.5 -0.5 1.5］)



运行程序，得到如图1010所示的积分区域。



图1010四分之一圆盘的区域






I = int(int(x*y, y, 0, sqrt(1-x^2)), x, 0, 1)



I= 18



I_appx = double(I)



I_appx = 0.1250
例1013计算以y=x2、z=3y、z=2+y三面为界的固体的体积。
解： 



syms x y z

f(x) = x^2;

% 这是两个平面，它们的交点是r(x)

plane1 = z == 3*y;

plane2 = z == 2+y; 

L = solve(［plane1 plane2］, ［y z］);

r(x) = ［x L.y L.z］; % 将这条直线投影到xOy平面上:y = 1, x自由

x_bounds = solve(x^2 == 1, x); % y_bounds: x^2 <= y <= 1。注意这里的y是非负的，且

% 0 <= y <= 1

% 注意平面2在平面1之上

V = int(int(2+y - 3*y, y, x^2, 1), x, -1, 1)



V =1615



V_appx = double(V)



V_appx = 1.0667



figure

fsurf(x, y, 2+y, ［-1 1 -0.25 1.25］, 'r', 'MeshDensity', 10); hold on % 平面 z = 2+y 以上

fsurf(x, y, 3*y, ［-1 1 -0.25 1.25］, 'g', 'MeshDensity', 10) % 平面 z = 3*y 以下

fsurf(x, x^2, z, ［-1.25 1.25 -1 4］, 'y', 'MeshDensity', 16)% 侧面为抛物柱面

xlabel('x');

ylabel('y');

zlabel('z');



运行程序，得到如图1011所示的几何体。



图1011y=x2、z=3y、z=2+y三个面所围成的几何体





figure

xs = linspace(-1,1); ys = xs.^2;

X = ［xs xs(1)］; Y = ［ys ys(1)］;

fill(X,Y,'y'); grid on; hold on

plot(X,Y,'b', 'LineWidth', 2)

plot(［-1 1］, ［0 0］, 'k')

plot(［0 0］, ［-0.5 1.5］, 'k')

axis equal; axis(［-1 1 -0.5 1.5］)

xlabel('x');

ylabel('y');

xticks(-1:1); yticks(0:1)



运行程序，得到如图1012所示的二维积分区域。


图1012二维积分区域



例1014确定两个椭圆抛物面与圆柱所围几何体的体积。
解： 



syms r x y z theta

LP = z == 2*x^2 + y^2; UP = z == 8 - x^2 - 2*y^2; % 上、下抛物面

f = rhs(LP); g = rhs(UP); s = sqrt(1-x^2); s83 = sqrt(8/3);

CC = x^2 + y^2 == 1; % 圆柱

V = int(int(g-f, y, -s, s), x, -1, 1)



V =13π2



V_appx = double(V)



V_appx = 20.4204



ct = cos(theta); st = sin(theta); xp = r*ct; yp = r*st;

fp = simplify(subs(f, ［x y］, ［xp yp］));

gp = simplify(subs(g, ［x y］, ［xp yp］));

figure

fsurf(xp, yp, gp, ［0 s83 0 2*pi］, 'r', 'MeshDensity', 20); hold on % 抛物面 z = 8 - x^2 - 2*
y^2 之上

fsurf(xp, yp, fp, ［0 s83 0 2*pi］, 'g', 'MeshDensity', 20) % 抛物面 z = 2*x^2 + y^2 之下

fsurf(ct, st, z, ［0 2*pi -1 8］, 'y', 'MeshDensity', 20) % 圆柱体 x^2 + y^2 = 1 的侧面

xlabel('x');

ylabel('y');

zlabel('z');

xticks(-1:1); yticks(-1:1); zticks(0:4:8)



运行程序，得到如图1013所示的几何体。


图1013两个椭圆抛物面与圆柱所围的几何体


例1015确定积分∫10∫π4arctanxf(x,y)dxdy改变积分顺序后的积分上、下限并绘制积分区域。
解： 



figure

x = linspace(0,1); x = tan(y);






X = ［x 0 0］; Y = ［y pi/4 0］;

fill(X,Y,'y'); grid on; hold on

plot(X,Y,'b', 'LineWidth', 2)

plot(［-0.5 1.5］, ［0 0］, 'k')

plot(［0 0］, ［-0.5 1.5］, 'k')

axis equal; axis(［-0.5 1.5 -0.5 1.5］)

xlabel('x');

ylabel('y');

xticks(0:1); yticks(0:1)

text(0.5, 0.4, 'y = arctan(x) or x = tan(y)', 'FontSize', 12)

text(0.4, 0.9, 'y = ＼pi/4', 'FontSize', 12)

text(-0.25, 0.5, 'x = 0', 'FontSize', 12)



运行程序，得到如图1014所示的二维积分区域。


图1014二维积分区域


例1016改变积分∫80∫2y13ex4dxdy的顺序并计算积分。
解： 



syms x y

I = int(int(exp(x^4), y, 0, x^3), x, 0, 2)



I = e164-14



I_appx = double(I)



I_appx = 2.2215e+06



figure

xs = linspace(0,2); ys = xs.^3;

X = ［xs 2 0］; Y = ［ys 0 0］;

fill(X,Y,'y'); grid on; hold on

plot(X,Y,'b', 'LineWidth', 2)

plot(［-2 4］, ［0 0］, 'k')

plot(［0 0］, ［-1 9］, 'k')

axis equal; axis(［-2 4 -1 9］)

xlabel('x');

ylabel('y');

xticks(0:2:2); yticks(0:2:8)

text(-1.25, 4.6, 'x = y^{1/3} 或 y = x^3', 'FontSize', 12)

text(0.6, -0.5, 'y = 0', 'FontSize', 12)

text(2.5, 3.8, 'x = 2', 'FontSize', 12)



运行程序，得到如图1015所示的积分区域。


图1015二维积分区域


例1017求平面x+y+z=1和曲面z=4-x2-y2所围几何体的体积。
解： 



syms r t x y z

plane = x + y + z == 1; paraboloid = z == 4 - x^2 - y^2;







eq1 = subs(paraboloid, z, 1-x-y); 

eq2 = lhs(eq1) - rhs(eq1) == 0 ; 

eq3 = eq2 + 1/4 + 1/4 + 3 ;      

c = (x-1/2)^2 + (y-1/2)^2 == 7/2;

xsp = 1/2 + r*cos(t); ysp = 1/2 + r*sin(t); %极坐标变换

f = 1 - xsp - ysp;

g = 4 - xsp^2 - ysp^2 ;

V = int(int((g-f)*r, r, 0, sqrt(7/2)), theta, 0, 2*pi)



V =49π8



V_appx = double(V) % 体积



V_appx = 19.2423



figure

fsurf(xsp, ysp, g, ［0 sqrt(7/2) 0 2*pi］, 'r', 'MeshDensity', 20); hold on % 抛物面上面

fsurf(xsp, ysp, f, ［0 sqrt(7/2) 0 2*pi］, 'g', 'MeshDensity', 10) % 平面下面

xlabel('x');

ylabel('y');

zlabel('z');

axis equal

view(-40,44)



运行程序，得到如图1016所示的积分区域。


图1016平面x+y+z=1和曲面z=4-x2-y2所围的几何体


10.5.2四维积分的计算
MATLAB中的integral求积分函数直接支持一维、二维和三维积分。然而，要求解四维和更高阶积分，需要嵌套对求解器的调用。接下来将介绍一个通过使用integral3和integral的嵌套调用来计算四维球体体积的例子。
已知半径为r的四维球体的体积为： 

V4(r)=∫2π0∫π0∫π0∫r0r3sin2(θ)sin()drdθddξ

现在需要计算该球体的数值解。
这是一个四维积分，要实现对它的求解，首先可以为被积函数创建函数句柄f(r,θ,,ξ)： 



f = @(r,theta,phi,xi) r.^3 .* sin(theta).^2 .* sin(phi);



接下来，创建一个函数句柄，它使用integral3计算三个积分： 



Q = @(r) integral3(@(theta,phi,xi) f(r,theta,phi,xi),0,pi,0,pi,0,2*pi);



最后，使用Q作为被积函数，用integral来对其进行积分。需要为半径r选择一个值，此处r=2： 



I = integral(Q,0,2,'ArrayValued',true)



I = 78.9568
得到的确切答案是π2r42Γ(2)，然后计算其数值： 



I_exact = pi^2*2^4/(2*gamma(2))



I_exact=78.9568
10.6上机实践
1. 计算I=D1-x2a2-y2b2dxdy，其中D为椭圆x2a2+y2b2=1所围成的闭区域。
2. 计算三重积分I=∫2-2dz∫3-3dy∫2-2(sinx2+z2cosy)dx的数值解并验证数值解的准确性。
3. 计算三重积分
I=Ωxyzdv
其中，Ω由曲面z=xy，x+y+z=1及z=0围成。
4. 计算三重积分

I=Ωx2+y22 dv 

其中，Ω由曲线y2=2z
x=0绕z轴旋转一周而成的曲面与平面z=8围成。