simpson.m

% 复合辛普森公式
% 注释的地方需要修改
clear

a=0; b=pi/4;          % a为积分下界，b为积分上界
n=6;                  % 分割的区间数

syms x
f=(4-sin(x)^2)^0.5;   % 函数表达式

h=(b-a)/n;

for k=0:n
    xk(k+1)=a+k*h;
end

for k=0:n-1
    xk1(k+1)=xk(k+1)+1/2*h;
end

fa=subs(f,x,xk(1));
fb=subs(f,x,xk(n+1));

Sn=h/6*(fa+fb);

for k=1:n 
    fxk1=subs(f,x,xk1(k));
    Sn=Sn+2*h/3*fxk1;
end

for k=2:n 
    fxk=subs(f,x,xk(k));
    Sn=Sn+h/3*fxk;
end

vpa(Sn,7)             % 结果有效数字位数