龙贝格求积

matlab代码：

function R=romberg(f, a, b, n,e)
R=zeros(1,4);
R(1, 1)= (b-a)/2* (feval(f,a)+feval(f,b));
for k=1:n-1 
    h=(b-a)/2^k;
    sum=0;
    for i=1:2^ (k-1)
        sum= sum+feval(f,a+(2*i-1)*h) ;
    end
    R(k+1, 1)=R(k, 1)/2+h* sum; 
end 
for j=1:3
    fac=1/(4^j-1);
    for k=1:n-j
        R(k+j, j+1)= (1+fac) *R(k+j, j)-fac*R(k+j-1,j);
    end
end
format long
format compact
e1=abs(R(n,4)-R(n-1,4)),
if e1>e
    fprintf('对分次数不够，达不到精度\n')
    return;
else
    fprintf ('对分次数能够达到给定精度\n')
    I=R(n, 4) ;
end
fprintf ('龙贝格积分过程如下: \n')
end

function y=f(x)
y=2.*exp(-x)./pi.^0.5;
end

运行结果：