Romberg求积

使用梯形求积公式进行数值积分，虽然算法简单，但收敛慢，精度低。而Romberg求积公式可以克服这个缺点。它是对近似值进行修正以后得到的更近似的公式，能自动改变积分步长，以使其相邻的两个值的绝对误差或相对误差小于预先设定的允许误差。Romberg求积公式为

\[{T}_{j+!}^{k}f=\frac{{4}^{j}{T}_{i}^{k+1}f-{T}_{j}^{k}f}{{4}^{j}-1}\]

据此，可以编写Romberg求积公式的M文件Romberg.m，如下所示。

code.matlab 
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;
      % n=n+1;
      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

【例25】用Romberg求积公式计算定积分\(I=\int_{0}^{1} {\operatorname{cos}{(}\frac{\pi}{2}x)dx}\)。

在MATLAB中按照如下步骤实现。

首先，编写求积函数ff1.m。

code.matlab 
function y=ff1(x)
  y=cos(pi*x/2);

然后，可以利用编写的Romberg求积函数进行求积，如下所示。

code.matlab 
>> Romberg(0,1)
ans =
    0.6366