romberg.m

数值分析各种公式的解法

% Romberg(龙贝格)求积公式
% Romberg求积公式是由对近似值进行修正而得到更近似的公式，它能自动改变积分步长，...
% 以使其相邻两次值的绝对误差或相对误差小于预先设定的允许误差。
function s= romberg(a,b,eps)
% Romberg求积法进行数值积分，其中a与b为积分区间，eps为允许的误差
if nargin==2
    eps=1.0e-6;
elseif nargin<2
    error
    return
end
t1=10000;
t2=-10000;
n=2;
t(1,1)=0.5*(b-a)*(ff1(a)+ff1(b));
while abs(t2-t1)>=eps
    area=0.0;
    h=(b-a)/2^(n-1);
    for i=1:(2^(n-1))
        area=area+0.5*h*(ff1(a+h*(i-1))+ff1(h*i+a));
    end
    t(n,1)=area;
    for j=2:n
        for i=1:(n-j+1)
            t(i,j)=(4^(j-1)*t(i+1,j-1)-t(i,j-1))/(4^(j-1)-1);
        end
    end
    t1=t(1,n)
    t2=t(1,n-1)
    n=n+1;
end
s=t1;
return;

% 首先编写ff1(x)的函数程序
function f=ff1(x)
% f=cos(pi*x/2);
f=exp(-x^2/2)/((2*pi)^(1/2));