一般区域二重、三重积分MATLAB计算方法.doc

一般区域二重、三重积分MATLAB计算方法这里讨论的计算方法指的是利用现有的MATLAB函数来求解，而不是根据具体的数值计算方法来编写相应程序。目前最新版的2009a有关于一般区域二重积分的计算函数quad2d（详细介绍见/viewthread.php?tid=873479），但没有一般区域三重积分的计算函数，而NIT工具箱似乎也没有一般区域三重积分的计算函数。 本文目的是介绍一种在7.X版本MATLAB（不一定是2009a）里求解一般区域二重三重积分的思路方法。需要说明的是，上述链接里已经讨论了一种求解一般区域二重三重积分的思路方法，就是将被积函数“延拓”到矩形或者长方体区域，但是这种方法不可避免引入很多乘0运算浪费时间。因此，新的思路将避免这些。由于是调用已有的MATLAB函数求解，在求一般区域二重积分时，效率和2009a的quad2d相比有一些差距，但是相对于延拓函数的做法，效率大大提高了。下面结合一些简单例子说明下计算方法。diff(x) 如果x是向量，返回向量x的差分；如果x是矩阵，则按各列作差分diff(x,k) k阶差分q=polyder(p) 求得由向量p表示的多项式导函数的向量表示qFx=gradient(F,x) 返回向量F表示的一元函数沿x方向的导函数F(x)，其中x是与F同维数的向量z=trapz(x,y) x表示积分区间的离散化向量；y是与x同维数的向量，表示被积函数；z返回积分的近似值z=guad(fun,a,b,tol) 自适应步长Simpson积分法求得Fun在区间a,b上的定积分，Fun为M文件函数句柄，tol为积分精度z=dblquad(fun,a,b,c,d,tol,method) 求得二元函数Fun(x,y)的重积分z=triplequad(fun,a,b,c,d,e,f,tol,method) 求得三元函数Fun(x,y,z)的重积分 譬如二元函数f(x,y) = x*y，y从sin(x)积分到cos(x)，x从1积分到2，这个积分可以很容易用符号积分算出结果 1. syms x y2. int(int(x*y,y,sin(x),cos(x),1,2) 结果是 -1/2*cos(1)*sin(1)-1/4*cos(1)2+cos(2)*sin(2)+1/4*cos(2)2 = -0.635412702399943如果你用的是2009a，你可以用 1. quad2d(x,y) x.*y,1,2,(x)sin(x),(x)cos(x),AbsTol,1e-12)得到上述结果。 如果用的不是2009a，那么你可以利用NIT工具箱里的quad2dggen函数。 那么我们如果既没有NIT工具箱用的也不是2009a，怎么办呢？答案是我们可以利用两次quadl函数，注意到quadl函数要求积分表达式必须写成向量化形式，所以我们构造的函数必须能接受向量输入。见如下代码 1. function IntDemo2. function f1 = myfun1(x)3. f1 = zeros(size(x);4. for k = 1:length(x)5. f1(k) = quadl(y) x(k)*y,sin(x(k),cos(x(k);6. end7. end8. y = quadl(myfun1,1,2)9. endmyfun1函数就是构造的原始被积函数对y积分后的函数，这时候是关于x的函数，要能接受向量形式的x输入，所以构造这个函数的时候考虑到x是向量的情况。利用arrayfun函数（7.X后的版本都有这个函数，不了解这个函数的朋友可以查看帮助文档，或者搜索本版）可以将IntDemo函数精简成一句代码： 1. quadl(x) arrayfun(xx) quadl(y) xx*y,sin(xx),cos(xx),x),1,2)上面这行代码体现了用MATLAB7.X求一般区域二重积分的一般方法。可以这么理解这句代码：首先 1. (x) arrayfun(xx) quadl(y) xx*y,sin(xx),cos(xx),x)定义了一个关于x的匿名函数，供quadl调用求最外重（x从1到2的）积分，这时候，x对于 1. arrayfun(xx) quadl(y) xx*y,sin(xx),cos(xx),x) 就是已知的了。而 1. (xx) quadl(y) xx*y,sin(xx),cos(xx)定义的是对于给定的xx,求xx*y关于y的积分函数，这就相当于数学上积完第一重y的积分后得到一个关于xx的函数而 1. arrayfun(xx) quadl(y) xx*y,sin(xx),cos(xx),x)只是对 1. (xx) quadl(y) xx*y,sin(xx),cos(xx)加了一个循环的壳，保证“积完第一重y的积分后得到一个关于xx的函数”能够接受向量化的xx的输入，从而能够被quadl调用。有了这个模板，我们可以方便求其他一般积分区域（上下限是函数）形式的二重积分，例如被积函数f = (x,y) exp(sin(x)*ln(y)，y从5*x积分到x2，x从10积分到20。用quad2d,上面介绍的方法，还有dblquad帮助文档里给的延拓函数的方法 1. tic,y1 = quad2d(x,y) exp(sin(x).*log(y),10,20,(x)5*x,(x)x.2),toc2. tic,y2 = quadl(x) arrayfun(x) quadl(y) exp(sin(x).*log(y),5*x,x.2),x),10,20),toc3. tic,y3 = dblquad(x,y) exp(sin(x).*log(y).*(y=5*x & y=x.2),10,20,50,400),toc4. y1 =5. 9.368671342614414e+0036. Elapsed time is 0.021152 seconds.7. y2 =8. 9.368671342161189e+0039. Elapsed time is 0.276614 seconds.10. y3 =11. 9.368671498376889e+00312. Elapsed time is 1.674544 seconds.可见上述方法在2009a以外的版本中不失为一种方法，起码效率高于dblquad帮助文档里推荐的做法。更重要的是，这给我们求解一般区域三重积分提供了一种途径。由于时间太晚了，求一般区域三重积分的方法留待下来有时间再写。 对于积分补充一点，有人经常会问对于有参数的积分怎么做。问题：y=sin(a*x)在不同的a参数情形下如何计算x从pi/4,pi/2的积分这里给出个matlab的数值积分函数的一个用法 1. a = 0.1:0.1:1;2. result = zeros(1,length(a);3. for i = 1: length(a)4. result (i)=quadl(x,a) sin(a*x),pi/4,pi/2,a(i)5. endquadl类函数在输入其自身的5个参数外，还可以将被积函数的除被积变量外的其他函数参数按照顺序放在函数参数后面。 上述用法需要注意quadl函数自身的参数需要写全，比如示例中的两个需要写上，其代表的参数含义可参考help文档。quadl(fun,a,b,tol,trace) with non-zero trace shows the values of fcnt a b-a q during the recursion. 上述用法quadl函数的help文档中没有给出，但在Anonymous Functions的介绍中有说明，包括rocwoods所说的多重匿名函数PS：另外这个问题当然也可以采用其他常规的方法，如应用符号工具箱及多重循环技术分 + 3 技术分 发短消息 加为好友 rocwoods 当前离线 帖子1183积分2865技术分150 技术分研学币1286 研学币学历研究方向职业所在地在线时间91 小时注册时间2005-9-13最后登录2009-8-30前面讨论的一般二重积分的例子： 1. quadl(x) arrayfun(xx) quadl(y) xx*y,sin(xx),cos(xx),x),1,2)写成模板即： 1. quadl(x) arrayfun(xx) quadl(y) 被积二元函数f(xx,y),y的积分下限表达式g1(xx),y的积分上限表达式g2(xx),x),x积分下限值,x积分上限值)现在来看一般区域三重积分的做法，有两种思路，一种是quadl+quad2d函数，这需要2009a来支持，另一个是用三个quadl函数。前者还可细分成先quad2d后quadl，以及先quadl后quad2d。我们可以得到三种模板，同时结合f(x,y,z) = x*y*z ;z从x*y到2*x*y积分,y从x到2*x积分，x从1到2积分这个简单的三重积分例子说明 1. 模板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积分上限值)2. 实例：quadl(x) arrayfun(xx) quad2d(y,z) xx.*y.*z,xx,2*xx,(y) xx*y,(y) 2*xx*y),x),1,2)3. 模板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)4. 实例：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)5. 模板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积分上限值)6. 实例：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)模板使用说明：x,y,z是积分变量，模板中除了用语言描述的参量用相应表达式替换掉外，其余结构保持不变。大家可以实际运行这三个实例，都比用triplequad 拓展函数法快得多，而且triplequad拓展函数得到的结果还不精确，在我的电脑上结果如下： 1. tic;quadl(x) arrayfun(xx) quad2d(y,z) xx.*y.*z,xx,2*xx,(y) xx*y,(y) 2*xx*y),x),1,2),toc2. tic;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),toc3. tic;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),toc4. tic;triplequad(x,y,z) x.*y.*z.*(z=x.*y&y=x),1,2,1,4,1,16),toc5. ans =6. 179.29697. Elapsed time is 0.037453 seconds.8. ans =9. 179.296910. Elapsed time is 0.223533 seconds.11. ans =12. 179.296913. Elapsed time is 0.090477 seconds.14. ans =15. 178.930116. Elapse