romberg.m

来自「用matlab解决一些数值分析中常用的算法」· M 代码 · 共 47 行

%龙贝格积分法
function [result,matrix]=Romberg(fun_in,a,b,e)
format short;
syms x;
fun=sym(fun_in);
N=1;
T=zeros(1);
x=a;if(x==0) x=eps;end
fa=eval(fun);
x=b;if(x==0) x=eps;end  
fb=eval(fun);
T(1,1)=(b-a)/2*(fa+fb) ;% 初值  
while 1
    N=N+1;
    T0=zeros(N);
    T0(1:N-1,1:N-1)=T; 
    T=T0; 
    for j=1:1:N
      for k=N-j:1:N-j   % k为分区间的次数,每次对角线元素
          if j==1       % 梯形递推求的第一行  
              if k>0    % 2^k区间段个数
                  T(j,k+1)=T(j,k)/2; 
                  for i=0:1:2^k-1  
                      x=a+(b-a)/2^k*(i+0.5);
                      %if(x==0) x=x+eps;end  
                      T(j,k+1)=T(j,k+1)+(b-a)/2^(k+1)*eval(fun);
                  end
              end
          else          % 龙贝格公式求解其余 
            T(j,k+1)=(4^(j-1)*T(j-1,k+2)-T(j-1,k+1))/(4^(j-1)-1);            
          end           % 结束龙贝格公式求解其余
      end               % 结束一行的生成 
     
    end                 % 结束矩阵的生成
             
    if (abs(T(N,1)-T(N-1,1))<e) %精度判断
       result=T(N,1);
       break;  
    end                % 结束精度判断   
    
    if N>7 break;
    end
end                     % 矩阵扩大结束
matrix=T