r_quad_iter.m

来自「matlab具有强大的数值计算能力」· M 代码 · 共 21 行

function I=R_quad_iter(fun,a,b,ep)
%Romberg求积公式，其中fun为被积函数；a,b为积分区间的端点，要求a<b;ep为精度要求。
if nargin< 4 ep = 1e-6; end;
m = 1;h = b - a;
I = h/2 * (feval(fun,a) + feval(fun,b));T(1,1) = I;
while 1
    N = 2^(m-1);h = h/2; I = I/2;
    for i = 1:N
        I = I + h * feval(fun,a + (2 * i - 1) * h);
    end
    T(m +1,1) = I; M = 2 * N; k = 1;
    while M > 1;
        T(m + 1,k + 1) = (4^k * T(m+1,k) - T(m,k))/(4^k-1);
        M = M / 2;k = k + 1;
    end
    if abs(T(k,k) - T(k - 1,k - 1)) < ep
        break;
    end
    m = m + 1;
end
I = T(k,k);