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

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

三角区域变换

设有三角形$△ABC$其中$A:(x_1,y_1),B:(x_2,y_2),C(x_3,y_3)$，现在我们将三角形以点$A$为定点，旋转三角形$△ABC$一定角度$\theta$，记旋转后的三角形为$△A^{\prime }{B}^{\prime }{C}^{\prime }$，使得边$A^{\prime }{B}^{\prime }$与 $y$ 轴平行（即$x_1^{\prime }=x_2^{\prime }$）且$y_1^{\prime }<y_2^{\prime }$，其中$△ABC\cong △A^{\prime }{B}^{\prime }{C}^{\prime },x_1=x_1^{\prime },y_1=y_1^{\prime }$。

• 目的：使得变换后的三角形从纵向上看，边$A^{\ast }{C}^{\ast }$为下边界，$B^{\ast }{C}^{\ast }$为上边界。

旋转变换：

设点$P:(x,y)$围绕点$A(x_1,y_2)$旋转角度$\theta$后得到$P^{\prime }:(x^{\prime },y^{\prime })$，则有

$\left[\begin{array}{c}{x}^{\prime }\\ y^{\prime }\end{array}\right]=\left[\begin{array}{cc}\cos \theta & -\sin \theta \\ \sin \theta & \cos \theta \end{array}\right]\left[\begin{array}{c}x-x_1\\ y-y_1\end{array}\right]+\left[\begin{array}{c}{x}_1\\ y_1\end{array}\right]$

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

$\left[\begin{array}{c}x\\ y\end{array}\right]=\left[\begin{array}{cc}\cos \theta & \sin \theta \\ -\sin \theta & \cos \theta \end{array}\right]\left[\begin{array}{c}{x}^{\prime }-x_1\\ y^{\prime }-y_1\end{array}\right]+\left[\begin{array}{c}{x}_1\\ y_1\end{array}\right]$

进而旋转三角形是简单的，只需以$A$点为定点，角度$\theta$由向量$\overrightarrow{AB}=(x_2-x_1,y_2-y_1)$与$y$的正半轴的夹角确定，$C$围绕$A$旋转$\theta$角度即可。

1. 当$x_2-x_1>0$时，$\alpha =\arctan \frac{y_2-y_1}{x_2-x_1}$是向量$\overrightarrow{AB}$与$x$正半轴的夹角，此时$\overrightarrow{AB}$需要旋转$\theta =\pi /2-\alpha$方能满足条件(条件：$x_1^{\prime }=x_2^{\prime }$且$y_1^{\prime }<y_2^{\prime }$)；
2. 当$x_2-x_1<0$时，$\alpha =\arctan \frac{y_2-y_1}{x_2-x_1}$是向量$\overrightarrow{AB}$与$x$负半轴的夹角，此时$\overrightarrow{AB}$需要旋转$\theta =-\pi /2-\alpha$方能满足条件(条件：$x_1^{\prime }=x_2^{\prime }$且$y_1^{\prime }<y_2^{\prime }$)；
3. 当$x_2-x_1=0$时，(1)若$y_1<y_2$，令$\theta =0$ ，即此时满足条件，无需选择；(2)若$y_1>y_2$，令$\theta =pi$ ，即将三角形以$A$为定点旋转角度$\pi$即可满足条件(条件：$x_1^{\prime }=x_2^{\prime }$且$y_1^{\prime }<y_2^{\prime }$)。

实验例子1

绘制$\overrightarrow{AB}=(x_2-x_1,y_2-y_1)$的四种情况的三角形以及旋转后的三角形

$x_2-x_1>0,y_2-y_1>0$	$x_2-x_1<0,y_2-y_1>0$

$x_2-x_1<0,y_2-y_1<0$	$x_2-x_1>0,y_2-y_1<0$

实验例子2

绘制$\overrightarrow{AB}=(x_2-x_1,y_2-y_1)$的两种一般情况的三角形以及旋转后的三角形

$x_2-x_1=0,y_2-y_1>0$	$x_2-x_1=0,y_2-y_1<0$

• 从图中易得对于变换后的三角形能够契合上述目的，即边$A^{\prime }{C}^{\prime }$为下边界，$B^{\prime }{C}^{\prime }$为上边界。

实例代码

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)>0
    sita = pi/2-atan(Y(1)/X(1));
elseif X(1)<0
    sita = -pi/2-atan(Y(1)/X(1));
else
    if Y(1)>0
        sita = 0;
    else
        sita = 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:3
    text(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)$ 在平面区域 $x_{min}\le x\le x_{max}$ 和 $y_{min(x)} ≤ y ≤ y_{max(x)} $上逼近函数 $z=fun(x,y)$的积分。

记旋转后的三角形为$△A^{\prime }{B}^{\prime }{C}^{\prime }$其中$A:(x_1^{\prime },y_1^{\prime }),B:(x_2^{\prime },y_2^{\prime }),C(x_3^{\prime },y_3^{\prime })$，设边$A^{\prime }{C}^{\prime }$所在直线方程为$y_1(x^{\prime })=a_1{x}^{\prime }+b_1$，$B^{\prime }{C}^{\prime }$所在直线方程为$y_2(x^{\prime })=a_2{x}^{\prime }+b_2$，根据上述三角形的变换，可以将变换后的三角区域书写为Ω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′)}

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

${\iint }_{\Omega }f(x,y)dxdy={\iint }_{{\Omega }_1}f(x(x^{\prime },y^{\prime }),y(x^{\prime },y^{\prime }))\left|\begin{array}{cc}\partial x/\partial x^{\prime }& \partial x/\partial y^{\prime }\\ \partial y/\partial x^{\prime }& \partial y/\partial y^{\prime }\end{array}\right|d{x}^{\prime }d{y}^{\prime }$

其中

$\left[\begin{array}{c}x\\ y\end{array}\right]=\left[\begin{array}{cc}\cos \theta & \sin \theta \\ -\sin \theta & \cos \theta \end{array}\right]\left[\begin{array}{c}{x}^{\prime }-x_1\\ y^{\prime }-y_1\end{array}\right]+\left[\begin{array}{c}{x}_1\\ y_1\end{array}\right]$，

因此

$\left|\begin{array}{cc}\partial x/\partial x^{\prime }& \partial x/\partial y^{\prime }\\ \partial y/\partial x^{\prime }& \partial y/\partial y^{\prime }\end{array}\right|=1$

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

$\{\begin{array}{c}xmin=min(x_1^{\prime },x_2^{\prime },x_3^{\prime })\\ xmax=max(x_1^{\prime },x_2^{\prime },x_3^{\prime })\\ ymin=y_1(x^{\prime })\\ ymax=y_2(x^{\prime })\\ fun(x,y)=fun(x(x^{\prime },y^{\prime }),y(x^{\prime },y^{\prime }))\end{array}$

实验例子：

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

实例代码

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)>0
    sita = pi/2-atan(Y(1)/X(1));
elseif X(1)<0
    sita = -pi/2-atan(Y(1)/X(1));
else
    if Y(1)>0
        sita = 0;
    else
        sita = 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)

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

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