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;
发布者：全栈程序员-站长，转载请注明出处：https://javaforall.net/225708.html原文链接：https://javaforall.net
