玩轉 matlab 之二維 gauss 數值積分公式使用及 matlab 原始碼(1)-常量區間
目錄
- 一. 標準區間
- 二. 一般區間
- 三. 數值實驗(沒有程式設計的計算是不完整的)
- 四. 總結和下節預告
- 五. matlab原始碼
繼續上一篇一維gauss積分的討論，本文討論二維gauss數值積分公式的使用，並給出數值實驗。
一. 標準區間
按照一維情況類似方式，首先給出如下二重積分公式：
\[
I = \int_{-1}^1dx\int_{-1}^1 f(x,y)dy=\Sigma_{k=1}^m\Sigma_{l=1}^mA_{kl}f(x_k,y_l).
\]
這是標準的二維gauss積分 \(m\times m\) (即 x 軸和 y 軸分別取 m 個點對應\(m^2\)個點)公式的表示式，其中，\(x_{kl}\) 稱作gauss積分點，\(A_{kl}\)稱作\(x_{kl}\) 點處的權重。\(m=1,2,3\) 時候草圖如下所示：
下面給出部分gauss積分點和權重對應列表
| gauss點個數 \(m^2(m\times m)\) | gauss 點 \((x_i,y_i)\) | 權重 \(A_i\) | 精度(2m-1) |
|---|---|---|---|
| 1 (\(1\times 1\)) | \(x_1=y_1=0\) | \(A_1\)=4 | 1 |
| 4 (\(2 \times 2\)) | \(x_{1}=-1/\sqrt{3},y_{1}=-1/\sqrt{3}\) \(x_{2}=1/\sqrt{3},y_{2}=-1/\sqrt{3}\) \(x_{3}=1/\sqrt{3},y_{3}=1/\sqrt{3}\) \(x_{4}=-1/\sqrt{3},y_{4}=1/\sqrt{3}\) | \(A_1=A_2=A_3=A_4=1\) | 3 |
| 9 (\(3 \times 3\)) 令\(gpt=\sqrt{3/5}\) |
\(x_1=-gpt,y_1=-gpt\) \(x_2=gpt,y_2=-gpt\) \(x_3=gpt,y_3=gpt\) \(x_4=-gpt,y_4=gpt\) \(x_5=0,y_5=-gpt\) \(x_6=gpt,y_6=0\) \(x_7=0,y_7=gpt\) \(x_8=-gpt,y_8=0\) \(x_9=0,y_9=0\) |
\(A_1=A_2=A_3=A_4=25/81\) \(A_5=A_6=A_7=A_8=40/81\) \(A_9=64/81\) |
5 |
一般二重積分近似值也就是使用 \(2\times 2,\quad 3\times 3\) 公式就完全足夠了。所以不需要太多的點，徒增計算量。因此就可以得到gauss積分的座標表示式為:
\[
I = \int_{-1}^1dx\int_{-1}^1 f(x,y)dy=\Sigma_{i=1}^{m^2}A_{i}f(x_i,y_i).
\]
二. 一般區間
對於一般區間，先考慮區間端點為常量情況（下一節介紹區間為變數情況），即
\[
I = \int_{a}^bdx\int_{c}^d f(x,y)dy.
\]
其中 \(a,b,c,d\) 都是已知常數。與一維情況類似，只需要做變數變換，於是\([s,t]\in[-1,1]\times [-1,1]\)（通俗講就是換元法）
\[
令 \quad x =\frac{b+a+(b-a)s}{2},\qquad 則\quad dx=\frac{b-a}{2}ds,\\
\quad y =\frac{d+c+(d-c)t}{2},\qquad 則\quad dy=\frac{d-c}{2}dt.\\
記\qquad jac = \frac{(b-a)(d-c)}{4},\qquad 則 \quad dxdy=jac\times dsdt
\]
於是二重積分就變成了
\[
I = \int_{a}^bdx\int_{c}^d f(x,y)dy
=jac\times\int_{-1}^1ds\int_{-1}^1f(\frac{b+a+(b-a)s}{2},\frac{d+c+(d-c)t}{2})dt.\\=jac\times\Sigma_{i=1}^{m^2}A_{i}f(\frac{b+a+(b-a)s_i}{2},\frac{d+c+(d-c)t_i}{2}).
\]
其中 \((s_i,t_i)\) 即為上表中的 gauss 節點，對應的權重因子為 \(A_i\)。
三. 數值實驗(沒有程式設計的計算是不完整的)
使用matlab2018a 計算結果，並且與matlab自帶函式 integral2 計算的結果進行比較給出誤差。
算例如下:
\[
計算定積分 I = \int_{a}^bdx\int_{c}^d f(x,y)dy
\]
其中，\(a=1.4,b=2,c=1,d=1.5,f(x,y)=ln(x+2*y), ln\)是以e為底對數函式。使用matlab的integral2 函式計算結果為\(I =0.429554527548275\).
自己程式設計計算結果如下：
| 高斯點數 \(m^2(m\times m)\) | 積分值 \(I_m\) | 誤差norm(\(I_m-I\)) |
|---|---|---|
| 4(\(2\times 2\)) | 0.429556088022242 | 1.56E-06 |
| 9(\(3\times 3\)) | 0.429554531152490 | 3.60E-09 |
四. 總結和下節預告
- 從實驗資料可以發現，二重gauss數值積分使用\(2\times 2\) 4個點的情況下，結果已經很準確了，達到了1e-6的誤差，所以在實際計算中，一般使用4點或者9點的gauss公式就可以滿足要求了。
- 下一節介紹當積分割槽間在變係數下的二重gauss公式的計算方法
- 歡迎與我交流，數值分析，矩陣計算，PDE數值解等，QQ群 315241287
五. matlab原始碼
clc;clear;
% compute int_a^b [int_c)^d f(x,y)]
% (x,y) \in [a,b] X [c,d]
%% setup the integral interval and gauss point and weight
a = 1.4; b = 2;
c = 1; d =1.5;
fun=@(x,y) log(x+2*y);
fprintf('***********************************************\n')
for gauss = 2:3 % m points rule in 2 dimensional case
if gauss == 2
fprintf('******* 2X2 points gauss rule result *******')
gpt=1/sqrt(3);
s(1) = -gpt; t(1) = -gpt;
s(2) = gpt; t(2) = -gpt;
s(3) = gpt; t(3) = gpt;
s(4) = -gpt; t(4) = gpt;
wt = [1 1 1 1];
elseif gauss == 3
gpt=sqrt(0.6);
fprintf('******* 3X3 points gauss rule *******')
s(1) = -gpt; t(1) = -gpt; wt(1)=25/81;
s(2) = gpt; t(2) = -gpt; wt(2)=25/81;
s(3) = gpt; t(3) = gpt; wt(3)=25/81;
s(4) = -gpt; t(4) = gpt; wt(4)=25/81;
s(5) = 0.0; t(5) = -gpt; wt(5)=40/81;
s(6) = gpt; t(6) = 0.0; wt(6)=40/81;
s(7) = 0.0; t(7) = gpt; wt(7)=40/81;
s(8) = -gpt; t(8) = 0.0; wt(8)=40/81;
s(9) = 0.0; t(9) = 0.0; wt(9)=64/81;
end
%% 區間變換到 [-1,1] X [-1,1]
jac = (b-a)*(d-c)/4;
x = (b+a+(b-a)*s)/2;
y = (d+c+(d-c)*t)/2;
f = fun(x,y);
comp = wt(:) .* f(:) .* jac;%無論一個向量是行還是列，寫成x(:)都會變成列向量
format long
comp = sum(comp)
exact = integral2(fun,a,b,c,d);
fprintf('the error is norm(comp-exact)=%10.6e\n\n',norm(comp-exact))
end
fprintf('******************************************\n')
fprintf('matlab built-in function ''integral2''\n')
exact
format short