【原创】一般区域二重、三重积分MATLAB计算方法

rocwoods

这里讨论的计算方法指的是利用现有的MATLAB函数来求解，而不是根据具体的数值计算方法来编写相应程序。目前最新版的2009a有关于一般区域二重积分的计算函数quad2d（详细
介绍见http://forum.simwe.com/viewthread.php?tid=873479），但没有一般区域三重积分的计算函数，而NIT工具箱似乎也没有一般区域三重积分的计算函数。
本贴的目的是介绍一种在7.X版本MATLAB（不一定是2009a）里求解一般区域二重三重积分的思路方法。需要说明的是，上述链接里已经讨论了一种求解一般区域二重三重积分的思
路方法，就是将被积函数“延拓”到矩形或者长方体区域，但是这种方法不可避免引入很多乘0运算浪费时间。因此，新的思路将避免这些。由于是调用已有的MATLAB函数求解，在
求一般区域二重积分时，效率和2009a的quad2d相比有一些差距，但是相对于"延拓"函数的做法，效率大大提高了。下面结合一些简单例子说明下计算方法。
譬如二元函数f(x,y) = x*y，y从sin(x)积分到cos(x)，x从1积分到2，这个积分可以很容易用符号积分算出结果

1. 
   
2. syms x y
   
3. 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. 
   
2. function IntDemo
   
3. function f1 = myfun1(x)
   
4.      f1 = zeros(size(x));
   
5.      for k = 1:length(x)
   
6.          f1(k) = quadl(@(y) x(k)*y,sin(x(k)),cos(x(k)));
   
7.      end
   
8.      end
   
9. y = quadl(@myfun1,1,2)
   
10. end
    

复制代码

myfun1函数就是构造的原始被积函数对y积分后的函数，这时候是关于
x的函数，要能接受向量形式的x输入，所以构造这个函数的时候考虑到x是向量的情况。
利用arrayfun函数（7.X后的版本都有这个函数，不了解这个函数的朋友可以查看帮助文档，或者搜索本版）可以将IntDemo函数精简成一句代码：

1. 
   
2. 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积分到x^2，x从10积分到20。
用quad2d,上面介绍的方法，还有dblquad帮助文档里给的延拓函数的方法

1. tic,y1 = quad2d(@(x,y) exp(sin(x)).*log(y),10,20,@(x)5*x,@(x)x.^2),toc
   
2. tic,y2 = quadl(@(x) arrayfun(@(x) quadl(@(y) exp(sin(x)).*log(y),5*x,x.^2),x),10,20),toc
   
3. tic,y3 = dblquad(@(x,y) exp(sin(x)).*log(y).*(y>=5*x & y<=x.^2),10,20,50,400),toc
   
4. y1 =
   
5.     9.368671342614414e+003
   
6. Elapsed time is 0.021152 seconds.
   
7. y2 =
   
8.     9.368671342161189e+003
   
9. Elapsed time is 0.276614 seconds.
   
10. y3 =
    
11.     9.368671498376889e+003
    
12. Elapsed time is 1.674544 seconds.
    

复制代码

可见上述方法在2009a以外的版本中不失为一种方法，起码效率高于dblquad帮助文档里推荐的做法。更重要的是，这给我们求解一般区域三重积分提供了一种途径。
由于时间太晚了，求一般区域三重积分的方法留待下来有时间再写。

TBE_Legend

1# rocwoods

说一点小事。

这种帖子我很喜欢，但是不便保存，因为拷贝文字时会有麻烦。
所以，建议rocwoods版主能够顺便给出一个word或pdf格式的文档放到附件中，保持就方便得多了。当然，论坛还是以帖子的形式讨论为主，但是可以通过在附件中加上较高的仿真币限制来实现这个。

吹毛求疵了，呵呵，谢谢！

rocwoods

前面讨论的一般二重积分的例子：

1. quadl(@(x) arrayfun(@(xx) quadl(@(y) xx*y,sin(xx),cos(xx)),x),1,2)

复制代码

写成模板即：

1. 
   
2. 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. 
   
2. 模板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积分上限值)
   
3. 实例：quadl(@(x) arrayfun(@(xx) quad2d(@(y,z) xx.*y.*z,xx,2*xx,@(y) xx*y,@(y) 2*xx*y),x),1,2)
   
4. 模板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))
   
5. 实例：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)
   
6. 模板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积分上限值)
   
7. 实例：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. 
   
2. tic;quadl(@(x) arrayfun(@(xx) quad2d(@(y,z) xx.*y.*z,xx,2*xx,@(y) xx*y,@(y) 2*xx*y),x),1,2),toc
   
3. 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),toc
   
4. 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),toc
   
5. tic;triplequad(@(x,y,z) x.*y.*z.*(z<=2*x.*y&z>=x.*y&y<=2*x&y>=x),1,2,1,4,1,16),toc
   
6. ans =
   
7.   179.2969
   
8. Elapsed time is 0.037453 seconds.
   
9. ans =
   
10.   179.2969
    
11. Elapsed time is 0.223533 seconds.
    
12. ans =
    
13.   179.2969
    
14. Elapsed time is 0.090477 seconds.
    
15. ans =
    
16.   178.9301
    
17. Elapsed time is 78.421721 seconds.

复制代码

可见如果大家用的是2009a，那么计算一般区域三重积分时可以用模板1，2009a以前的版本（当然都是7.X版本）可以用模板3，而模板2效率比较低（似乎是更加频繁调用函数增加

了系统开销）。

以上二重三重积分模板还可以推广应用范围，譬如计算int(1/int(y,x,x^2),10,100)就不能套用dblquad或者quad2d，但是用一般二重积分模板稍作变形，可以这样求解：

1. 
   
2. quadl(@(x) 1./arrayfun(@(xx)quadl(@(y)y,xx,xx^2),x),10,100)
   

复制代码

PS：收藏帖子的时候可以保存成mht格式的，这样一来方便知道出处，二来将来看的时候还可以顺便看看是否有后续讨论。还有代码都写成代码格式了，这样方便copy，所以暂不上传pdf格式的了。

bainhome

建议rocwoods版主能够顺便给出一个word或pdf格式的文档放到附件

How about this?

不建议仿真币的限制。
最后，建议本贴加二级精华，个人认为MATLAB只是个工具，既然是工具，那么但凡做出了论坛前人未做之事，就是对软件使用、论坛总信息库的有益补充，那么就可以视条件加入精华，这是对作者在论坛之内最高的荣誉奖励了，精华贴的数目现在是唯一可以表达一个ID含金量的尺度了，真诚希望它再不要像积分一样被滥用了。以后加入精华建议有一个评定回帖的方式，就像版主申请通过一样，要多少多少人认同并回帖再加精华，大家看如何？
最后，如果有时间+帖子我自己看着爽，我会尽量把它们用LaTeX变成pdf格式附加上去，但是希望有更多的人参与这个工作，我认为这样做是十分有意义的，非常便于成型文章的保存，只要最后加上水印就可以了。

bainhome

另外，如下这段代码：

1. function IntDemo
   
2. 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.      end
   
7.      end
   
8. y = quadl(@myfun1,1,2)
   
9. end

复制代码

是嵌套的nested function吗？不用行不行？

messenger

4# bainhome

觉得加点限制也好，比如有些人在这边免费下载以后，转手拿这个去其它论坛赚钱，这种不劳而获总让人心理不舒服。

觉得投票就没必要了，现在本来也没几篇精华。

另外，想问一下，一级精华、二级精华、三级精华倒底有什么区别呀？点精华链接时，不是不能区分几级精华吗？

如果精华可以区分级别，可以二级及以上精华大家投票，一级精华大家看着办就行了，不用太严格。

bainhome

无所谓，反正以大家能学到东西为主要目的，作者真正珍惜的知识一般也不会在论坛上放，既然都是基础问题、典型问题的解决，知识无国界，大家都学到也不是坏事，转贴注明是谁原创即可。至于有人拿着去赚什么其他论坛的币，随他们吧，也防不了，像百思还是什么论坛(记不清了)大批量整个copy盗窃其他论坛的帖子，同时不加任何说明，一个论坛都能无耻至斯，除了鄙视之外我看最好的办法就是干脆不要写，它就抄不了。上回还看有个帖子里整篇抄junziyang关于MATLAB认不到编译器的文章，一句不提junziyang，还要junziyang自己去半开玩笑质问那位老兄，就这样还嘴硬说是只“参考”了一些，才勉强“感谢”了一下junziyang，这样的同学，你还能怎样——树不要皮，必死无疑，人不要皮，天下无敌:lol所以做好自己那部分求个心安即可。

至于精华贴，我的意思不是投票，而是得有个别几个人，他们中的多数人短信或其他方式认同精华再加精，投票是最扯淡的事情，我看论坛搞了那么多的投票没什么大用处，总之你们版主说了算，我一个打酱油的，随口一说，你们参考，hoho...
ps：题已经偏了，到此为止吧，大家还是集中在讨论积分的问题上，用脚后跟也能看出：这个帖子是非常有价值的。

rocwoods

非常感谢bainhome兄花时间整理这个帖子并用latex排版！下面代码确实是nested function结构

1. 
   
2. function IntDemo
   
3. function f1 = myfun1(x)
   
4.      f1 = zeros(size(x));
   
5.      for k = 1:length(x)
   
6.          f1(k) = quadl(@(y) x(k)*y,sin(x(k)),cos(x(k)));
   
7.      end
   
8.      end
   
9. y = quadl(@myfun1,1,2)
   
10. end
    

复制代码

不用nested，用子函数也可以。如下：

1. 
   
2. function IntDemo
   
3. y = quadl(@myfun1,1,2)
   
4. function f1 = myfun1(x)
   
5.      f1 = zeros(size(x));
   
6.      for k = 1:length(x)
   
7.          f1(k) = quadl(@(y) x(k)*y,sin(x(k)),cos(x(k)));
   
8.      end
   

复制代码

用nested function主要考虑就是对一些复杂的问题传递参数方便。而且Mathworks也推荐传递参数的时候用nested function结构。子函数同样也可以达到目的。
其实对于第一重积分有解析原函数的一般区域二重、三重积分，也可以先对第一重积分进行符号积分，得到表达式后再构造函数句柄进行数值积分。这个可以和上面的方法权衡使用。当然对于连第一重积分都没有解析原函数的积分问题，符号积分就行不通了。
被积函数表达式简单的可以用quadl+arrayfun来方便实现，对于复杂表达式的就得借助nested function或者子函数来解决了。

TBE_Legend

How about this?
206453
至于仿真币什么的，说实话别的地方我不知道，反正在仿真MATLAB版我现在最不喜欢就是那些在资料上加限制的帖子(尽管好像我以前在其他论坛也干过几次:lol)，这就好比妓女从良，花了大价钱在 ...
bainhome 发表于 2009-6-17 11:32

这个pdf做得漂亮，专业，看起来就和读纸版的书一样，很舒服。多谢bainhome。

可能是每天读文献读的，养成的看pdf的毛病。呵呵

yiyeguzhou206

这个比较精辟了.
但是不知道2009里面解决符号二重积分问题是不是有很好的办法。

sjtleo

赞一个
做学问先做人
只有诚信而务实的环境才能让众人的力量得到最大的发挥

oyanglove5212

这个帖子很使用，写的很好，收益匪浅。bainhome版主的精彩见解以及对无私的精神，让我很受感动，如果在各论坛，在学术界都有像bainhome版主一样的想法，我们国家的教育及学术会更上一层楼。

lutheo

好文章啊，看了大有启发。

ChaChing

转贴完才想到还未经rocwoods的同意就转载了, 太失礼了!
http://www.chinavib.com/forum/vi ... =page%3D1#pid460477
请rocwoods谅解, 并请跟帖好给应有的奖励!

5623728

感谢版主的无私的精神

shookerpear

受益匪浅，向rocwoods和bainhome学习。

upcjacob

真是好东西，谢谢楼主

qiaobonba

谢谢了 我太喜欢这个帖子了 希望我以后也能帮助更多像学东西的人！

zczyxmj

从微博看到的，好帖加好的推广才能让更多人知道

LuckyAngel

很好的资料，更感谢几位无私的精神，学习了、