rombergaccelerate.m

来自「研一刚上完数值分析」· M 代码 · 共 31 行

function[s,n]=RombergAccelerate(fun,a,b,eps)
%龙贝格加速算法
%a,b为求积上下界
%eps为精度
if nargin<4,eps=1e-6;end
s=10;
s0=0;
k=2;
t(1,1)=(b-a)*(fun(a)+fun(b))/2;
while (abs(s-s0)>eps)
    h=(b-a)/2^(k-1);
    w=0;
    if(h~=0)
        for i=1:(2^(k-1)-1)
            w=w+fun(a+i*h);
        end
        t(k,1)=h*(fun(a)/2+w+fun(b)/2);
        for l=2:k
            for i=1:(k-l+1)
                t(i,l)=(4^(l-1)*t(i+1,l-1)-t(i,l-1))/(4^(l-1)-1);
            end
        end
        s=t(1,k);
        s0=(t(1,k-1));
        k=k+1;
        n=k;
    else s=s0;
        n=-k;
    end
end