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

ChaChing

转载自http://forum.simwe.com/thread-885049-1-1.html

这里讨论的计算方法指的是利用现有的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. syms x y
   
2. int(int(x*y,y,sin(x),cos(x)),1,2) ]
   
3. 结果是 -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 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

复制代码

myfun1函数就是构造的原始被积函数对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积分到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帮助文档里推荐的做法。更重要的是，这给我们求解一般区域三重积分提供了一种途径。
由于时间太晚了，求一般区域三重积分的方法留待下来有时间再写。

[ 本帖最后由 ChaChing 于 2010-3-2 10:10 编辑 ]

ChaChing

转贴完才想到还未经rocwoods的同意就转载了, 太失礼了!
请rocwoods谅解, 并请跟帖好给应有的奖励!

rocwoods

呵呵，Chaching兄太客气了。我这半年多一直潜水写书，现在接近尾声了，书的内容就是有关MATLAB高效编程以及一些计算技巧的书，里面收录了这些年自己帮助别人以及自己研究的一些内容，详细讨论了数值积分、Fredholm、Volterra积分方程、各类型的微分方程等解法，也有一些优化方面的、遗传算法、非线性方程求解、支持向量机、层次分析法等案例！
上述方法我也在书籍中予以详细讨论，并将该方法推广到了任意区域n重积分情况，并有详细源代码。
书籍预计6月前后出版，欢迎大家关注。
最近等书籍定稿了，我会把目录发上来。

[ 本帖最后由 rocwoods 于 2010-3-2 11:36 编辑 ]

ChaChing

希望有网购服务, 届时再拜读学习
期待中!

maigicku

书名是什么啊？希望机会拜读。。。

xiezhh

恭喜rocwoods 的书即将完成！期待早日拜读！

rocwoods

感谢大家的支持，书名暂定为：《MATLAB高效编程技巧与应用：25个案例分析》
由北航出版社出版

[ 本帖最后由 rocwoods 于 2010-3-3 10:56 编辑 ]

yudailin

学习到了很多东西，牛人！！
期待拜读你的著作~

ChaChing

刚才发现早请rocwoods转帖过来了
http://forum.vibunion.com/forum/ ... 3426&highlight=
个人记性真差

liwang11

希望有机会拜读。。。

ChaChing

原帖由 liwang11 于 2010-4-4 21:21 发表

                               
登录/注册后可看大图

希望有机会拜读。。。

不是贴好了!?
不能读吗?

rocwoods

呵呵，我想10楼的朋友是想等拙作出来之后阅读吧。谢谢支持！

ChaChing

喔! 没错, 又误会了!

maigicku

个人感觉还是用匿名函数比较容易。。三重积分在二重外再加一个用匿名函数也是可以的吧。。
强烈期待LZ的新书。。