函数在任意三角区域二重积分的计算

三角区域变换

设有三角形△ABC\triangle ABC△ABC其中A:(x1,y1),B:(x2,y2),C(x3,y3)A:(x_1,y_1),B:(x_2,y_2),C(x_3,y_3)A:(x1,y1),B:(x2,y2),C(x3,y3)，现在我们将三角形以点AAA为定点，旋转三角形△ABC\triangle ABC△ABC一定角度θ\thetaθ，记旋转后的三角形为△A′B′C′\triangle A'B'C'△A′B′C′，使得边A′B′A'B'A′B′与 yyy 轴平行（即x1′=x2′x_1'=x_2'x1′=x2′）且y1′<y2′y_1'<y_2'y1′<y2′，其中△ABC≅△A′B′C′,x1=x1′,y1=y1′\triangle ABC\cong\triangle A'B'C',x_1=x_1',y_1=y_1'△ABC≅△A′B′C′,x1=x1′,y1=y1′。

目的：使得变换后的三角形从纵向上看，边A∗C∗A^*C^*A∗C∗为下边界，B∗C∗B^*C^*B∗C∗为上边界。

旋转变换：

设点P:(x,y)P:(x,y)P:(x,y)围绕点A(x1,y2)A(x_1,y_2)A(x1,y2)旋转角度θ\thetaθ后得到P′:(x′,y′)P':(x',y')P′:(x′,y′)，则有

[x′y′]=[cos⁡θ−sin⁡θsin⁡θcos⁡θ][x−x1y−y1]+[x1y1]\begin{bmatrix}x' \\y' \end{bmatrix}=\begin{bmatrix}\cos\theta&-\sin\theta \\\sin\theta &\cos\theta\end{bmatrix}\begin{bmatrix}x-x_1 \\y-y_1 \end{bmatrix}+\begin{bmatrix}x_1 \\y_1 \end{bmatrix}[x′y′]=[cosθsinθ−sinθcosθ][x−x1y−y1]+[x1y1]

对应的逆变换为（这个后面求JacobiJacobiJacobi行列式时会使用到，其实这个变换的JacobiJacobiJacobi行列式为111）：

[xy]=[cos⁡θsin⁡θ−sin⁡θcos⁡θ][x′−x1y′−y1]+[x1y1]\begin{bmatrix}x \\y \end{bmatrix}=\begin{bmatrix}\cos\theta&\sin\theta \\-\sin\theta &\cos\theta\end{bmatrix}\begin{bmatrix}x'-x_1 \\y'-y_1 \end{bmatrix}+\begin{bmatrix}x_1 \\y_1 \end{bmatrix}[xy]=[cosθ−sinθsinθcosθ][x′−x1y′−y1]+[x1y1]

进而旋转三角形是简单的，只需以AAA点为定点，角度θ\thetaθ由向量AB⃗=(x2−x1,y2−y1)\vec{AB}=(x_2-x_1,y_2-y_1)AB=(x2−x1,y2−y1)与yyy的正半轴的夹角确定，CCC围绕AAA旋转θ\thetaθ角度即可。

当x2−x1>0x_2-x_1>0x2−x1>0时，α=arctan⁡y2−y1x2−x1\alpha = \arctan\frac{y_2-y_1}{x_2-x_1}α=arctanx2−x1y2−y1是向量AB⃗\vec{AB}AB与xxx正半轴的夹角，此时AB⃗\vec{AB}AB需要旋转θ=π/2−α\theta=\pi/2-\alphaθ=π/2−α方能满足条件(条件：x1′=x2′x_1'=x_2'x1′=x2′且y1′<y2′y_1'<y_2'y1′<y2′)；
当x2−x1<0x_2-x_1<0x2−x1<0时，α=arctan⁡y2−y1x2−x1\alpha = \arctan\frac{y_2-y_1}{x_2-x_1}α=arctanx2−x1y2−y1是向量AB⃗\vec{AB}AB与xxx负半轴的夹角，此时AB⃗\vec{AB}AB需要旋转θ=−π/2−α\theta=-\pi/2-\alphaθ=−π/2−α方能满足条件(条件：x1′=x2′x_1'=x_2'x1′=x2′且y1′<y2′y_1'<y_2'y1′<y2′)；
当x2−x1=0x_2-x_1=0x2−x1=0时，(1)若y1<y2y_1<y_2y1<y2，令θ=0\theta=0θ=0 ，即此时满足条件，无需选择；(2)若y1>y2y_1>y_2y1>y2，令θ=pi\theta=piθ=pi ，即将三角形以AAA为定点旋转角度π\piπ即可满足条件(条件：x1′=x2′x_1'=x_2'x1′=x2′且y1′<y2′y_1'<y_2'y1′<y2′)。

实验例子1

绘制AB⃗=(x2−x1,y2−y1)\vec{AB}=(x_2-x_1,y_2-y_1)AB=(x2−x1,y2−y1)的四种情况的三角形以及旋转后的三角形

x2−x1>0,y2−y1>0x_2-x_1>0,y_2-y_1>0x2−x1>0,y2−y1>0	x2−x1<0,y2−y1>0x_2-x_1<0,y_2-y_1>0x2−x1<0,y2−y1>0

x2−x1<0,y2−y1<0x_2-x_1<0,y_2-y_1<0x2−x1<0,y2−y1<0	x2−x1>0,y2−y1<0x_2-x_1>0,y_2-y_1<0x2−x1>0,y2−y1<0

实验例子2

绘制AB⃗=(x2−x1,y2−y1)\vec{AB}=(x_2-x_1,y_2-y_1)AB=(x2−x1,y2−y1)的两种一般情况的三角形以及旋转后的三角形

x2−x1=0,y2−y1>0x_2-x_1=0,y_2-y_1>0x2−x1=0,y2−y1>0	x2−x1=0,y2−y1<0x_2-x_1=0,y_2-y_1<0x2−x1=0,y2−y1<0

从图中易得对于变换后的三角形能够契合上述目的，即边A′C′A'C'A′C′为下边界，B′C′B'C'B′C′为上边界。

实例代码

clc,clear
% x = [1 3 4 1];
% y = [1 2 -1 1];
% x = [1 -3 -4 1];
% y = [1 2 -1 1];
% x = [1 -2 -4 1];
% y = [1 -4 3 1];
% x = [1 3 -4 1];
% y = [1 -2 1 1];
% x = [1 1 3 1];
% y = [1 2 2 1];
x = [1,1,3,1];
y = [1,-2,2,1];
X = x(2:end)-x(1);
Y = y(2:end)-y(1);
if X(1)>0sita = pi/2-atan(Y(1)/X(1));
elseif X(1)<0sita = -pi/2-atan(Y(1)/X(1));
elseif Y(1)>0sita = 0;elsesita = pi;end
end
A = [cos(sita) -sin(sita)sin(sita) cos(sita)];
P =A*([x;y]-[x(1);y(1)])+[x(1);y(1)];
plot(x,y)
hold on
plot(P(1,:),P(2,:))
tag1 = {'A','B','C'};
tag2 = ["A'","B'","C'"];
for i = 1:3text(x(i),y(i),tag1{i},...'HorizontalAlignment','right',...'Color','black')text(P(1,i),P(2,i),tag2{i},...'HorizontalAlignment','left',...'Color','red')
end
axis equal

三角区域的二重积分

使用到的内置函数：integral2

函数基本用法及所需参数：I=integral2(fun,xmin,xmax,ymin,ymax)I = integral2(fun,xmin,xmax,ymin,ymax)I=integral2(fun,xmin,xmax,ymin,ymax) 在平面区域 xmin≤x≤xmaxx_{min} ≤ x ≤ x_{max}xmin≤x≤xmax 和 $y_{min(x)} ≤ y ≤ y_{max(x)} $上逼近函数 z=fun(x,y)z = fun(x,y)z=fun(x,y)的积分。

记旋转后的三角形为△A′B′C′\triangle A'B'C'△A′B′C′其中A:(x1′,y1′),B:(x2′,y2′),C(x3′,y3′)A:(x_1',y_1'),B:(x_2',y_2'),C(x_3',y_3')A:(x1′,y1′),B:(x2′,y2′),C(x3′,y3′)，设边A′C′A'C'A′C′所在直线方程为y1(x′)=a1x′+b1y_1(x')=a_1x'+b_1y1(x′)=a1x′+b1，B′C′B'C'B′C′所在直线方程为y2(x′)=a2x′+b2y_2(x')=a_2x'+b_2y2(x′)=a2x′+b2，根据上述三角形的变换，可以将变换后的三角区域书写为Ω1={(x′,y′)∣min(x1′,x2′,x3′)≤x′≤max(x1′,x2′,x3′),y1(x′)≤y′≤y2(x′)}\Omega_1=\{(x',y')|min(x_1',x_2',x_3')≤x' ≤ max(x_1',x_2',x_3'),y_1(x') ≤ y' ≤ y_2(x')\}Ω1={(x′,y′)∣min(x1′,x2′,x3′)≤x′≤max(x1′,x2′,x3′),y1(x′)≤y′≤y2(x′)}

根据积分的变量替换规则有

∬Ωf(x,y)dxdy=∬Ω1f(x(x′,y′),y(x′,y′))∣∂x/∂x′∂x/∂y′∂y/∂x′∂y/∂y′∣dx′dy′\iint_{\Omega}f(x,y)dxdy=\iint_{\Omega_1}f(x(x',y'),y(x',y'))\begin{vmatrix} \partial x/\partial x'& \partial x/\partial y'\\ \partial y/\partial x'& \partial y/\partial y' \end{vmatrix}dx'dy'∬Ωf(x,y)dxdy=∬Ω1f(x(x′,y′),y(x′,y′))∣∣∣∣∂x/∂x′∂y/∂x′∂x/∂y′∂y/∂y′∣∣∣∣dx′dy′

其中

因此

∣∂x/∂x′∂x/∂y′∂y/∂x′∂y/∂y′∣=1\begin{vmatrix} \partial x/\partial x'& \partial x/\partial y'\\ \partial y/\partial x'& \partial y/\partial y' \end{vmatrix}=1∣∣∣∣∂x/∂x′∂y/∂x′∂x/∂y′∂y/∂y′∣∣∣∣=1

所以三角区域上f(x,y)f(x,y)f(x,y)二重积分为I=integral2(fun,xmin,xmax,ymin,ymax)I = integral2(fun,xmin,xmax,ymin,ymax)I=integral2(fun,xmin,xmax,ymin,ymax) ，其中：

{xmin=min(x1′,x2′,x3′)xmax=max(x1′,x2′,x3′)ymin=y1(x′)ymax=y2(x′)fun(x,y)=fun(x(x′,y′),y(x′,y′))\left\{\begin{matrix} xmin= min(x_1',x_2',x_3')\\ xmax = max(x_1',x_2',x_3')\\ ymin = y_1(x')\\ymax = y_2(x')\\fun(x,y)=fun(x(x',y'),y(x',y'))\end{matrix}\right.⎩⎪⎪⎪⎪⎨⎪⎪⎪⎪⎧xmin=min(x1′,x2′,x3′)xmax=max(x1′,x2′,x3′)ymin=y1(x′)ymax=y2(x′)fun(x,y)=fun(x(x′,y′),y(x′,y′))

实验例子：

这里我们就随意举个例子，算三角形的面积，即fun(x,y)=1fun(x,y)=1fun(x,y)=1，三角形的三个顶点为A:(1,1),B:(3,4),C(5,2)A:(1,1),B:(3,4),C(5,2)A:(1,1),B:(3,4),C(5,2)，该三角形的面积为555，下面是代码跑出来的结果：

实例代码

clc,clear
% 该函数用于求三角区域上函数的积分
point = [1,3,5;1,4,2];
x = point(1,:);
y = point(2,:);
X = x(2:3) - x(1);
Y = y(2:3) - y(1);
if X(1)>0sita = pi/2-atan(Y(1)/X(1));
elseif X(1)<0sita = -pi/2-atan(Y(1)/X(1));
elseif Y(1)>0sita = 0;elsesita = pi;end
end
A = [cos(sita) -sin(sita)sin(sita) cos(sita)];
P =A*([x;y]-[x(1);y(1)])+[x(1);y(1)];
xp = P(1,:);
yp = P(2,:);
f = @(t,s) t-t+s-s+1;
% 本来上面的f应该写为 f = @(t,s) 1的，但是为了保证输入和输出具有相同的维度，我们做了如上的处理
% 不然在调用integral2函数时会报输入与输出大小不匹配的错误，导致程序崩溃。
fun = @(t,s) f(cos(sita)*(t-x(1))+sin(sita)*(s-y(1))+x(1),-sin(sita)*(t-x(1))+cos(sita)*(s-y(1))+y(1));
y_down = @(t) (yp(3)-yp(1))/(xp(3)-xp(1)).*(t-xp(1)) + yp(1);
y_up = @(t) (yp(3)-yp(2))/(xp(3)-xp(2)).*(t-xp(2)) + yp(2);
I = integral2(fun,min(xp),max(xp),y_down,y_up);
disp('积分值')
disp(I)

最后，附上代码中关于注释所说的错误，我们仅仅将fff写为f=@(t,s)1f=@(t,s) 1f=@(t,s)1;

有兴趣的小伙伴可以自己修改了观察。当然，如果你给出的f(x,y)f(x,y)f(x,y)不是常数，也不是单变量函数(即f(x,y)=f(x)f(x,y) = f(x)f(x,y)=f(x)或者f(x,y)=f(y)f(x,y)=f(y)f(x,y)=f(y)的情形)，那么是不存在上述的错误类型的。