romberg.m

数值分析---部分算法 希望对需要的人有一点点用！

% -------------------------------------------------------------------------
% ----------------------------龙贝格积分算法--------------------------------
% 函数  In=Romberg(func,a,b)
% 输入参数: func   被积函数的名称,如'sin' 
%          a,b    积分区间
% 输出结果: In     积分结果,在命令窗口显示ans=____
%    >> I = Romberg( 'sin', 0, pi/2)
% 编程思路:1.设计一个矩阵rb中,根据算法得出矩阵元素之间的关系
%         2.积分限由(a,b)--->(0,1)  便于寻找规律进行编程
%         3.也有可能会出现计算很浅（行数），就终止的情况，这会出现很大误差！！！
%           譬如：f(0.5)=0.5*[f(0)+f(1)],第二行就停了。
%           为避免此种情况，让程序先自动运行jingdu行，再判断误差，用大计算量换取可靠。
% -------------------------------------------------------------------------
function In=Romberg(func,a,b)
wucha=1e-8;
jingdu=6;
rb(1,1)=(b-a)*(feval(func,a)+feval(func,b))/2;  
i=1;
while i<=jingdu||abs(rb(i,i)-rb(i,i-1))>wucha  %第一步判断条件时i＝1，但是执行不到abs(rb(i,i)-rb(i,i-1))>wucha!!!
    i=i+1;
    temp=0;
    for k=1:2:2^(i-1)
        temp=temp+(b-a)*feval(func,a+(b-a)*k/(2^(i-1)));
    end
    rb(i,1)=rb(i-1,1)/2+temp/(2^(i-1));%T
    for j=2:i
        rb(i,j)=(-rb(i-1,j-1)+(4^(j-1))*rb(i,j-1))/(4^(j-1)-1);%S T C R
    end
end
In=rb(i,i);