simpson法求积分

求积分方程的离散解并比较它与真解之间的误差。用n?4的复化Simpson公式求积分方程4x3?5x2?2x?5?8(x?1)21?(01?x)y(t)dt?y(x)1?t1比较。2(x?1)

的近似解，然后用三次样条函数近似y(x)的离散近似解，并在区间?0,1?的一些点上与真解

y(x)?主程序为： 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;

本文来源：https://www.bwwdw.com/article/wwix.html