梯形法误差

用梯形法计算函数有y = sqrt(x) 在区间（1,2）上的积分。

clear all; close all; clc; 

Nlist = 2.^(2:14);

slist = zeros(1, length(Nlist));

elist = zeros(1, length(Nlist)); 

for s = 1 : length(Nlist)

    N = Nlist(s);

    h = 1 /N ;

    f = sqrt((N: 2*N)./N);

    slist(s) = h * (0.5*(f(1) + f(end)) + sum(f(2: end-1)));

    elist(s) = abs(slist(s) - 2/3*(2^1.5-1));

end

h0 = figure;

plot(log2(Nlist), log2(elist),"-*")

可以看到误差随步长h的平方减小。

但是如果积分区间是（0,1）呢？

Nlist = 2.^(2:14);

slist = zeros(1, length(Nlist));

elist = zeros(1, length(Nlist)); 

for s = 1 : length(Nlist)

    N = Nlist(s);

    h = 1 /N ;

    f = sqrt((0: N)./N);

    slist(s) = h * (0.5*(f(1) + f(end)) + sum(f(2: end-1)));

    elist(s) = abs(slist(s) - 2/3);

end

h1 = figure;

plot(log2(Nlist), log2(elist),"-*")

我们看到误差随步长h的3/2次方减小。这个相对缓慢的幂次来自被积函数在区间端点x=0处的奇异性。

实际工作中遇到这种奇异性时，我们可以考虑解析地扣除之。比如考虑函数y= sinx/sqrt(x)在区间（0,1）上的积分。这个函数在x=0处有一个形如sqrt（x）的奇异部分。我们可以扣除之，将函数y分解成y1=sqrt（x）和y2=(sinx-x)/sqrt(x)之和。对y1我们可以解析地做，对y2我们可以用梯形法。这样可以重新获得梯形法的平方收敛性。

Nlist = 2.^(2:14);

slist2 = zeros(1, length(Nlist));

for s = 1 : length(Nlist)

    N = Nlist(s);

    h = 1 /N ;

    x = (1: N)./N ;

    f =  [0,(sin (x) - x)./ sqrt(x)];

    slist2(s) = h * (0.5*(f(1) + f(end)) + sum(f(2: end-1))) + 2/3 ;

end

elist2 = abs(slist2(1:end-1) - slist2(end));

h2 = figure;

plot(log2(Nlist(1: end-1)), log2(elist2),"-*")