# 數值計算：三角形積分

陳橙橙發表於2022-03-21

數值計算：三角形積分

書接上回《高斯-勒朗德積分公式》

需求：在給定空間三角形\(\Delta ABC\)中，\(A(x_1,y_1,z_1),B(x_2,y_2,z_2),C(x_3,y_3,z_3)\)，已知函式\(f(x,y,z)\)，求利用數值方法求解積分：\(\iint_{\Delta ABC}f(x,y,z)\text dS\)。

解決方法：參考triangle_lyness_rule給出的積分方法，具體細節也不是太懂，但是思路上與高斯積分類似，計算平面上的積分點與係數權重進行積分

三角形積分點與積分權重計算

Triangle Llyness Rule

triangle_lyness_rule中給出了不同階數下，在標準三角形中的係數點位置與權重係數。例如下圖中展示\(Rule=10\)時的積分點位置與權重係數。下表中顯示了不同\(Rule\)下的積分精度\(Precision\)積分點數目\(order\)以及積分點是否包含三角形中心\(center\)。

image

Rule Order Precision Center
0 1 1 YES
1 3 2 NO
2 4 2 YES
3 4 3 YES
4 7 3 YES
5 6 4 NO
6 10 4 YES
7 9 4 NO
8 7 5 YES
9 10 5 YES
10 12 6 NO
11 16 6 YES
12 13 6 YES
13 13 7 YES
14 16 7 YES
15 16 8 YES
16 21 8 NO
17 16 8 YES
18 19 9 YES
19 22 9 YES
20 27 11 NO
21 28 11 YES

座標變換

image

採用與之前文章中《高斯-勒朗德積分公式》形函式方式計算座標轉換關係。得到三節點形函式為，剩下步驟與《高斯-勒朗德積分公式》中類似。

\[\begin{cases} N_1(s,t)=1-s-t\\ N_2(s,t)=s\\ N_3(s,t)=t\\ \end{cases}\\ \]

改進

在KY師兄指點下，以上步驟可以進一步簡化。原因在於三角形座標變換的形函式簡單，可以直接進行座標運算，Jacobi係數直接等於三角形面積，具體見程式碼。

積分測試

以下為測試積分函式，其中\(LYNESS_RULE.txt\)儲存的資料太長了，就放到Gitee:連結待更新倉庫了。

%% 測試三角形積分
clc;clear;
global TriCoeff
% 匯入積分系數
TriCoeff=loadLynessFromTxT("LYNESS_RULE.txt");

P1=[0,0,0];
P2=[2,0,0];
P3=[0,3,0];
% 積分函式
func=@(x,y,z) (x^6+y^3+1);

count=1;
for rule=0:1:21
    [P_W] = getTrianglePoints([P1;P2;P3],rule);
    [N,~]=size(P_W);
    res=0;
    for i=1:1:N
        res=res+func(P_W(i,1),P_W(i,2),P_W(i,3))*P_W(i,4);
    end
    resA(count,1)=res;
    resA(count,2)=rule;
    resA(count,3)=N;
    count=count+1;
end

%% matlab 自帶積分函式
pfun = @(x,y) (x.^6+y.^3+1);
xmin = 0;
xmax = 2;
ymin = 0;
ymax = @(x) -3/2*x+3;
r = integral2(pfun,xmin,xmax,ymin,ymax);

%% plot
figure(22)
plot(resA(:,2),resA(:,1),'r-o');hold on;
plot(resA(:,2),r*ones(22,1),'b-');grid on;
xticks([0:2:22]);
xlim([0,22]);
legend("TRIANGLE LYNESS RULE 積分","Matlab integral2積分");
text(10,14,"積分函式：(x^6+y^3+1)")
text(10,12,"積分割槽域：(0,0,0),(2,0,0),(0,3,0)");
xlabel("Lyness Rule");
ylabel("積分數值");

積分結果對比

image

程式碼

getTrianglePoints.m

function [P_W] = getTrianglePoints(Triangle,Rule)
% getTrianglePoints 三角形面元積分
% https://people.sc.fsu.edu/~jburkardt/cpp_src/triangle_lyness_rule/triangle_lyness_rule.html
%   輸入：
%   Triangle(3,3)：三角形面元三個點
%   Rule：triangle_lyness_rule
%   輸出：
%   P_W(:,4)：P_W(:,1:3)積分點、P_W(:,4)權重係數

%% 任意空間三角形 =》平面直角三角形 座標轉換
% 形函式
N1=@(s,t) -s-t+1;
N2=@(s,t) s;
N3=@(s,t) t;
N1_s=@(s,t) -1;
N2_s=@(s,t) 1;
N3_s=@(s,t) 0;
N1_t=@(s,t) -1;
N2_t=@(s,t) 0;
N3_t=@(s,t) 1;

P1=Triangle(1,:);
P2=Triangle(2,:);
P3=Triangle(3,:);

global TriCoeff;
data=TriCoeff{Rule+1,1};
[order,~]=size(data);
P_W=zeros(order,4);

for i=1:1:order
    P_W(i,1:3)=Loc2Glo(data(i,1:2));
    P_W(i,4)=data(i,3)*Jacobi(data(i,1:2));
end
    function Pglobal=Loc2Glo(loc)
        %loc(1,2)
        Pglobal=N1(loc(1),loc(2))*P1+...
            N2(loc(1),loc(2))*P2+...
            N3(loc(1),loc(2))*P3;
    end

    function J=Jacobi(Loc)
        s=N1_s(Loc(1),Loc(2))*P1+...
            N2_s(Loc(1),Loc(2))*P2+...
            N3_s(Loc(1),Loc(2))*P3;
        t=N1_t(Loc(1),Loc(2))*P1+...
            N2_t(Loc(1),Loc(2))*P2+...
            N3_t(Loc(1),Loc(2))*P3;
        %三角形，這裡多除了一個2
        J=norm(cross(s,t))/2;
    end
end

getTrianglePointsSimplified.m

function [P_W] = getTrianglePointsSimplified(Triangle,Rule)
    global TriCoeff;
    points = TriCoeff{Rule+1};
    weights = points(:,3)';
    points(:,3) = 1-points(:,1)-points(:,2);
    P_W = zeros(size(points,1),4);
    P_W(:,1:3)=points*Triangle;
    area = 0.5*norm(cross(Triangle(1,:)-Triangle(2,:),Triangle(1,:)-Triangle(3,:)));
    P_W(:,4)=weights*area;
end

loadLynessFromTxT.m

function TriCoeff = loadLynessFromTxT(filename)
%LOADLYNESSFROMTXT 載入係數
TriCoeff=cell(22,1);
fp=fopen(filename,'r');
data=textscan(fp,"%f,%f,%f");
fclose(fp);
ALL=[data{1,1},data{1,2},data{1,3}];
[N,~]=size(ALL);
i=1;
while i<N
   rule=ALL(i,1);
   order=ALL(i,2);
   TriCoeff{rule+1,1}=ALL(i+1:i+order,:);
   i=i+order+1;
end
end

相關文章