第10章 数值积分计算

2023/2/2数值积分计算吴鹏（rocwoods）rocwoods@126.comMATLAB从零到进阶2023/2/2主要内容矩形区域积分以及离散数据积分

含参数的积分一般区域二重、三重积分一般区域n重积分

蒙特卡洛法计算n重积分

2023/2/2第一节矩形区域积分以及离散数据积分2023/2/2一、矩形区域积分积分上下限是常数，分单重、二重、三重三种情况1.单重quad(自适应Simpson积分)quadl(自适应Gauss-Lobatto积分，最常用)quadgk(自适应Gauss-Kronrod积分，尤其适合震荡积分、含奇点的积分；从R2007b开始支持)

具体用法参见书中【例10.1-1】2.二重情形一般二重积分可以用下面两个函数解决：dblquadquad2d(R2009a开始支持，不仅可以求矩形区域二重积分，还可以求一般区域二重积分)具体用法参见书中【例10.1-2】2023/2/2一、矩形区域积分3.三重情形MATLAB中可以利用triplequad函数求解长方体区域的三重积分。具体用法参见书中【例10.1-3】4.向量化积分所谓向量化积分，是指被积函数含有参数，需要对参数的一系列值求出相应的积分。向量话积分可以用quadv函数。具体用法参见书中【例10.1-4】2023/2/2二、离散数据积分离散数据积分的函数MATLAB只有trapz函数，是针对一重情形的，如果是二重情形，需要自己编写，其中会用到trapz函数。用法示例参见书中【例10.1-5】以及【例10.1-6】2023/2/2第二节含参数的积分2023/2/2

在解决实际问题的时候我们常要求解一些带参数的积分问题，这涉及到参数的传递，总的说来有以下几种方法可以求解：用inline+num2str方法(MATLAB7.0以前版本使用，7.0以后不推荐使用)利用匿名函数实现利用嵌套函数实现利用积分函数本身传递参数 这几种方法的具体实现办法可以参见【例10.2-1】。2023/2/2第三节一般区域二重、三重积分2023/2/2一、概要

一般区域上的二重积分，低版本的MATLAB向来支持不好，7.X版本之前不能通过简单的形式直接求取一般区域上的的二重积分，往往要借助广泛流传的NIT工具箱（数值积分工具箱）来实现。后来到了7.X版本，MATLAB引入了匿名函数结构，dblquad的被积函数可以是匿名函数的形式，利用匿名函数结构，通过适当改写被积函数。dblquad可以求解一般区域上的二重积分。但这种改写就是将积分区域表示成逻辑表达式形式并与被积函数相乘，本质上还是计算的矩形区域的积分。 仔细分析dblquad的帮助文档中求一般区域二重积分的例子就会发现，这种办法是把原被积函数外推到了一个矩形区域上来实现的，该矩形区域包含了被积区域，在被积区域上，外推函数取值和原函数一样，而在矩形区域内被积区域外的那部分，外推函数取值都为0.这样外推函数在矩形区域上积分的结果就等于原被积函数在被积区域的结果。这种数学上简单的变换导致的结果就是，按这种办法，计算量大增，特别是被积区域非常不规则并布满整个矩形区域的时候。2023/2/2一、概要

从MATLABR2009a版本起，MATLAB终于有了专门求解一般区域二重积分的函数——quad2d，该函数采用的自适应积分算法基于L.F.Shampine的文章："VectorizedAdaptiveQuadratureinMATLAB,"JournalofComputationalandAppliedMathematics,211,2008。该函数求解一般区域二重积分的效率要远高于上面提到的dblquad的方法。 虽然quad2d可以求解一般区域二重积分，但还是不能直接求解一般区域三重积分，而NIT工具箱也没有一般区域三重积分的计算函数。这些需要灵活运用MATLAB相关函数来构建。2023/2/2二、一般区域二重积分的计算dblquad求解一般区域二重三重积分的思路方法，就是将被积函数“延拓”到矩形或者长方体区域，但是这种方法不可避免引入很多乘0运算浪费时间。而quad2d虽然效率很高，但是R2009a之前没有这个函数。因此，要想在2009a之前的MATLAB7版里里计算一般区域二重三重积分，需要开辟新的方法。新的方法是调用已有的MATLAB函数求解，在求一般区域二重积分时，效率和quad2d相比有一些差距，但是相对于“延拓”函数的做法，效率大大提高了

。【例10.3-2】和【例10.3-3】给出了这一方法的应用实例，可以看出该方法在R2009a以前的版本中不失为一种方法，效率要明显高于dblquad帮助文档里推荐的做法。更重要的是，这给我们求解一般区域三重积分提供了一种途径。2023/2/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积分上限值)模板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))模板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积分上限值)模板使用说明：x,y,z是积分变量，模板中除了用语言描述的参量用相应表达式替换掉外，其余结构保持不变。具体使用示例可以参考【例10.3-4】

2023/2/2第四节一般区域n重积分的计算2023/2/2介绍求解一般区域n重积分的函数——nIntegratef=nIntegrate(fun,low,up)

f为函数的返回值是n重积分积分结果。fun是被积函数字符串形式，不同的变量依次以x1,x2,...xn表示，（需要注意的是，必须以x1,x2,...xn这种形式表示，其余像y1,y2,...yn或是其他表示方法都不行）。low和up都是长度为n的cell数组，low存储从外到内各重积分的积分下限函数，up存储从外到内各重积分的积分上限函数（都是字符串形式）。low和up内的函数表示都要遵循一些原则，这些原则在程序注释里进行了说明.nIntegrate函数源代码以及使用示例参考书中相关内容以及【例10.4-1】至【例10.4-4】2023/2/2第五节蒙特卡洛法计算n重积分2023/2/2 一些积分重数比较高的积分，按照传统方法进行计算的话，运算量会非常大，求解时间随着重数的增加而迅速增加。在10.4节虽然给出了一般区域n重积分的通用程序，但是在实际计算中，n不会很大，一般n超过5后，求解时间就会长到无法接受。对于更高维的积分，一个非常有效的求解方法就是蒙特卡洛积分法。 蒙特卡洛方法又称随机抽样法或统计试验方法，它用于求积分时，与积分重数无关，这点非常重要。虽然四维以下的积分用蒙特卡洛法效率可能不如传统的一些数值积分方法，但是维数高的时候，蒙特卡洛法比传统方法要有效的多，而且实现起来也非常容易。可以说，计算高维积分是蒙特卡洛方法最成功和典型的应用。 实际应用中，有多种蒙特卡洛方法可以计算n重积分，比较常用的是基本蒙特卡洛法和等分布序列的蒙特卡洛法。2023/2/2