%Compound trapezoid formula复化梯形法
function y = traint(a,b,n,f)
h = (b - a) / n;
x = linspace(a,b,n+1);
y1 = h * feval(f,x);
y1(1) = y1(1) / 2;
y1(n+1) = y1(n+1) / 2;
y = sum(y1);
%compound Simpson formula复化辛普森公式
function y = sraint(a,b,n,f)
h = (b - a) / n;
x = linspace(a,b,2*n+1);
y1 = feval(f,x);
y1(2:2:2*n) = 4 * y1(2:2:2*n);
y1(3:2:2*n-1) = 2 * y1(3:2:2*n-1);
y = h / 6 * sum(y1);
>>f=inline('4./(1+x.*x)');
>> traint(-1,1,8,f)
ans =
6.2624
>> sraint(-1,1,8,f)
ans =
6.2832
