龙贝格算法 MATLAB实现

提示：扫一扫查出行
【扫一扫了解最新限行尾号】
复制提示

龙贝格算法主要是不断递推和加速，直到满足精度要求

递推：

加速：

得到T表：

MATLAB代码：

function I = Romberg(f, a, b, epsilon)
I = 0;
h = b-a;
k = 0;
m = 0;
T = zeros(5);
%下标转换：T^(k)_0 => T(k+1,1)
T(1,1) = h/2 * (subs(f,a) + subs(f,b));%即T^(0)_0
delta = 2*epsilon;
while delta > epsilon
k = k + 1;
h = h / 2;
tmpSum = 0; j = a + h;
while j < b
tmpSum = tmpSum + subs(f, j);
j = j + 2*h;
end
T(k+1, 1) = 0.5*T(k, 1) + h*tmpSum;
j = 2;
while j <= k+1
T(k+1,j) = (4^(j-1)*T(k+1,j-1) - T(k, j-1)) / (4^(j-1)-1);
j = j + 1;
end
delta = abs(T(k+1, k+1) - T(k, k));
end
I = T(k+1,k+1);
T
end