《matlab二重积分.docx》由会员分享，可在线阅读，更多相关《matlab二重积分.docx（6页珍藏版）》请在三一办公上搜索。
1、matlab二重积分一 使用两次一重积分 %二重积分f= (x,y)exp(sin(x)*ln(y)，y从5*x积分到x2，x从10积分到20 1 (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,bc,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) ki I
3、= 1.57449318974494 ki = 1038 来源 精通MATLAB科学计算 一书王正林,龚纯,何倩编写,电子工业出版社 三 复合辛普森公式(矩形积分区域) 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*s
4、ubs(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+(
5、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),
6、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; 四 复合梯形公式(矩形积分区域) 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),find
7、sym(sym(f),a,B)+. subs(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矩阵 end F=zeros(n+1,m+1); q=0; for i=0:n for j=0:m x=a+i*(A-a)/n; y=b+j*(B-b)/m; F(i+1,j+1
8、)=subs(sym(f),findsym(sym(f),x,y); q=q+F(i+1,j+1)*C(i+1,j+1); end end q=(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) %用指定的精度t
9、ol代替缺省精度10-6，再进行计算。 q = dblquad(fun,xmin,xmax,ymin,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);%当然可以使用匿
10、名函数或者M函数 xlower = inline(-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(xx,y),x),x积分下限值,x
11、积分上限值) 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(xx,yy), %z积分上限z2(xx,yy),y),y积分下限y1(xx),y积分上限y2(xx),x),x积分下限值,x积分上限值) quadl(x) arrayfun(xx) quadl(y) arrayfun(yy) quadl(z) xx.*yy.*z,xx*yy,2*xx*yy),y),xx,2*xx),x),1,2) 看下这个就可以了，祝好运
