romberg.m

来自「数值分析中的插值计算」· M 代码 · 共 33 行

function Romberg
%龙贝格积分算法

syms x;
f=input('被积函数式f(x)=');
a=input('积分下限a=');
b=input('积分下限b=');
yita=input('误差精度=');

h=b-a;
R=zeros(1,1);
start=subs(f,x,a);
bottom=subs(f,x,b);
R(1,1)=(b-a)*(start+bottom)/2;
k=0;
d=1;
while d>yita
    k=k+1
    sum=0;
    h=h/2;
    for i=1:2^k
        sum=sum+subs(f,x,a+(2*i-1)*h)
    end
    R(k+1,1)=sum*h+R(k,1)/2;
    for m=1:k
        R(k+1,m+1)=R(k+1,m)+[R(k+1,m)-R(k,m)]/(4^m-1);
    end
    d=abs(R(k+1,m+1)-R(k,m));
    end
end
R
d