dblsimpson.m

来自「这些主要是讲解数值积分算法的...通过MATLAB编程来实现这些算法」· M 代码 · 共 37 行

function q=DblSimpson(f,a,A,b,B,m,n)
if(m==1 && n==1)           %辛普森公式
    q=((B-b)*(A-a)/9)*(subs(sym(f),findsym(sym(f)),{a,b})+...
        subs(sym(f),findsym(sym(f)),{a,B})+...
        subs(sym(f),findsym(sym(f)),{A,b})+...
        subs(sym(f),findsym(sym(f)),{A,B})+...
        4*subs(sym(f),findsym(sym(f)),{(A-a)/2,b})+...
        4*subs(sym(f),findsym(sym(f)),{(A-a)/2,B})+...
        4*subs(sym(f),findsym(sym(f)),{a,(B-b)/2})+...
        4*subs(sym(f),findsym(sym(f)),{A,(B-b)/2})+...
        16*subs(sym(f),findsym(sym(f)),{(A-a)/2,(B-b)/2}));
    
else                       %复合辛普森公式
    q=0;
    for i=0:n-1
        for j=0:m-1
            x=a+2*i*(A-a)/2/n;
            y=b+2*j*(B-b)/2/m;            
            x1=a+(2*i+1)*(A-a)/2/n;
            y1=b+(2*j+1)*(B-b)/2/m;
            x2=a+2*(i+1)*(A-a)/2/n;
            y2=b+2*(j+1)*(B-b)/2/m;
            
            q=q+subs(sym(f),findsym(sym(f)),{x,y})+...
                subs(sym(f),findsym(sym(f)),{x,y2})+...
                subs(sym(f),findsym(sym(f)),{x2,y})+...
                subs(sym(f),findsym(sym(f)),{x2,y2})+...
                4*subs(sym(f),findsym(sym(f)),{x,y1})+...
                4*subs(sym(f),findsym(sym(f)),{x2,y1})+...
                4*subs(sym(f),findsym(sym(f)),{x1,y})+...
                4*subs(sym(f),findsym(sym(f)),{x1,y2})+...
                16*subs(sym(f),findsym(sym(f)),{x1,y1});
        end
   end
end

q=((B-b)*(A-a)/36/m/n)*q;