demotrapsimpsimp4.asv

来自「梯形与simpson求积公式 分别是梯形,simpson,和3/8simpso」· ASV 代码 · 共 43 行

function demoAll
%demoAll Use all 3 ways to integrate x*exp(-x) on [0,5]
a = 0; b = 5; Iexact = -exp(-b)*(1+b) + exp(-a)*(1+a);
fprintf('\n\tIexact = %14.9f\n', Iexact);

fprintf('\ntrapezoid:')
fprintf('\n   n    h     i      error     alpha\n');

errold = [];
for np = [2 4 8 16 32 64 128 256]
    I = trapezoid('xemx',a,b,np);
    err = I - Iexact;
    n = np + 1; h = (b-a)/(n-1);
    fprintf(' %4d %9.5d %14.9f %14.9f',n,h,I,err);
    
    if ~isempty(errold)
        fprintf(' %9.5f\n', log(err/errold)/log(h/hhold));
    else
        fprintf('\n');
    end
    hhold = h; errold = err;
end

%----- Simpson method
fprintf('\nsimpson:')
fprintf('\n   n    h     i      error     alpha\n');

errold = [];
for np = [2 4 8 16 32 64 128 256]
    I = simpson('xemx',a,b,np);
    err = I - Iexact;
    n = np + 1; h = (b-a)/(n-1);
    fprintf(' %4d %9.5d %14.9f %14.9f',n,h,I,err);
    
    if ~isempty(errold)
        fprintf(' %9.5f\n', log(err/errold)/log(h/hhold));
    else
        fprintf('\n');
    end
    hhold = h; errold = err;
end