一般区域二重三重积分MATLAB计算方法.docx
《一般区域二重三重积分MATLAB计算方法.docx》由会员分享，可在线阅读，更多相关《一般区域二重三重积分MATLAB计算方法.docx（12页珍藏版）》请在冰豆网上搜索。
一般区域二重三重积分MATLAB计算方法
[【原创】]一般区域二重、三重积分MATLAB计算方法
这里讨论的计算方法指的是利用现有的MATLAB函数来求解，而不是根据具体的数值计算方法来编写相应程序。
目前最新版的2009a有关于一般区域二重积分的计算函数quad2d（详细
介绍见
本贴的目的是介绍一种在7.X版本MATLAB（不一定是2009a）里求解一般区域二重三重积分的思路方法。
需要说明的是，上述链接里已经讨论了一种求解一般区域二重三重积分的思
路方法，就是将被积函数“延拓”到矩形或者长方体区域，但是这种方法不可避免引入很多乘0运算浪费时间。
因此，新的思路将避免这些。
由于是调用已有的MATLAB函数求解，在
求一般区域二重积分时，效率和2009a的quad2d相比有一些差距，但是相对于"延拓"函数的做法，效率大大提高了。
下面结合一些简单例子说明下计算方法。
譬如二元函数f（x,y）=x*y，y从sin（x）积分到cos（x），x从1积分到2，这个积分可以很容易用符号积分算出结果
1.symsxy
2.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.functionIntDemo
2.functionf1=myfun1（x）
3. f1=zeros（size（x））;
4. fork=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.Elapsedtimeis0.021152seconds.
7.y2=
8. 9.368671342161189e+003
9.Elapsedtimeis0.276614seconds.
10.y3=
11. 9.368671498376889e+003
12.Elapsedtimeis1.674544seconds.
复制代码
可见上述方法在2009a以外的版本中不失为一种方法，起码效率高于dblquad帮助文档里推荐的
做法。
更重要的是，这给我们求解一般区域三重积分提供了一种途径。
由于时间太晚了，求一般区域三重积分的方法留待下来有时间再写。
回复引用
订阅报告道具TOP
∙发短消息
∙加为好友
spirit3772（蓝精灵）当前离线
Matlab
帖子
1658
积分
2678
技术分
92技术分
研学币
1772研学币
学历
在读博士生
研究方向
信息与通信工程：
信道均衡
职业
学生
所在地
北京海淀
在线时间
50小时
注册时间
2004-5-11
最后登录
2009-8-29
2#
发表于2009-6-1617:
21|只看该作者
赞rocwoods的原创！
对于积分补充一点，有人经常会问对于有参数的积分怎么做。
问题：
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.fori=1:
length（a）
4. result（i）=quadl（@（x,a）sin（a*x）,pi/4,pi/2,[],[],a（i））
5.end
复制代码
quadl类函数在输入其自身的5个参数外，还可以将被积函数的除被积变量外的其他函数参数按照顺序放在函数参数后面。
上述用法需要注意quadl函数自身的参数需要写全，比如示例中的两个[]需要写上，其代表的参数含义可参考help文档。
quadl（fun,a,b,tol,trace）withnon-zerotraceshowsthevaluesof[fcntab-aq]duringtherecursion.
上述用法quadl函数的help文档中没有给出，但在AnonymousFunctions的介绍中有说明，包括rocwoods所说的多重匿名函数
PS：
另外这个问题当然也可以采用其他常规的方法，如应用符号工具箱及多重循环
[本帖最后由spirit3772于2009-6-1618:
52编辑]
1
评分次数
∙ 技术分+3技术分
rocwoods
仁者乐山，智者乐水，儒者乐研，愚者乐学，知足常乐，乐在山水，乐在研学，其乐无穷！
欢迎到工程软件区进行交流！
回复引用
评分报告道具TOP
rocwoods
∙发短消息
∙加为好友
rocwoods当前离线
帖子
1183
积分
2865
技术分
150技术分
研学币
1286研学币
学历
研究方向
职业
所在地
在线时间
91小时
注册时间
2005-9-13
最后登录
2009-8-30
版主
3#
发表于2009-6-1617:
58|只看该作者
spirit3772版主给出的方法不错
更多带参数积分的解决方法可以看这个链接
上面列出的各种方法中有利用到inline结构的，不过使用MATLAB7以后版本的建议不要再使用inline了，效率比匿名函数低很多。
回复引用
评分报告道具TOP
rocwoods
∙发短消息
∙加为好友
rocwoods当前离线
帖子
1183
积分
2865
技术分
150技术分
研学币
1286研学币
学历
研究方向
职业
所在地
在线时间
91小时
注册时间
2005-9-13
最后登录
2009-8-30
版主
4#
发表于2009-6-1709:
53|只看该作者
前面讨论的一般二重积分的例子：
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）,toc
2.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
3.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
4.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
5.ans=
6. 179.2969
7.Elapsedtimeis0.037453seconds.
8.ans=
9. 179.2969
10.Elapsedtimeis0.223533seconds.
11.ans=
12. 179.2969
13.Elapsedtimeis0.090477seconds.
14.ans=
15. 178.9301
16.Elapsedtimeis78.421721seconds.
复制代码
可见如果大家用的是2009a，那么计算一般区域三重积分时可以用模板1，2009a以前的版本（当然都是7.X版本）可以用模板3，而模板2效率比较低（似乎是更加频繁调用函数增加
了系统开销）。
以上二重三重积分模板还可以推广应用范围，譬如计算int（1/int（y,x,x^2）,10,100）就不能套用dblquad或者quad2d，但是用一般二重积分模板稍作变形，可以这样求解：
1.quadl（@（x）1./arrayfun（@（xx）quadl（@（y）y,xx,xx^2）,x）,10,100）
升级会员 