利用复化辛普森求积公式计算∫abf(x)dx\int _ { a } ^ { b } f ( x ) d x∫abf(x)dx的近似值
辛普森求积
function s=simpson( f_name,a, b, n)%复化辛普生公式% f_name为要求的定函数y=f(x)所在的程序文件名% a为积分下限% b为积分上限% n为积分区间[a,b]划分成小区间的等份数h=(b-a)/n;fa=feval(f_name,a);fb=feval(f_name,b);s=fa-fb;x=a;for i=1:nx=x+h/2;s=s+4*feval(f_name,x);x=x+h/2;s=s+2*feval(f_name,x);ends=s*h/6;s=vpa(s,6);end
柯特斯求积
function t=trapz(fname,a,b,n)% fname为被积函数，a,b分别为下界和上界，n为等分数h=(b-a)/n;fa=feval(fname,a);fb=feval(fname,b);f=feval(fname,a+h:h:b-h+0.001*h);t=h*(0.5*(fa+fb)+sum(f));
龙贝格求积
function r=agui_rbg(fname,a,b)%fname为被积函数，a,b,分别为上界和下界e=1e-6;i=1;j=1;h=b-a;T(i,1)=h/2*(feval(fname,a)+feval(fname,b));T(i+1,1)=T(i,1)/2+sum(feval(fname,a+h/2:h:b-h/2+0.001*h))*h/2;T(i+1,j+1)=4^j*T(i+1,j)/(4^j-1)-T(i,j)/(4^j-1);while abs(T(i+1,i+1)-T(i,i))>ei=i+1;h=h/2;T(i+1,i)=T(i,1)/2+sum(feval(fname,a+h/2:b-h/2+0.001*h))*h/2;for j=1:iT(i+1,j+1)=4^j*T(i+1,j)/(4^j-1)-T(i,j)/(4^j-1);endendTr=(T(i+1,j+1));
