Question

matlab中如何用复合辛普森公式求二重积分

Answer 1

代码如下：
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;
叫我雷锋！