数值积分（辛普森求积、柯特斯求积、龙贝格求积）

利用复化辛普森求积公式计算${\int }_a^{b}f(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:n
    x=x+h/2;s=s+4*feval(f_name,x);
    x=x+h/2;s=s+2*feval(f_name,x);
end
s=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))>e
    i=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:i
        T(i+1,j+1)=4^j*T(i+1,j)/(4^j-1)-T(i,j)/(4^j-1);
    end
end
T
r=(T(i+1,j+1));