php画梯形,科学网—梯形法误差 - 张江敏的博文

用梯形法计算函数有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),'-*')

转载本文请联系原作者获取授权，同时请注明本文来自张江敏科学网博客。

链接地址：http://blog.sciencenet.cn/blog-100379-1206773.html

上一篇：pi by monte carlo

下一篇：对称或者厄米矩阵的对角化