demotrapsimpsimp4.m

梯形与simpson求积公式 分别是梯形,simpson,和3/8simpson的matlab代码,计算的区间被划分成为了又2到256个,精度逐渐增加.

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

%----- Simpson4 method
fprintf('\nsimpson4:')
fprintf('\n    n       h         i             error           alpha\n');

errold = [];
for np = [2 4 8 16 32 64 128 256]
    I = simpson4('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