在數值分析中，尤其是有限元剛度矩陣、質量矩陣等的計算中，必然要求如下定積分:
\[I=\int_a^b f(x)dx $$學好**gauss**積分也是學好**有限元**的重要基礎，學過高等數學的都知道，手動積分能把人搞死（微笑臉），而且有些函數還不存在原函數，使用原始的手動算出原函數幾乎是不現實的。因此非常有必要學習數值積分，簡單講就是近似計算，只要這個近似值精確度高和穩定性好就行。Gauss積分公式就是這么一個非常好用的工具。本文介紹高斯積分公式的使用以及簡單的數值算例。 # 標准區間 先考慮特殊情況，對於一般區間呢？待會會處理這個問題。 \]
I=\int_{-1}^1 f(x)dx
\[不加證明的直接給出gauss公式如下：詳情參閱任何一本數值分析書都有詳細的證明過程： \]
I=\int_{-1}^1 f(x)dx=\Sigma_{i=1}^n A_if(x_i)
\[其中$A_i$稱作**權**，$x_i$稱作 **gauss 點**。 下面的問題就是如何選擇$n,A_i,x_i$。 理論表明**n**個點的Gauss公式代數精度為$2n-1$，其選擇如下表，（這里僅僅舉1-4個點情況，實際使用的時候一般2點或者3點的精度已經完全夠了）更多積分點可參考 [gauss表](http://blog.sina.com.cn/s/blog_67c5d9870102xafw.html). | gauss點個數 $n$ | gauss 點 $x_i$ | 權重 $A_i$ | 精度 | | ----------------- | ------------------------------------------------------------ | ------------------------------------------------------------ | ---- | | 1 | $x_1$=0 | $A_1$=2 | 1 | | 2 | $x_{1,2}=\pm1/\sqrt{3}$ | $A_1=A_2=1$ | 3 | | 3 | $x_1=-\sqrt{3/5}$<br>$x_2=0$<br>$x_3=\sqrt{3/5}$ | $A_1=5/9$<br>$A_2=8/9$<br>$A_3=5/9$ | 5 | | 4 | $x_{1}=-\sqrt{\dfrac{15+2\sqrt{30}}{35}}$<br>$x_{2}=-\sqrt{\dfrac{15-2\sqrt{30}}{35}}$<br>$x_{3}=\sqrt{\dfrac{15-2\sqrt{30}}{35}}$<br>$x_{4}=\sqrt{\dfrac{15+2\sqrt{30}}{35}}$ | $A_1=\frac{90-5\sqrt{30}}{180}$<br>$A_2=\frac{90+5\sqrt{30}}{180}$<br>$A_3=\frac{90+5\sqrt{30}}{180}$<br>$A_4=\frac{90-5\sqrt{30}}{180}$ | 7 | # 一般區間 \]
I=\int_a^b f(x)dx
\[ 根據上面的討論情況，可知只要做變換(相當於換元積分一樣) \]
令\quad x=\frac{b+a+(b-a)s}{2},\
則\quad dx = \frac{b-a}{2}ds.
\[那么有$s\in[-1,1]$，於是即可使用標准區間公式如下： \]
I = \int_a^bf(x)dx=\int_{-1}^1f(\frac{b+a+(b-a)s}{2})\times\frac{b-a}{2}ds\
=\frac{b-a}{2}\Sigma_{i=1}^mA_if(\frac{b+a+(b-a)s_i}{2})
\[上述公式中的$A_i$即為表格中的權重，$s_i$即為上表中對應的gauss點，代入公式即可計算積分值。 # 數值實驗 [toc] 所有實驗在MATLAB2018a版本下完成。（建議安裝新版本，因為很多函數在新版中已經優化了或者改名字了，比如老版本積分函數**quad** 新版已經改為**integral**，只不過目前quad函數還是可以使用的，將來會被刪除）。 我們取2個函數做實驗，分別計算出其gauss積分值再與matlab自帶的函數 **integral** 計算結果作比較，實驗模型是： \]
計算 \quad I= \int_1^2 f(x)dx
\[ ## 實驗一 取函數 \]
f(x)=lnx, \quad 即自然對數函數以e為底.
\[使用matlab函數 integral 計算得到: $I= 0.386294361119891$。 使用gauss積分的matlab計算結果為： | 高斯點數 m | 積分值 $I_m$ | 誤差norm($I_m-I$) | | ---------- | ----------------- | ----------------- | | 2 | 0.386594944116741 | 3.01E-04 | | 3 | 0.386300421584011 | 6.06E-06 | | 4 | 0.386294496938714 | 1.36E-07 | | 5 | 0.386294364348948 | 3.23E-09 | ## 實驗二 取函數 \]
f(x)=\dfrac{x^2+2x+1}{1+(1+x)^4},
\[使用matlab函數 integral 計算得到: $I= 0.161442165779443$。 使用gauss積分的matlab計算結果為： | 高斯點數 m | 積分值 $I_m$ | 誤差norm($I_m-I$) | | ---------- | ----------------- | ----------------- | | 2 | 0.161394581386268 | 4.76E-05 | | 3 | 0.161442818737102 | 6.53E-07 | | 4 | 0.161442196720137 | 3.09E-08 | | 5 | 0.161442166345131 | 5.66E-10 | # 總結 1. 隨着gauss點m的個數增多，精度在逐漸提高，但是要注意的是，gauss點取得多的話，計算量也會增大，只是因為我們計算的問題規模比較小，所以感覺不到而已。 2. 另外可以看到2點3點的gauss公式的精度已經很高了，說明並不需要取太多的點，而在實際計算中，比如有限元的計算中，也僅僅取2點或者3點gauss積分就完全足夠。 # 下節預告 下次介紹gauss積分的二維公式使用以及matlab數值實驗，歡迎有問題與我交流，偏微分方程，矩陣計算，數值分析等問題，我的qq 群 315241287 # matlab代碼 ```matlab clc;clear; % using 2 3 4 5 points compute the integral % x \in [a,b] % %% test a=1; b = 2; fun = @(x) log(x); % fun = @(x) 2*x./(1+x.^4); % fun = @(x) exp(-x.^2/2); % fun = @(x) (x.^2+2*x+1)./(1+(1+x).^4); %% setup the gauss data for gauss = 2:5 if gauss == 2 s=[-1 1]/sqrt(3); wt=[1 1]; fprintf('*************************** 2 points gauss *******') elseif gauss == 3 s = [-sqrt(3/5) 0 sqrt(3/5)]; wt = [5 8 5]/9; fprintf('*************************** 3 points gauss *******') elseif gauss == 4 fprintf('*************************** 4 points gauss *******') s = [ -sqrt((15+2*sqrt(30))/35), -sqrt((15-2*sqrt(30))/35), ... sqrt((15-2*sqrt(30))/35), sqrt((15+2*sqrt(30))/35)]; wt = [ (90-5*sqrt(30))/180, (90+5*sqrt(30))/180,... (90+5*sqrt(30))/180, (90-5*sqrt(30))/180]; elseif gauss == 5 fprintf('*************************** 5 points gauss *******') s(1)=.906179845938664 ; s(2)=.538469310105683; s(3)=.0; s(4)=-s(2) ; s(5)=-s(1); wt(1)=.236926885056189 ; wt(2)=.478628670499366; wt(3)=.568888888888889 ; wt(4)=wt(2) ; wt(5)=wt(1); end %% % 區間變換到 s \in[-1,1] s = (b-a)/2*s+(b+a)/2; jac = (b-a)/2;% dx = jac * ds f = fun(s); f = wt.* f .* jac; format long exact = integral(fun,a,b); comp = sum(f) fprintf('the error is norm(comp-exact)=%10.6e\n\n\n',norm(comp-exact)) end fprintf('\n\n********* matlab built-in function ''integral''*********\n') exact = integral(fun,a,b) format short ```\]