文章目录

标准区间
一般区间
数值实验
- 实验一
- 实验二
总结
下节预告
matlab代码

在数值分析中，尤其是 有限元刚度矩阵、质量矩阵等的计算中，必然要求如下定积分:
I=∫abf(x)dxI=\int_a^b f(x)dx I=∫abf(x)dx学好 gauss积分也是学好 有限元的重要基础，学过高等数学的都知道，手动积分能把人搞死（微笑脸），而且有些函数还不存在原函数，使用原始的手动算出原函数几乎是不现实的。因此非常有必要学习数值积分，简单讲就是近似计算，只要这个近似值精确度高和稳定性好就行。Gauss积分公式就是这么一个非常好用的工具。本文介绍高斯积分公式的使用以及简单的数值算例。

标准区间

先考虑特殊情况，对于一般区间呢？待会会处理这个问题。
I=∫−11f(x)dxI=\int_{-1}^1 f(x)dx I=∫−11f(x)dx
不加证明的直接给出gauss公式如下：详情参阅任何一本数值分析书都有详细的证明过程：
I=∫−11f(x)dx=Σi=1nAif(xi)I=\int_{-1}^1 f(x)dx=\Sigma_{i=1}^n A_if(x_i) I=∫−11f(x)dx=Σi=1nAif(xi)
其中AiA_iAi称作权，xix_ixi称作 gauss 点。
下面的问题就是如何选择n,Ai,xin,A_i,x_in,Ai,xi。

理论表明n个点的Gauss公式代数精度为2n−12n-12n−1，其选择如下表，（这里仅仅举1-4个点情况，实际使用的时候一般2点或者3点的精度已经完全够了）更多积分点可参考 gauss表.

gauss点个数 nnn	gauss 点 xix_ixi	权重 AiA_iAi	精度
1	x1x_1x1=0	A1A_1A1=2	1
2	x1,2=±1/3x_{1,2}=\pm1/\sqrt{3}x1,2=±1/3	A1=A2=1A_1=A_2=1A1=A2=1	3
3	x1=−3/5x_1=-\sqrt{3/5}x1=−3/5 x2=0x_2=0x2=0 x3=3/5x_3=\sqrt{3/5}x3=3/5	A1=5/9A_1=5/9A1=5/9 A2=8/9A_2=8/9A2=8/9 A3=5/9A_3=5/9A3=5/9	5
4	x1=−15+23035x_{1}=-\sqrt{\dfrac{15+2\sqrt{30}}{35}}x1=−3515+230 x2=−15−23035x_{2}=-\sqrt{\dfrac{15-2\sqrt{30}}{35}}x2=−3515−230 x3=15−23035x_{3}=\sqrt{\dfrac{15-2\sqrt{30}}{35}}x3=3515−230 x4=15+23035x_{4}=\sqrt{\dfrac{15+2\sqrt{30}}{35}}x4=3515+230	A1=90−530180A_1=\frac{90-5\sqrt{30}}{180}A1=18090−530 A2=90+530180A_2=\frac{90+5\sqrt{30}}{180}A2=18090+530 A3=90+530180A_3=\frac{90+5\sqrt{30}}{180}A3=18090+530 A4=90−530180A_4=\frac{90-5\sqrt{30}}{180}A4=18090−530	7

一般区间

I=∫abf(x)dxI=\int_a^b f(x)dx I=∫abf(x)dx
根据上面的讨论情况，可知只要做变换(相当于换元积分一样)
令x=b+a+(b−a)s2,则dx=b−a2ds.令\quad x=\frac{b+a+(b-a)s}{2},\\ 则\quad dx = \frac{b-a}{2}ds. 令x=2b+a+(b−a)s,则dx=2b−ads.
那么有s∈[−1,1]s\in[-1,1]s∈[−1,1]，于是即可使用标准区间公式如下：
I=∫abf(x)dx=∫−11f(b+a+(b−a)s2)×b−a2ds=b−a2Σi=1mAif(b+a+(b−a)si2)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}) I=∫abf(x)dx=∫−11f(2b+a+(b−a)s)×2b−ads=2b−aΣi=1mAif(2b+a+(b−a)si)
上述公式中的AiA_iAi即为表格中的权重，sis_isi即为上表中对应的gauss点，代入公式即可计算积分值。

数值实验

所有实验在MATLAB2018a版本下完成。（建议安装新版本，因为很多函数在新版中已经优化了或者改名字了，比如老版本积分函数quad 新版已经改为integral，只不过目前quad函数还是可以使用的，将来会被删除）。

我们取2个函数做实验，分别计算出其gauss积分值再与matlab自带的函数 integral 计算结果作比较，实验模型是：
计算I=∫12f(x)dx计算 \quad I= \int_1^2 f(x)dx 计算I=∫12f(x)dx

实验一

取函数
f(x)=lnx,即自然对数函数以e为底.f(x)=lnx, \quad 即自然对数函数以e为底. f(x)=lnx,即自然对数函数以e为底.
使用matlab函数 integral 计算得到: I=0.386294361119891I= 0.386294361119891I=0.386294361119891。
使用gauss积分的matlab计算结果为：

高斯点数 m	积分值 ImI_mIm	误差norm(Im−II_m-IIm−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)=x2+2x+11+(1+x)4,f(x)=\dfrac{x^2+2x+1}{1+(1+x)^4}, f(x)=1+(1+x)4x2+2x+1,
使用matlab函数 integral 计算得到: I=0.161442165779443I= 0.161442165779443I=0.161442165779443。
使用gauss积分的matlab计算结果为：

高斯点数 m	积分值 ImI_mIm	误差norm(Im−II_m-IIm−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

总结

随着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:5if gauss == 2s=[-1 1]/sqrt(3);wt=[1 1];fprintf('***************************  2     points gauss  *******')elseif gauss == 3s = [-sqrt(3/5) 0 sqrt(3/5)];wt = [5 8 5]/9;fprintf('***************************    3   points gauss  *******')elseif gauss == 4fprintf('***************************    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 == 5fprintf('***************************    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 * dsf = fun(s);f = wt.* f .* jac;format longexact = 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

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

玩转matlab之一维 gauss 数值积分公式及matlab源代码
目录标准区间一般区间数值实验实验一实验二总结下节预告 matlab代码在数值分析中,尤其是有限元刚度矩阵.质量矩阵等的计算中,必然要求如下定积分: \[ I=\int_a^b f(x) ...
matlab一维矩形积分,玩转matlab之一维 gauss 数值积分公式及matlab源代码
在数值分析中,尤其是有限元刚度矩阵.质量矩阵等的计算中,必然要求如下定积分: $$ I=\int_a^b f(x)dx $$学好gauss积分也是学好有限元的重要基础,学过高等数学的都知道,手动积分能 ...
matlab的积分公式,玩转matlab之一维 gauss 数值积分公式及matlab源代码
释放双眼,带上耳机,听听看~! 目录在数值分析中,尤其是有限元刚度矩阵.质量矩阵等的计算中,必然要求如下定积分: \\[ I=\\int_a^b f(x)dx \\]学好gauss积分也是学好有限元 ...
matlab圆柱内导热分离变量法,一维热传导方程数值解法及matlab实现分离变量法和有限差分法...
一维热传导方程数值解法及matlab实现分离变量法和有限差分法一维热传导方程的Matlab解法分离变量法和有限差分法问题描述实验原理分离变量法实验原理有限差分法实验目的利用分离变量法和有限差分法解热 ...
matlab数值分析拟合实例,数值分析函数拟合matlab代码.doc
数值分析函数拟合matlab代码.doc 第一题MATLAB代码用SPLINE作图XI0204060810YI098092081064038X10012Y1NEWTON3XI,YI,X源代码见M文件Y ...
matlab圆柱内导热分离变量法,一维热传导方程数值解法及matlab实现分离变量法和有限差分法.doc...
一维热传导方程的Matlab解法分离变量法和有限差分法问题描述实验原理分离变量法实验原理有限差分法实验目的利用分离变量法和有限差分法解热传导方程问题利用matlab进行建模构建图形研究 ...
分段二次插值的matlab程序,一维优化方法之二次插值法matlab程序
程序为: a=0; b=pi/2; a1=a; a3=b; a2=a; a4=b; while a3-a1>0.01 f1=-atan(2*sin(2*a1)/(5-cos(2*a1))); f ...
抛物线求积公式求积分算法matlab,基于Matlab的数值积分公式问题.doc
基于Matlab的数值积分公式问题数值分析学号: 学生姓名 :教师 : 教师 2数值积分算法介绍............................................. ...
abel数值反演的matlab实现,abel变换数值反演的积分算子方法.pdf
abel变换数值反演的积分算子方法.pdf 还剩 4页未读, 继续阅读下载文档到电脑,马上远离加班熬夜! 亲,喜欢就下载吧,价低环保! 内容要点: V01．29(2009)No．3数学杂志J．of ...

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