matlab一维数据扩充,玩转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表.

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}\\)

\\(x_2=0\\)

\\(x_3=\\sqrt{3/5}\\)

\\(A_1=5/9\\)

\\(A_2=8/9\\)

\\(A_3=5/9\\)

5

4

\\(x_{1}=-\\sqrt{\\dfrac{15+2\\sqrt{30}}{35}}\\)

\\(x_{2}=-\\sqrt{\\dfrac{15-2\\sqrt{30}}{35}}\\)

\\(x_{3}=\\sqrt{\\dfrac{15-2\\sqrt{30}}{35}}\\)

\\(x_{4}=\\sqrt{\\dfrac{15+2\\sqrt{30}}{35}}\\)

\\(A_1=\\frac{90-5\\sqrt{30}}{180}\\)

\\(A_2=\\frac{90+5\\sqrt{30}}{180}\\)

\\(A_3=\\frac{90+5\\sqrt{30}}{180}\\)

\\(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点，代入公式即可计算积分值。

数值实验

所有实验在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