MATLAB实现Newton Cotes求积

clear
syms t x;
 
f(x)=input('请输入函数f(x)=');
a=input('请输入区间下限:  ');
b=input('请输入区间上限:  ');
n=input('请输入n（即将区间几等分）: ');
h=(b-a)/n;
 
for i=1:n+1   %将区间n等分后，计算各节点的横坐标
    xdata(i)=a+(i-1)*h;
end
ydata=subs(f,'x',xdata);   
 
for k=1:n+1    %计算Cotes系数
    l=1;
    for j=1:n+1
        if(j~=k)
            l=l*(t-j+1)/(k-j);
        end
    end
    s(k)=int(l,0,n)/n;
end
 
yy1=sum(s.*ydata);
yy2=int(f,a,b);    %直接用matlab求积
disp('精确值');
vpa(yy2,8)
 
If=(b-a)*yy1;   %用Newton_Cptes求积公式求得
disp('用Newton_Cptes求积公式求得：');
vpa(If,8)
 
wucha=yy2-If;   %计算误差
 
if (wucha<0)
    wucha=-wucha;
end
 
disp('误差为： ')
vpa(wucha,8)

将区间2等分，即为梯形公式公式，n取1

显然误差相对较大

将区间3等分，即Simpson公式，n取2