高斯-勒朗德積分公式
高斯-勒朗德積分原理參考《數值分析》第五版P188
需求：給定空間平面\(S\)四個點的座標\(Q_1(x,y,z),Q_2(x,y,z),Q_3(x,y,z),Q_4(x,y,z)\)，已知函式\(f(x,y,z)\)，求利用數值方法求解積分：\(\iint_Sf(x,y,z)\text dS\)。
解決方法：採用高斯-勒朗德積分方法進行求解。
計算步驟：
- 轉換座標至引數座標系，計算高斯點
- 計算積分
座標變換
\((x,y,z)\)與\((s,t,0)\)分別表示總體t座標系與轉換之後的引數座標系，通過在單元內引入合適的形函式（又稱插值函式）\(N_k(s, t)(k=1, 2,...,NE)\)並進行如下座標變換：
\[\begin{align}(x,y,z)=\sum_{k=1}^{NE}N_k(s,t)(x_k,y_k,z_k) \quad s,t\in[-1,1]\end{align}
\]
可以將物理空間\((x,y,z)\)的空間四邊形單元轉化為引數空間\((s,t,0)\)中的正方形單元，\(NE\)與\((x_k,y_k,z_k)\)分別代表單元節點數與第\(k\)個節點的座標。在單元內，形函式\(N_k(s,t)\)應該滿足以下兩個條件：
\[\begin{align}
N_k(s,t)&=\begin{cases}
1&在節點k上\\
0&在其他節點上
\end{cases}\\
\sum_{k=1}^{NE}N_k(s,t)&=1
\end{align}
\]
待定係數法，求解得到四節點形函式為
\[\begin{cases}
N_1(s,t)=0.25(1+s)(1+t)\\
N_2(s,t)=0.25(1-s)(1+t)\\
N_3(s,t)=0.25(1-s)(1-t)\\
N_4(s,t)=0.25(1+s)(1-t)
\end{cases}\\
\]
四節點等參單元的形函式是關於 $s, t \(的線性函式，而九節點等參單元的形函式是關於\) s, t $的二次函式。
積分座標系轉換
\[\begin{equation} \begin{aligned}
\iint_{S_{}}f(x,y,z)\text dS(x,y,z)&=\iint_{S_{0}}f(s,t,0)|J(s,t)|\text dS(s,t,0)\\
&=\sum_{i=0}^{n}\sum_{j=0}^{n}A_iA_jf(s_i,t_j,0)|J(s_i,t_j)|\\
\end{aligned} \end{equation}
\]
其中\(|J(s,t)|\)為\(Jacobi\)行列式
\[|J(s,t)|=\left| \frac{\boldsymbol r}{\partial s}\cross\frac{\boldsymbol r}{\partial t}\ \right |
\]
式中\(\boldsymbol r=(x,y,z)\)與\(s,t\)的關係可以寫作
\[\boldsymbol r=(x,y,z)=\sum_{k=1}^{NE}N_k(s,t)(x_k,y_k,z_k)
\]
至此已經完成高斯勒朗德積分法的求解主要步驟。
另外，可將係數與函式值分離（避免多次計算高斯點），預先求解每個空間四邊形面的高斯點位置與係數（係數與高斯點對應係數與座標變換Jacobi係數矩陣相關）
測試例項
例項1
某二重積分
\[\int_{1.4}^{2}\int_{1}^{1.5} \text{ln}(x+2y)dydx\approx 0.42955453
\]
變換為區域\(R=\{(s,t)|-1\le s,t\le 1\}\)
\[s=\frac{1}{0.6}(2x-3.4),t=\frac{1}{0.5}(2y-2.5)\\
x=0.3s+1.7,y=0.25t+1.25
\]
\(n=2\)時的高斯積分節點與係數
| \(i\) | 0 | 1 | 2 |
|---|---|---|---|
| \(s_i,t_i\) | \(-0.774596662\) | \(0\) | \(0.774596662\) |
| \(A_i\) | \(\frac{5}{9}\) | \(\frac{8}{9}\) | \(\frac{5}{9}\) |
計算結果
%demo
clc;clear
Q=[1.4 1 0;
2 1 0;
2 1.5 0;
1.4 1.5 0];
n=2;
[Gauss_P_A] = GetGaussPoint(Q,n);
P0=[0,0,3];
Func=@(Q) log(Q(1)+2*Q(2));
res=0;
for i=1:1:n*n
res=res+Func(Gauss_P_A(i,1:3))*Gauss_P_A(i,4)*Gauss_P_A(i,5);
end
res
例項2
計算如下積分式：
clc;clear
Q=[0 0 0;
1 0 0;
1 1 0;
0 1 0];
for n=1:1:4
% n=2;
[Gauss_P_A] = GetGaussPoint(Q,n);
P0=[0,0,1];
Func=@(Q) norm(Q-P0);
res=0;
for i=1:1:n*n
res=res+Func(Gauss_P_A(i,1:3))*Gauss_P_A(i,4)*Gauss_P_A(i,5);
end
re(n)=res;
end
plot(1:1:n,re,'o');
| n | 積分數值 |
|---|---|
| 0 | 1.224744871391589 |
| 1 | 1.280924113061923 |
| 2 | 1.280797365089580 |
| 3 | 1.280788916595401 |
例項3
\[\frac{1}{r(P,Q)}=\frac{1}{\sqrt{(x_P-x_Q)^2+(y_P-y_Q)^2+(z_P-z_Q)^2}}
\]
對\(P\)的全導數
\[\begin{align}
\nabla{G(P,Q)}&=\nabla{\frac{1}{r(P,Q)}}\\&=-\frac{1}{r^3}(x_P-x_Q,y_P-y_Q,z_P-z_Q)\\&=-\frac{1}{r^3}(P-Q)
\end{align}
\]
GetGaussPoint.m
function [Gauss_P_A] = GetGaussPoint(Q,n)
%GETGAUSSPOINT 計算高斯點
%輸入:Q面元點四個點(4*3),高斯點階數n
%輸出:Gauss_P_A:[高斯點位置,高斯點係數];
Nk1=@(s,t)0.25*(1+s)*(1+t);
Nk2=@(s,t)0.25*(1-s)*(1+t);
Nk3=@(s,t)0.25*(1-s)*(1-t);
Nk4=@(s,t)0.25*(1+s)*(1-t);
Nk1_s=@(s,t)0.25*( 1)*(1+t);
Nk2_s=@(s,t)0.25*(-1)*(1+t);
Nk3_s=@(s,t)0.25*(-1)*(1-t);
Nk4_s=@(s,t)0.25*( 1)*(1-t);
Nk1_t=@(s,t)0.25*(1+s)*( 1);
Nk2_t=@(s,t)0.25*(1-s)*( 1);
Nk3_t=@(s,t)0.25*(1-s)*(-1);
Nk4_t=@(s,t)0.25*(1+s)*(-1);
%高斯點對應係數
Gauss_Loc{1,1}=[+0,2];
Gauss_Loc{2,1}=[+0.577350300000000,1;
-0.577350300000000,1];
Gauss_Loc{3,1}=[+0.774596700000000,0.555555600000000;
-0, 0.888888900000000;
-0.774596700000000,0.555555600000000];
Gauss_Loc{4,1}=[+0.861136300000000,0.347854800000000;
+0.339981000000000,0.652145200000000;
-0.339981000000000,0.652145200000000;
-0.861136300000000,0.347854800000000];
The_Gauss=Gauss_Loc{n,1};
[Num,~]=size(The_Gauss);
Gauss_P_A=zeros(Num*Num,5);
count=1;
for i=1:1:Num
for j=1:1:Num
P_Loc=[The_Gauss(i,1),The_Gauss(j,1)];
Gauss_P_A(count,1:3)=L2G(P_Loc);
J=Jacobi(P_Loc);
Gauss_P_A(count,4)=The_Gauss(i,2)*The_Gauss(j,2);
%Jacobi係數
Gauss_P_A(count,5)=J;
count=count+1;
end
end
function Glo=L2G(Loc)
Glo=Nk1(Loc(1),Loc(2))*Q(1,:)+...
Nk2(Loc(1),Loc(2))*Q(2,:)+...
Nk3(Loc(1),Loc(2))*Q(3,:)+...
Nk4(Loc(1),Loc(2))*Q(4,:);
end
function J=Jacobi(Loc)
s=Nk1_s(Loc(1),Loc(2))*Q(1,:)+...
Nk2_s(Loc(1),Loc(2))*Q(2,:)+...
Nk3_s(Loc(1),Loc(2))*Q(3,:)+...
Nk4_s(Loc(1),Loc(2))*Q(4,:);
t=Nk1_t(Loc(1),Loc(2))*Q(1,:)+...
Nk2_t(Loc(1),Loc(2))*Q(2,:)+...
Nk3_t(Loc(1),Loc(2))*Q(3,:)+...
Nk4_t(Loc(1),Loc(2))*Q(4,:);
J=norm(cross(s,t));
end
end