第10章 數(shù)值積分計算

2023/2/2數(shù)值積分計算吳鵬（rocwoods）rocwoods@126.comMATLAB從零到進階2023/2/2主要內(nèi)容矩形區(qū)域積分以及離散數(shù)據(jù)積分

含參數(shù)的積分一般區(qū)域二重、三重積分一般區(qū)域n重積分

蒙特卡洛法計算n重積分

2023/2/2第一節(jié)矩形區(qū)域積分以及離散數(shù)據(jù)積分2023/2/2一、矩形區(qū)域積分積分上下限是常數(shù)，分單重、二重、三重三種情況1.單重quad(自適應(yīng)Simpson積分)quadl(自適應(yīng)Gauss-Lobatto積分，最常用)quadgk(自適應(yīng)Gauss-Kronrod積分，尤其適合震蕩積分、含奇點的積分；從R2007b開始支持)

具體用法參見書中【例10.1-1】2.二重情形一般二重積分可以用下面兩個函數(shù)解決：dblquadquad2d(R2009a開始支持，不僅可以求矩形區(qū)域二重積分，還可以求一般區(qū)域二重積分)具體用法參見書中【例10.1-2】2023/2/2一、矩形區(qū)域積分3.三重情形MATLAB中可以利用triplequad函數(shù)求解長方體區(qū)域的三重積分。具體用法參見書中【例10.1-3】4.向量化積分所謂向量化積分，是指被積函數(shù)含有參數(shù)，需要對參數(shù)的一系列值求出相應(yīng)的積分。向量話積分可以用quadv函數(shù)。具體用法參見書中【例10.1-4】2023/2/2二、離散數(shù)據(jù)積分離散數(shù)據(jù)積分的函數(shù)MATLAB只有trapz函數(shù)，是針對一重情形的，如果是二重情形，需要自己編寫，其中會用到trapz函數(shù)。用法示例參見書中【例10.1-5】以及【例10.1-6】2023/2/2第二節(jié)含參數(shù)的積分2023/2/2

在解決實際問題的時候我們常要求解一些帶參數(shù)的積分問題，這涉及到參數(shù)的傳遞，總的說來有以下幾種方法可以求解：用inline+num2str方法(MATLAB7.0以前版本使用，7.0以后不推薦使用)利用匿名函數(shù)實現(xiàn)利用嵌套函數(shù)實現(xiàn)利用積分函數(shù)本身傳遞參數(shù) 這幾種方法的具體實現(xiàn)辦法可以參見【例10.2-1】。2023/2/2第三節(jié)一般區(qū)域二重、三重積分2023/2/2一、概要

一般區(qū)域上的二重積分，低版本的MATLAB向來支持不好，7.X版本之前不能通過簡單的形式直接求取一般區(qū)域上的的二重積分，往往要借助廣泛流傳的NIT工具箱（數(shù)值積分工具箱）來實現(xiàn)。后來到了7.X版本，MATLAB引入了匿名函數(shù)結(jié)構(gòu)，dblquad的被積函數(shù)可以是匿名函數(shù)的形式，利用匿名函數(shù)結(jié)構(gòu)，通過適當(dāng)改寫被積函數(shù)。dblquad可以求解一般區(qū)域上的二重積分。但這種改寫就是將積分區(qū)域表示成邏輯表達式形式并與被積函數(shù)相乘，本質(zhì)上還是計算的矩形區(qū)域的積分。 仔細分析dblquad的幫助文檔中求一般區(qū)域二重積分的例子就會發(fā)現(xiàn)，這種辦法是把原被積函數(shù)外推到了一個矩形區(qū)域上來實現(xiàn)的，該矩形區(qū)域包含了被積區(qū)域，在被積區(qū)域上，外推函數(shù)取值和原函數(shù)一樣，而在矩形區(qū)域內(nèi)被積區(qū)域外的那部分，外推函數(shù)取值都為0.這樣外推函數(shù)在矩形區(qū)域上積分的結(jié)果就等于原被積函數(shù)在被積區(qū)域的結(jié)果。這種數(shù)學(xué)上簡單的變換導(dǎo)致的結(jié)果就是，按這種辦法，計算量大增，特別是被積區(qū)域非常不規(guī)則并布滿整個矩形區(qū)域的時候。2023/2/2一、概要

從MATLABR2009a版本起，MATLAB終于有了專門求解一般區(qū)域二重積分的函數(shù)——quad2d，該函數(shù)采用的自適應(yīng)積分算法基于L.F.Shampine的文章："VectorizedAdaptiveQuadratureinMATLAB,"JournalofComputationalandAppliedMathematics,211,2008。該函數(shù)求解一般區(qū)域二重積分的效率要遠高于上面提到的dblquad的方法。 雖然quad2d可以求解一般區(qū)域二重積分，但還是不能直接求解一般區(qū)域三重積分，而NIT工具箱也沒有一般區(qū)域三重積分的計算函數(shù)。這些需要靈活運用MATLAB相關(guān)函數(shù)來構(gòu)建。2023/2/2二、一般區(qū)域二重積分的計算dblquad求解一般區(qū)域二重三重積分的思路方法，就是將被積函數(shù)“延拓”到矩形或者長方體區(qū)域，但是這種方法不可避免引入很多乘0運算浪費時間。而quad2d雖然效率很高，但是R2009a之前沒有這個函數(shù)。因此，要想在2009a之前的MATLAB7版里里計算一般區(qū)域二重三重積分，需要開辟新的方法。新的方法是調(diào)用已有的MATLAB函數(shù)求解，在求一般區(qū)域二重積分時，效率和quad2d相比有一些差距，但是相對于“延拓”函數(shù)的做法，效率大大提高了

?！纠?0.3-2】和【例10.3-3】給出了這一方法的應(yīng)用實例，可以看出該方法在R2009a以前的版本中不失為一種方法，效率要明顯高于dblquad幫助文檔里推薦的做法。更重要的是，這給我們求解一般區(qū)域三重積分提供了一種途徑。2023/2/2三、一般區(qū)域三重積分的計算模板1：quadl(@(x)arrayfun(@(xx)quad2d(被積函數(shù)f(xx,y,z)關(guān)于y,z變量的函數(shù)句柄,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(被積函數(shù)f(xx,yy,z)關(guān)于z變量的函數(shù)句柄,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(被積函數(shù)f(xx,yy,z)關(guān)于z變量的函數(shù)句柄,z積分下限z1(xx,yy),z積分上限z2(xx,yy)),y),y積分下限y1(xx),y積分上限y2(xx)),x),x積分下限值,x積分上限值)模板使用說明：x,y,z是積分變量，模板中除了用語言描述的參量用相應(yīng)表達式替換掉外，其余結(jié)構(gòu)保持不變。具體使用示例可以參考【例10.3-4】

2023/2/2第四節(jié)一般區(qū)域n重積分的計算2023/2/2介紹求解一般區(qū)域n重積分的函數(shù)——nIntegratef=nIntegrate(fun,low,up)

f為函數(shù)的返回值是n重積分積分結(jié)果。fun是被積函數(shù)字符串形式，不同的變量依次以x1,x2,...xn表示，（需要注意的是，必須以x1,x2,...xn這種形式表示，其余像y1,y2,...yn或是其他表示方法都不行）。low和up都是長度為n的cell數(shù)組，low存儲從外到內(nèi)各重積分的積分下限函數(shù)，up存儲從外到內(nèi)各重積分的積分上限函數(shù)（都是字符串形式）。low和up內(nèi)的函數(shù)表示都要遵循一些原則，這些原則在程序注釋里進行了說明.nIntegrate函數(shù)源代碼以及使用示例參考書中相關(guān)內(nèi)容以及【例10.4-1】至【例10.4-4】2023/2/2第五節(jié)蒙特卡洛法計算n重積分2023/2/2 一些積分重數(shù)比較高的積分，按照傳統(tǒng)方法進行計算的話，運算量會非常大，求解時間隨著重數(shù)的增加而迅速增加。在10.4節(jié)雖然給出了一般區(qū)域n重積分的通用程序，但是在實際計算中，n不會很大，一般n超過5后，求解時間就會長到無法接受。對于更高維的積分，一個非常有效的求解方法就是蒙特卡洛積分法。 蒙特卡洛方法又稱隨機抽樣法或統(tǒng)計試驗方法，它用于求積分時，與積分重數(shù)無關(guān)，這點非常重要。雖然四維以下的積分用蒙特卡洛法效率可能不如傳統(tǒng)的一些數(shù)值積分方法，但是維數(shù)高的時候，蒙特卡洛法比傳統(tǒng)方法要有效的多，而且實現(xiàn)起來也非常容易?？梢哉f，計算高維積分是蒙特卡洛方法最成功和典型的應(yīng)用。 實際應(yīng)用中，有多種蒙特卡洛方法可以計算n重積分，比較常用的是基本蒙特卡洛法和等分布序列的蒙特卡洛法。2023/2/2