目录
标准区间
一般区间
数值实验
实验一
实验二
总结
下节预告
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 表.
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}\)
| \(A_1=5/9\)
| 5 |
4 | \(x_{1}=-\sqrt{\dfrac{15+2\sqrt{30}}{35}}\)
| \(A_1=\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 点, 代入公式即可计算积分值.
数值实验
所有实验在 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 |
总结
随着 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
来源: https://www.cnblogs.com/sunzhenwei/p/10847059.html