Simpson公式与复合Simpson的MATLAB实现

function [S1,S2]=Simpson(f,a,b,n)
%f为插值函数
%a为积分下限
%b为积分上限
%n公式阶数
format long
%%%%%复合Simpson%%%%%
fa=feval(f,a);
fb=feval(f,b);
h=(b-a)/n;
for i=1:n-1
    x1(i)=a+i*h;
    x2(i)=x1(i)+0.5*h;
end
fc=feval(f,a+0.5*h);
f1=feval(f,x1);
f2=feval(f,x2);
%%%%%Simpson%%%%%
S1=(h/6)*(fa+fb+2*sum(f1)+4*(fc+sum(f2)));
fd=feval(f,(a+b)/2);
S2=(b-a)*(fa+fb+4*fd)/6;
end

 来看例题

[S1,S2]=Simpson(inline('exp(x)'),1.1,1.5,4)

S1 =

   1.477523097679250

S2 =

   1.477536117650764
此程序在应用时，若积分中有瑕点，责瑕点可用很小的数来代替，如1/sinx，积分下限0，0为其瑕点，用eps表示很小的数。