MATLAB学习之积分（三）_matlab三维区域积分

积分

1.定积分与广义积分

定积分定义：设函数f(x) 在区间[a,b]上连续，将区间[a,b]分成n个子区间[x0,x1], (x1,x2], (x2,x3], …, (xn-1,xn]，其中x0=a，xn=b。可知各区间的长度依次是：△x1=x1-x0，在每个子区间(xi-1,xi]中任取一点ξi（1,2,…,n），作和式
$$\sum _{i=1}^{n}f({\xi }_i)\delta x_i$$
。该和式叫做积分和，设λ=max{△x1, △x2, …, △xn}（即λ是最大的区间长度），如果当λ→0时，积分和的极限存在，则这个极限叫做函数f(x) 在区间[a,b]的定积分，记为
$${\int }_a^{b}f(x)dx$$
并称函数f(x)在区间[a,b]上可积。 [1] 其中：a叫做积分下限，b叫做积分上限，区间[a, b]叫做积分区间，函数f(x)叫做被积函数，x叫做积分变量，f(x)dx 叫做被积表达式，∫ 叫做积分号。

之所以称其为定积分，是因为它积分后得出的值是确定的，是一个常数， 而不是一个函数。

**广义积分定义：**反常积分又叫广义积分，是对普通定积分的推广，指含有无穷上限/下限，或者被积函数含有瑕点的积分，前者称为无穷限广义积分，后者称为瑕积分（又称无界函数的反常积分）。

MATLAB中使用int命令可以求积分值，调用格式如下：

syms x;
f = x^3/(exp(x)-1);
int(f,0,inf)

ans =

pi^4/15
1
2
3
4
5
6
7

求广义积分只需要将积分区间相应的改为无穷即可。

2.不定积分

MATLAB中，不定积分也可以用int函数求解，调用格式如下：

syms x y;
f = (sin(x)+cos(y))/x;
int(f,0,1)

ans =

piecewise(in(y/pi - 1/2, 'integer'), sinint(1), ~in(y/pi - 1/2, 'integer'), Inf*cos(y) + sinint(1))
1
2
3
4
5
6
7

3.二重积分

MATLAB中用来进行二重积分数值计算命令的是dblquad,其调用格式如下：

进行二重积分计算时，要注意积分区间的范围，比如计算${\int }_0^{\pi }{\int }_{\pi }^{2\pi }(y\sin x+x\cos y)\mathrm{d}x\mathrm{d}y$

先建立函数M文件my2int.m

function z=my2int(x,y)
global k;
k = k+1;   %定义一个全局变量，计算被积函数调用积分命令的次数
z = y*sin(x)+x*cos(y);
1
2
3
4

然后在命令行窗口输入命令进行计算

dblquad(@my2int,pi,2*pi,0,pi)

ans =

   -9.8696

k

k =

   493
1
2
3
4
5
6
7
8
9
10
11

如果利用int进行二重积分的计算，需要注意积分区间的求取，求得积分区间之后才能进行二重积分计算。

4.三重积分

计算三重积分由于积分域的原因，计算会比二重积分计算更加复杂。

下面通过计算椭球体的积分进行说明：

计算${\iiint }_V\left(x^2+y^2+z^2\right)\mathrm{d}x\mathrm{d}y\mathrm{d}z$,其中V是由椭球体$x^2+\frac{y^2}{4}+\frac{z^2}{9}=1$围成的内部区域。

x = -1:2/50:1;
y = -2:4/50:2;
z = -3;
for i=1:51
    for j = 1:51
        z(j,i) = (1-x(i)^2-y(j)^2/4)^0.5;
        if imag(z(j,i))<0
            z(j,i) = nan;
        end
        if  imag(z(j,i))>0
            z(j,i) = nan;
        end
    end
end
mesh(x,y,z)
hold on
mesh(x,y,-z)
mesh(x,-y,z)
mesh(x,-y,-z)
mesh(-x,y,z)
mesh(-x,y,-z)
mesh(-x,-y,-z)
mesh(-x,-y,z)
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23

三维积分区域如图：

确定积分限如下：

elloposd2int
view(0,90)
title('沿x轴侧视')
view(90,0)
title('沿y轴侧视')
1
2
3
4
5

由椭球面的性质可得：

syms x y z
f = x^2 + y^2 + z^2; %被积函数
a1 = -sqrt(1-(x^2));  %积分区域
a2 = sqrt(1-(x^2));
b1 = -3*sqrt(1-x^2-(y/2)^2);
b2 = 3*sqrt(1-x^2-(y/2)^2);
fdz = int(f,z,b1,b2);
fdzdy = int(fdz,y,a1,a2);
fdzdydx = int(fdzdy,1,2);
simplify(fdzdydx)

ans =

int(int(((- 4*x^2 - y^2 + 4)^(1/2)*(- 12*x^2 + (3*y^2)/2 + 18))/2, y, -(1 - x^2)^(1/2), (1 - x^2)^(1/2)), x, 1, 2)
1
2
3
4
5
6
7
8
9
10
11
12
13
14