标签为'数值积分'的全部存档

复化Gauss-Legendre求积

发表于： 2010年05月28日 分类：数字之美 0条评论 标签: , , , .

高斯型求积公式是数值积分中一个比较成熟的想法. 速度快, 精度高, 公式优雅.

一般Gauss-Legendre求积公式

对于一般的Gauss型求积公式

\int_a^b f(x) \rho (x)\,dx \approx \sum_{k=0}^{n} A_k f(x_k)

取权函数\rho (x)=1, 积分区间设定为\left[-1, 1\right], 则得到一般的Gauss-Legendre求积公式
\int_{-1}^{1} f(x) dx \approx \sum_{k=0}^{n} A_k f(x_k).

取Legendre多项式的零点作为Gauss点, 通过计算得到Gauss点个数n=2和n=3时的求积公式
\int_{-1}^{1} f(x)\,dx \approx f(-\frac{1}{\sqrt{3}})+f(\frac{1}{\sqrt{3}})

以及
\int_{-1}^{1} f(x)\,dx \approx \frac{5}{9} f(-\frac{\sqrt{15}}{5})+\frac{8}{9}f(0)+\frac{5}{9}f(\frac{\sqrt{15}}{5}).

变量替换方法

上述分析中, 积分区间固定为\left[-1, 1\right], 实际应用时做变量替换

x=\frac{a+b}{2}+\frac{b-a}{2}t.

将被积区间\left[a, b\right]化为\left[-1, 1\right].

复化Gauss-Legendre求积公式

将被积区间m等分, 记h=\frac{b-a}{m}, x_k=a+kh, k=0, 1, 2 \ldots, m. 作变换

x=\frac{x_k+x_{k-1}}{2}+\frac{h}{2}\,t,

在每个小区间上应用Gauss-Legendre公式, 累加即得复化Gauss-Legendre求积公式

\int_a^b f(x)\,dx=\frac{h}{2}\sum\sum A_k f(x_k).

不妨设

G(t)=f(\frac{x_k+x_{k+1}}{2}+\frac{h}{2}\,t),

则有:

Gauss点个数n=2时,

\int_a^b f(x)\,dx\approx\frac{h}{2}\sum G(-\frac{1}{\sqrt{3}})+G(\frac{1}{\sqrt{3}}),

Gauss点个数n=3时,

\int_a^b f(x)\,dx\approx\frac{h}{2}\sum \frac{5}{9}G(-\frac{\sqrt{15}}{5})+\frac{8}{9}G(0)+\frac{5}{9}G(\frac{\sqrt{15}}{5}).

总结复化Gauss-Legendre求积过程如下:

1. 分割区间, 记录区间端点值；
2. 通过查表或求解非线性方程组, 在所有小区间上, 将Gauss系数和Gauss点的值代入变量替换后的公式；
3. 将所有区间的结果累加, 即得到整个区间上的积分近似值.

针对Gauss点个数n=2和n=3的复化Gauss-Legendre求积公式编写的一个简单的MATLAB函数 compgauss() 如下:

?Download compgauss.m
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34

function [ ] = compgauss(a, b, n)
% Composite Gauss Integration
% Equation Type: n=2, n=3
% Coded by Nan.Xiao 2010-05-25
% Step.1 Divide Interval
% Step.2 Calculate
% Step.3 Sum Results
format long
f = @(x) exp(x).*sin(x);
h=(b-a)/n;
xk=zeros(n+1,1);
xk(1,1)=a;
xk(n+1,1)=b;
fk1=zeros(n,1);
fk2=zeros(n,1);
for i=1:n-1
    xk(i+1,1)=a+h*i;
end
for j=1:n
    fk1(j)=f((xk(j)+xk(j+1))/2+(h/2)*(-1/sqrt(3)))+...
        f((xk(j)+xk(j+1))/2+(h/2)*(1/sqrt(3)));
end
for r=1:n
    fk2(r)=(5/9)*f((xk(r)+xk(r+1))/2+(h/2)*(-sqrt(15)/5))+...
        (8/9)*f((xk(r)+xk(r+1))/2+(h/2)*(0))+...
        (5/9)*f((xk(r)+xk(r+1))/2+(h/2)*(sqrt(15)/5));
end
mysum1=h*sum(fk1)/2;
mysum2=h*sum(fk2)/2;
disp('Result of 2 Nodes:')
disp(mysum1);
disp('Result of 3 Nodes:')
disp(mysum2);
end

总结

1. Gauss求积公式较一般的机械求积公式的进步之处, 在于其针对插值型求积公式进行了改进. Gauss公式将插值节点设为未知, 成功地将代数精度由n次提高到2n+1次. 计算过程简单, 速度快, 达到要求精度所需步骤较少.

2. (复化) Gauss 求积的一个问题在于, 增加Gauss点个数, 继续求解Gauss系数和Gauss点值时, 需要解一系列非线性方程组. 其解析解较难求得. 而取数值解作为公式中的参数, 不如解析解理想.