正面教材分享 http://blog.sciencenet.cn/u/wdlang 70%的以色列人是无神论者,不过他们都相信上帝给了他们那块土地。这个世界经不起思考

博文

梯形法误差

已有 1920 次阅读 2019-11-19 19:50 |个人分类:计算方法|系统分类:教学心得

用梯形法计算函数有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的平方减小。

untitled11.jpg

但是如果积分区间是(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),'-*')

untitled11.jpg

我们看到误差随步长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),'-*')

untitled11.jpg




http://blog.sciencenet.cn/blog-100379-1206773.html

上一篇:pi by monte carlo
下一篇:对称或者厄米矩阵的对角化

2 尤明庆 杨正瓴

该博文允许注册用户评论 请点击登录 评论 (0 个评论)

数据加载中...

Archiver|手机版|科学网 ( 京ICP备14006957 )

GMT+8, 2019-12-9 13:51

Powered by ScienceNet.cn

Copyright © 2007- 中国科学报社

返回顶部