simpson法求积分.doc

主程序为：
clear;clc;
n=4;a=0;b=1;
syms x t
f=(4*x^3+5*x^2-2*x+5)/(8*(x+1)^2);
g=1/(1+t);sou=1/(x+1)^2;
h=(b-a)/n;l=2*n+1;xi=a:h/2:b;
c=repmat([2,4],1,n-1);c=[1,c,2,1]; %Simpson插值的各项系数
u=c.*subs(g,xi);u=repmat(u,l,1);
v=xi'*c;
u=h*(u-v)/6-eye(l);
w=-(subs(f,xi))';
[s,yi]=columngauss(u,w); %调用列主元Gauss消去法解方程组
if ~s
disp('Error!There is something wrong.');
else
disp('When x=');disp(xi);
disp('The discrete solution of y(x) is:');disp(yi');
end
sp=myspline(xi,yi'); %调用三次样条插值函数
xk=0:0.001:1;
sk=subs(sou,xk);
for i=1:length(xk)
r=find(xi(2:l)>=xk(i),1,'first');
yk(i)=subs(sp(r),xk(i)); %计算S(xk)
end
es=max(abs(yk-sk));
disp(['The maximum error is: ',num2str(es)]);
plot(xk,sk,'r',xk,yk); %画图
title('Solve equation by numerical intergration');
legend('Acurate solution','Approximate solution');
xlabel('x');ylabel('y(x)');
grid on;