1、一 使用兩次一重積分%二重積分f= (x,y)exp(sin(x)*ln(y)，y從5*x積分到x2，x從10積分到201 (7.X后版本才有此函數quad2d) y1=quad2d(x,y) exp(sin(x).*log(y),10,20,(x)5*x,(x)x.2)2 y2 =quadl(x) arrayfun(x) quadl(y)exp(sin(x).*log(y),5*x,x.2),x),10,20)3 y3 = dblquad(x,y)exp(sin(x).*log(y).*(y>=5*x & y<=x.2),10,20,50,400)詳細請看吳鵬老師的文章二

2、 使用dblquad函數q = dblquad(fun,xmin,xmax,ymin,ymax,tol,method)該函數求f(x,y)在a,b×c,d區域上的二重定積分。參數tol，trace的用法與函數quad完全相同。 例8-5 計算二重定積分(1) 建立一個函數文件fxy.m：function f=fxy(x,y)global ki;ki=ki+1; %ki用于統計被積函數的調用次數f=exp(-x.2/2).*sin(x.2+y);(2) 調用dblquad函數求解。global ki;ki=0;I=dblquad('fxy',-2,2,-1,1)kiI

3、= 1.57449318974494ki = 1038來源 精通MATLAB科學計算 一書王正林,龔純,何倩編寫,電子工業出版社1 / 4三 復合辛普森公式(矩形積分區域)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、.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-1for j=0:m-1x=a+2*i*(A-a)/2/n;y=b+2*j*(B-b)/2/m; x1=a+(2*i+1)

5、*(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(

6、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);endendendq=(B-b)*(A-a)/36/m/n)*q;四 復合梯形公式(矩形積分區域)function q=DblTraprl(f,a,A,b,B,m,n)if(m=1 && n=1) %梯形公式q=(B-b)*(A-a)/4)*(subs(sym(f),findsym(sym(f),a,b)+.subs(sym(f),findsym(sym(f),a,B)+.s

7、ubs(sym(f),findsym(sym(f),A,b)+.subs(sym(f),findsym(sym(f),A,B);else %復合梯形公式C=4*ones(n+1,m+1);C(1,:)=2;C(:,1)=2;C(n+1,:)=2;C(:,m+1)=2;C(1,1)=1;C(1,m+1)=1;C(n+1,1)=1;C(n+1,m+1)=1; %C矩陣endF=zeros(n+1,m+1);q=0;for i=0:nfor j=0:mx=a+i*(A-a)/n;y=b+j*(B-b)/m;F(i+1,j+1)=subs(sym(f),findsym(sym(f),x,y);q=q+

8、F(i+1,j+1)*C(i+1,j+1);endendq=(B-b)*(A-a)/4/m/n)*q;*以上為矩形區域求解*五 變量區域二重積分quad2dggen功能 任意區域上二元函數的數值積分格式 q = quad2dggen(fun,xlower,xupper,ymin,ymax) %在由xlower,xupper, ymin,ymax指定的區域上計算二元函數z=f(x,y)的二重積分。q = dblquad(fun,xlower,xupper,ymin,ymax,tol) %用指定的精度tol代替缺省精度10-6，再進行計算。q = dblquad(fun,xmin,xmax,ymi

9、n,ymax,tol,method) %用指定的算法method代替缺省算法。method的取值有缺省算法或用戶指定的、與缺省命令有相同調用次序的函數句柄。q=dblquad(fun,xlower,xupper,ymin,ymax,tol,method,p1,p2,) %將可選參數 p1,p2,.等傳遞給函數fun(x,y,p1,p2,)。若tol=，method=，則使用缺省精度和算法。例子參考 f = inline('exp(-x.2/2).*sin(x.2+y)','x','y');%當然可以使用匿名函數或者M函數xlower = inli

10、ne('-sqrt(1-y.2)','y'); xupper = inline('sqrt(1-y.2)','y');Q = quad2dggen(fun,xlower,xupper,-1,1,1e-4)Q = 0.5368603818 function sanchongjifen % %三重積分積分限可以是函數 % 模板1：quadl(x) arrayfun(xx) quad2d(被積函數f(xx,y,z)關于y,z變量的函數句柄, %y積分下限y1(xx),y積分上限y2(xx),z積分下限z1(xx,y),z積分上限z2(x

11、x,y),x),x積分下限值,x積分上限值) quadl(x) arrayfun(xx) quad2d(y,z) xx.*y.*z,xx,2*xx,(y) xx*y,(y) 2*xx*y),x),1,2) % 模板2：quad2d(x,y) arrayfun(xx,yy) quadl(被積函數f(xx,yy,z)關于z變量的函數句柄,z積分下限z1(xx,yy), %z積分上限z2(xx,yy),x,y),x積分下限值,x積分上限值,y積分下限y1(x),y積分上限y2(x) quad2d(x,y) arrayfun(xx,yy) quadl(z) xx.*yy.*z,xx*yy,2*xx*yy),x,y),1,2,(x)x,(x)2*x) % 模板3：quadl(x) arrayfun(xx) quadl(y) arrayfun(yy) %quadl(被積函數f(xx,yy,z)關于z變量的函數句柄,z積分下限z1