romberg.m

来自「《MATLAB实用指南》系列丛书DE源代码」· M 代码 · 共 32 行

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