变步长梯形求积公式C语言_复合梯形求积公式：算例

本文展示复合梯形求积公式的数值积分效果。 对于闭区间上的一般函数，利用复合梯形求积公式具有二阶精度。 但是对于周期函数，特别是解析函数，以及

上迅速衰减的函数，复合梯形公式具有几何收敛阶。

实轴上的积分

我们首先考虑定义在实轴

上的函数

的积分

如果

是光滑，并且衰减的，那么使用复合梯形公式积分将有几何阶的收敛性。

具体来说，数值积分公式为

引理： Suppose

is analytic in the strip

for some

. Suppose further that

uniformly as

in the strip, and for some

, it satisfies

for all

. Then, for any

,

exists and satisfies

见《The Exponentially Convergent Trapezoidal Rule》。

我们首先来测试一下求积公式对“不那么光滑”函数的积分效果

例子一

考虑有限区间上的函数

积分的精确值为

积分节点为等距节点

步长为

。 复合梯形积分公式为（因为函数在端点处值相同（周期））

Matlab程序

%% 复合梯形公式求积函数
%  4-x^2, x in [-2,2]

clear all
close all
p_set = 2.^(4:10);
w = @(z) 4-z.^2;
ErrInt = zeros(size(p_set));
% 精确值
ExactInt = integral(w,-2,2);

figure; set(gcf,'unit','centimeters','position',[10,10,24,12]);
for ii = 1:length(p_set)
    p = p_set(ii);
    xk = linspace(-2,2,p+1);
    xk = xk(2:end);
    ErrInt(ii) = abs(4*mean(w(xk))-ExactInt);
end
loglog(p_set,ErrInt,'-.o','LineWidth',1.5);
hold on;
legend('2阶')

横坐标是节点个数

，在loglog图像下，数值积分误差随着节点增加而二阶收敛。

例子二

我们选择

复合梯形公式的收敛为

• 左侧的图中，数值积分误差随着节点增加几何收敛。此时节点数目较少。
• 右侧图中，数值积分误差随着节点增加
  阶收敛。此时节点数目较多。

虽然上面例子中的函数都是给定区间上的周期函数，但它们周期化之后并不是解析函数。另一个角度则是，积分收敛或与它们的傅里叶系数的衰减有关。

例子三

被积函数选择为

其中

和

是参数。精确的积分值为

考虑

，

的参数情形。

之前对于这类无穷区间上的积分，我们采用的是先截断到固定区间上的做法。 此处，我们将数值积分的参数设置为步长

以及节点数

。 那么之前固定区间就等价于

。

下面是复合梯形公式的误差结果关于

的收敛，左右分别是不同的截断

的选取方法。

有以下两点观察：

• 误差在
  达到
  左右剧烈减小
• 和
  没有明显区别

第一点可以通过分析

的Fourier变换的带宽来解释。

是个Gauss函数，带宽为

Fourier变换

依旧是Gauss函数，带宽为

如果复合梯形公式的步长

满足

， 即

， 则无穷项的复合梯形公式几乎是精确的。

因此，在这种情况下，有限截断的复合梯形公式的误差主要来自于截断项。 但是被积函数

的带宽

非常小，因此截断的误差几乎没有影响。

最后，在右侧的图像中，从右往左看，倒数第二个点没有画出来，这是因为在Matlab中计算得到的误差为

的原因。

围道积分的应用

给定函数

延拓到复平面，它有两个单极点

和

。

我们希望已知极点，求出它们对应的系数（留数）。

在Matlab中，我们可以通过符号计算函数 limit 来执行这个操作。通过帮助文件可得

limit(F,x,a) takes the limit of the symbolic expression F as x -> a.

具体举例如下：

syms x
w = 1/(x^2-x-2);
coef1 = limit(w*(x-2),x,2);
disp(coef1)

借助于复合梯形公式，我们也可以利用围道积分来计算

复合梯形公式得到的近似为

Matlab测试程序

p_set = 4:10;
w = @(z) 1./((z-2).*(z+1));
ErrInt = zeros(size(p_set));
r_set = [1,0.5,0.25];

figure; set(gcf,'unit','centimeters','position',[10,10,24,12]);
for jj = 1:length(r_set)
    r = r_set(jj);
    for ii = 1:length(p_set)
        p = p_set(ii);
        k = 1:p;
        ErrInt(ii) = abs(real(mean(r*w(2+r*exp(2*pi*1i*k/p)).*exp(2*pi*1i*k/p)))-1/3);
    end
    semilogy(p_set,ErrInt,'LineWidth',1.5);
    hold on;
end
legend('r=1','r=0.5','r=0.25')

分别测试了三种半径对应的围道，其中

情况下，

个采样点（意味着数值积分中是10次乘法和加法）的精度达到了

。

相比来说，符号运算 limit 函数不需要我们费心。但是如果计算资源受限，围道积分是一种选择。