本文展示复合梯形求积公式的数值积分效果。 对于闭区间上的一般函数,利用复合梯形求积公式具有二阶精度。 但是对于周期函数,特别是解析函数,以及
上迅速衰减的函数,复合梯形公式具有几何收敛阶。
实轴上的积分
我们首先考虑定义在实轴
上的函数
的积分
如果
是光滑,并且衰减的,那么使用复合梯形公式积分将有几何阶的收敛性。
具体来说,数值积分公式为
引理: 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图像下,数值积分误差随着节点增加而二阶收敛。
例子二
我们选择
复合梯形公式的收敛为
- 左侧的图中,数值积分误差随着节点增加几何收敛。此时节点数目较少。
- 右侧图中,数值积分误差随着节点增加 阶收敛。此时节点数目较多。
变步长梯形求积公式C语言_复合梯形求积公式:算例
例子三
被积函数选择为
其中
和
是参数。精确的积分值为
考虑
,
的参数情形。
之前对于这类无穷区间上的积分,我们采用的是先截断到固定区间上的做法。 此处,我们将数值积分的参数设置为步长
以及节点数
。 那么之前固定区间就等价于
。
下面是复合梯形公式的误差结果关于
的收敛,左右分别是不同的截断
的选取方法。
有以下两点观察:
- 误差在 达到
变步长梯形求积公式C语言_复合梯形求积公式:算例 左右剧烈减小变步长梯形求积公式C语言_复合梯形求积公式:算例 - 和
变步长梯形求积公式C语言_复合梯形求积公式:算例 没有明显区别变步长梯形求积公式C语言_复合梯形求积公式:算例
第一点可以通过分析
的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
函数不需要我们费心。但是如果计算资源受限,围道积分是一种选择。