龙贝格积分MATLAB

% 龙贝格一维/二维积分公式。积分值=rombeg（‘函数‘XminXmax[YminYmax][精度隐含精度1e-8]）
% 二维情况下可做变上下限积分（必须是匿名函数形式）.此程序积分区间和函数可为矩阵
% 但两者不能同时为矩阵。变上下限积分时只能y变量是矩阵计算效果等同x有矩阵情况
function I=rombeg（fXminXmaxYminYmaxtol）
if ~isa（f‘function_handle‘）&&~isa（f‘inline‘）
   f=inline（f）;
end
switch nargin
    case 3 %一维积分
         I=jf1（fXminXmax）; %龙贝格一维积分计算函数
    case 4%龙贝格一维积分计算函数
        I=jf1（fXminXmaxYmin）;
    case 5  %二维积分
    h=Xmax-Xmin;
    F1=@（y） f（Xminy）;
    F2=@（y） f（Xmaxy）;
     if isa（Ymin‘function_handle‘）
       c1=Ymin（Xmin）;
       c2=Ymin（Xmax）;
     else c1=Ymin;c2=Ymin;
     end
     if isa（Ymax‘function_handle‘）  
         d1=Ymax（Xmin）;
         d2=Ymax（Xmax）;
     else d1=Ymax;d2=Ymax;
     end
    s1=jf1（F1c1d1）;
    s2=jf1（F2c2d2）;
    V（::1）=（s1+s2）.*h/2;
    n=1;
    k=1;
    while k<=10
        n=2*n;
        h=h/2;
        k=k+1;
        sum=zeros（size（s1））;
        for p=1:n/2
            F=@（y） f（Xmin+（2*p-1）*hy）;
             if isa（Ymin‘function_handle‘）
            c=Ymin（Xmin+（2*p-1）*h）;
             else c=Ymin;
             end
             if isa（Ymax‘function_handle‘）  
            d=Ymax（Xmin+（2*p-1）*h）;
             else d=Ymax;
             end
            s=jf1（Fcd）;
            sum=sum+s;
        end
        V（::k）=V（::k-1）/2+（sum.*h）;
        m=k;
          while 1%龙贝格递推
            V1（::m-1）=（4*V（::m）-V（::m-1））/3;
            m=m-1;
             if m==1
              break
             else V2（::m-1）=（8*V1（::m）-V1（::m-1））/7;
              m=m-1;
             end
             if m==1
                 break
               else V3（::m-1）=（16*V2（::m）-V2（::m-1））/15;
                     m=m-1;
                     break
             end
          end
              if k==4&&（max（（max（abs（real（V3（::1）-V2（::2）））））））<=1e-8&&（max（（max（abs（imag（V3（::1）-V2（::2）））））））<=1e-8%精度控制
                 I=V3（::m）;
                break
             elseif k>4&&（max（max（abs（real（V3（::m）-V3（::m-1））））））<=1e-8&&（max（（max（abs（imag（V3（::1）-V2（::2）））））））<=1e-8
                 I=V3（::m）;
              break
              end
    end
case 6
    h=Xmax-Xmin;
    F1=@（y） f（Xminy）;
    F2=@（y） f（Xmaxy）;
     if isa（Ymin‘function_handle‘）
       c1=Ymin（Xmin）;
       c2=Ymin（Xmax）;
     else c1=Ymin;c2=Ymin;
     end
     if isa（Ymax‘function_handle‘）  
         d1=Ymax（Xmin）;
         d2=Ymax（Xmax）;
     else d1=Ymax;d2=Ymax;
     end
    s1=jf1（F1c1d1）;
    s2=jf1（F2c2d2）;
    V（::1）=（s1+s2）.*h/2;
    n=1;
    k=1;
    while k<=10
        n=2*n;
        h=h/2;
        k=k+1;