matlab二重积分

一使用两次一重积分%%二重积分f=@(x,y)exp(sin(x))*ln(y)，y从5*x积分到x^2，x从10积分到201(7.X后版本才有此函数quad2d)y1=quad2d(@(x,y)exp(sin(x)).*log(y),10,20,@(x)5*x,@(x)x.^2)2y2=quadl(@(x)arrayfun(@(x)quadl(@(y)exp(sin(x)).*log(y),5*x,x.^2),x),10,20)3y3=dblquad(@(x,y)exp(sin(x)).*log(y).*(y=5*x&y=x.^2),10,20,50,400)详细请看吴鹏老师的文章二使用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：functionf=fxy(x,y)globalki;ki=ki+1;%ki用于统计被积函数的调用次数f=exp(-x.^2/2).*sin(x.^2+y);(2)调用dblquad函数求解。globalki;ki=0;I=dblquad('fxy',-2,2,-1,1)kiI=1.57449318974494ki=1038来源精通MATLAB科学计算一书王正林,龚纯,何倩编写,电子工业出版社三复合辛普森公式(矩形积分区域)functionq=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;fori=0:n-1forj=0:m-1x=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});endendendq=((B-b)*(A-a)/36/m/n)*q;四复合梯形公式(矩形积分区域)functionq=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})+...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矩阵endF=zeros(n+1,m+1);q=0;fori=0:nforj=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+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,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=[]，则使用缺省精度和算法。例子参考=764&page=1f=inline('exp(-x.^2/2).*sin(x.^2+y)','x','y');%当然可以使用匿名函数或者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.5368603818functionsanchongjifen%%三重积分积分限可以是函数%模板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积分上限值)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)看下这个就可以了，祝好运