在数值分析中，尤其是有限元刚度矩阵、质量矩阵等的计算中，必然要求如下定积分:
\[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表.
一般区间
\[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点，代入公式即可计算积分值。
数值实验
所有实验在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计算结果为：
实验二
取函数
\[f(x)=\dfrac{x^2+2x+1}{1+(1+x)^4},\]
使用matlab函数 integral 计算得到: \(I= 0.161442165779443\)。
使用gauss积分的matlab计算结果为：
总结
- 随着gauss点m的个数增多，精度在逐渐提高，但是要注意的是，gauss点取得多的话，计算量也会增大，只是因为我们计算的问题规模比较小，所以感觉不到而已。
- 另外可以看到2点3点的gauss公式的精度已经很高了，说明并不需要取太多的点，而在实际计算中，比如有限元的计算中，也仅仅取2点或者3点gauss积分就完全足够。
下节预告
下次介绍gauss积分的二维公式使用以及matlab数值实验，欢迎有问题与我交流，偏微分方程，矩阵计算，数值分析等问题，我的qq 群 315241287
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