下面是小编为大家整理的一般区域二重三重积分MATLAB计算方法（完整文档）,供大家参考。
[ [ 【原创】] ]
一般区域二重、三重积分 B MATLAB 计算方法
这里讨论的计算方法指的是利用现有的 MATLAB 函数来求解，而不是根据具体的数值计算方法来编写相应程序。目前最新版的 2009a 有关于一般区域二重积分的计算函数 quad2d（详细
介绍见），但没有一般区域三重积分的计算函数，而 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.
复制代码
如果你用的是 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帮助文档里推荐的做法。更重要的是，这给我们求解一般区域三重积分提供了一种途径。
由于时间太晚了，求一般区域三重积分的方法留待下来有时间再写。
订阅
发表于赞对于积分补充一点，有人经常会问对于有参数的积分怎么做。问题：y=sin(a*x)下如何计算这里给出个1.3. 4.sin(a*x),pi/4,pi/2,[],[],a
quadl 类函数在输入其自身的 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：另外这个问题当然也可以采用其他常规的方法，如应用符号工具箱及多重循环
[
本帖最后由 spirit3772 于 2009-6-16 18:52 编辑 ] 仁者乐山，智者乐水，儒者乐研，愚者乐学，知足常乐，乐在山水，乐在研学，其乐无穷！欢迎到工程软件区进行交流！
评分
rocwoods
版主
3#
发表于 2009-6-16 17:58 | 只看该作者
spirit3772 版主给出的方法不错 更多带参数积分的解决方法可以看这个链接
上面列出的各种方法中有利用到inline结构的，不过使用MATLAB7以后版本的建议不要再使用 inline 了，效率比匿名函数低很多。
回复 引用
评分 报告 道具 TOP
rocwoods
版主
4#
发表于 2009-6-17 09: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. Elapsed time is 0.037453 seconds. 8. ans = 9.
179.2969 10. Elapsed time is 0.223533 seconds. 11. ans = 12.
179.2969 13. Elapsed time is 0.090477 seconds. 14. ans = 15.
178.9301 16. 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. quadl(@(x) 1./arrayfun(@(xx)quadl(@(y)y,xx,xx^2),x),10,100)