matlab实现复合梯形法则

复合梯形法则：

function int_f = CompoundEchelon( f, a, b, m )
% input :   f : function handler
%           a : the lower limit of integral
%           b : the upper limit of integral
%           m : cut integral area into m peace
% output :  int_f : the answer of the integral
    h = (b - a) / m;
    int_f = 0;
    if m >= 2
        for i = 1 : m-1
           int_f = int_f + 2 * f(a + h * i); 
        end
    end
    int_f = int_f + f(a) + f(b);
    int_f = int_f * h / 2;
end

例子：

clear all
format long
clc
%% (a)
fprintf(' (a) 
')
f = @(x) x./((x.^2+9).^0.5);
int1_16 = CompoundEchelon(f, 0, 4, 16);
int1_32 = CompoundEchelon(f, 0, 4, 32);
correct_int1 = quadgk(f, 0, 4);
error1_16 = abs(correct_int1 - int1_16);
error1_32 = abs(correct_int1 - int1_32);
fprintf('int1_16 = %g
', int1_16);
fprintf('int1_32 = %g
', int1_32);
fprintf('correct_int1 = %g
', correct_int1);
fprintf('error1_16 = %g
', error1_16);
fprintf('error1_32 = %g
', error1_32);
%% (b)
fprintf(' (b) 
')
f = @(x) (x.^3)./(x.^2+1);
int2_16 = CompoundEchelon(f, 0, 1, 16);
int2_32 = CompoundEchelon(f, 0, 1, 32);
correct_int2 = quadgk(f, 0, 1);
error2_16 = abs(correct_int2 - int2_16);
error2_32 = abs(correct_int2 - int2_32);
fprintf('int2_16 = %g
', int2_16);
fprintf('int2_32 = %g
', int2_32);
fprintf('correct_int2 = %g
', correct_int2);
fprintf('error2_16 = %g
', error2_16);
fprintf('error2_32 = %g
', error2_32);
%% (c)
fprintf(' (c) 
')
f = @(x) x.*exp(x);
int3_16 = CompoundEchelon(f, 0, 1, 16);
int3_32 = CompoundEchelon(f, 0, 1, 32);
correct_int3 = quadgk(f, 0, 1);
error3_16 = abs(correct_int3 - int3_16);
error3_32 = abs(correct_int3 - int3_32);
fprintf('int3_16 = %g
', int3_16);
fprintf('int3_32 = %g
', int3_32);
fprintf('correct_int3 = %g
', correct_int3);
fprintf('error3_16 = %g
', error3_16);
fprintf('error3_32 = %g
', error3_32);
%% (d)
fprintf(' (d) 
')
f = @(x) (x.^2).*(log(x));
int4_16 = CompoundEchelon(f, 1, 3, 16);
int4_32 = CompoundEchelon(f, 1, 3, 32);
correct_int4 = quadgk(f, 1, 3);
error4_16 = abs(correct_int4 - int4_16);
error4_32 = abs(correct_int4 - int4_32);
fprintf('int4_16 = %g
', int4_16);
fprintf('int4_32 = %g
', int4_32);
fprintf('correct_int4 = %g
', correct_int4);
fprintf('error4_16 = %g
', error4_16);
fprintf('error4_32 = %g
', error4_32);
%% (e)
fprintf(' (e) 
')
f = @(x) (x.^2).*(sin(x));
int5_16 = CompoundEchelon(f, 0, pi, 16);
int5_32 = CompoundEchelon(f, 0, pi, 32);
correct_int5 = quadgk(f, 0, pi);
error5_16 = abs(correct_int5 - int5_16);
error5_32 = abs(correct_int5 - int5_32);
fprintf('int5_16 = %g
', int5_16);
fprintf('int5_32 = %g
', int5_32);
fprintf('correct_int5 = %g
', correct_int5);
fprintf('error5_16 = %g
', error5_16);
fprintf('error5_32 = %g
', error5_32);
%% (f)
fprintf(' (f) 
')
f = @(x) (x.^3)./((x.^4-1).^0.5);
int6_16 = CompoundEchelon(f, 2, 3, 16);
int6_32 = CompoundEchelon(f, 2, 3, 32);
correct_int6 = quadgk(f, 2, 3);
error6_16 = abs(correct_int6 - int6_16);
error6_32 = abs(correct_int6 - int6_32);
fprintf('int6_16 = %g
', int6_16);
fprintf('int6_32 = %g
', int6_32);
fprintf('correct_int6 = %g
', correct_int6);
fprintf('error6_16 = %g
', error6_16);
fprintf('error6_32 = %g
', error6_32);
%% (g)
fprintf(' (g) 
')
f = @(x) 1./((x.^2+4).^0.5);
int7_16 = CompoundEchelon(f, 0, 2*3^0.5, 16);
int7_32 = CompoundEchelon(f, 0, 2*3^0.5, 32);
correct_int7 = quadgk(f, 0, 2*3^0.5);
error7_16 = abs(correct_int7 - int7_16);
error7_32 = abs(correct_int7 - int7_32);
fprintf('int7_16 = %g
', int7_16);
fprintf('int7_32 = %g
', int7_32);
fprintf('correct_int7 = %g
', correct_int7);
fprintf('error7_16 = %g
', error7_16);
fprintf('error7_32 = %g
', error7_32);
%% (h)
fprintf(' (h) 
')
f = @(x) x./((x.^4+1).^0.5);
int8_16 = CompoundEchelon(f, 0, 1, 16);
int8_32 = CompoundEchelon(f, 0, 1, 32);
correct_int8 = quadgk(f, 0, 1);
error8_16 = abs(correct_int8 - int8_16);
error8_32 = abs(correct_int8 - int8_32);
fprintf('int8_16 = %g
', int8_16);
fprintf('int8_32 = %g
', int8_32);
fprintf('correct_int8 = %g
', correct_int8);
fprintf('error8_16 = %g
', error8_16);
fprintf('error8_32 = %g
', error8_32);