1. 程式人生 > >玩轉 matlab 之二維 gauss 數值積分公式使用及 matlab 原始碼(1)-常量區間

玩轉 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

四. 總結和下節預告

  1. 從實驗資料可以發現,二重gauss數值積分使用\(2\times 2\) 4個點的情況下,結果已經很準確了,達到了1e-6的誤差,所以在實際計算中,一般使用4點或者9點的gauss公式就可以滿足要求了。
  2. 下一節介紹當積分割槽間在變係數下的二重gauss公式的計算方法
  3. 歡迎與我交流,數值分析,矩陣計算,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