玩转matlab之一维 gauss 数值积分公式及matlab源代码

在数值分析中，尤其是有限元刚度矩阵、质量矩阵等的计算中，必然要求如下定积分:

\[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 ```\]