數(shù)值積分Matlab實(shí)現(xiàn)研究

1、數(shù)值積分的數(shù)值積分的MatlabMatlab實(shí)現(xiàn)研究實(shí)現(xiàn)研究一、數(shù)值積分的一、數(shù)值積分的MatlabMatlab指令指令( (一一).).一元函數(shù)數(shù)值積分指令一元函數(shù)數(shù)值積分指令指令名稱指令名稱指令格式指令格式采用的算法采用的算法特點(diǎn)特點(diǎn)quadquad(fun,xmin, xmax,tol)自適應(yīng)自適應(yīng)simpsonsimpson算法算法閉型積分閉型積分, ,精精度一般度一般quadlquadl(fun,xmin, xmax,tol)自適應(yīng)自適應(yīng)LobattontLobattont算算法法閉型積分閉型積分, ,精精度一般度一般quadgkquadgk(fun,xmin, xmax)自適應(yīng)自適

2、應(yīng)Gass-Gass-KronrodKronrod算法算法開型積分開型積分, ,精精度較高度較高, ,支持支持無窮區(qū)間積分無窮區(qū)間積分, ,提供提供AbsTolAbsTol等選項(xiàng)等選項(xiàng)說明:1 1、funfun可以是匿名函數(shù)可以是匿名函數(shù), ,字符串形式函數(shù)字符串形式函數(shù), ,內(nèi)嵌函數(shù)內(nèi)嵌函數(shù),M,M文文件的函數(shù)句柄等件的函數(shù)句柄等2 2、函數(shù)表達(dá)式要遵循、函數(shù)表達(dá)式要遵循”數(shù)組運(yùn)算數(shù)組運(yùn)算”3 3、xmin,xmaxxmin,xmax分別為積分變量的上下限分別為積分變量的上下限4 4、toltol為要求的計(jì)算結(jié)果絕對誤差為要求的計(jì)算結(jié)果絕對誤差, ,缺省值為缺省值為1e-61e-6例題例題:

3、2122)(604. 0)9 . 0(101. 0) 3 . 0(1)(:dxxfxxxf計(jì)算已知計(jì)算步驟如下計(jì)算步驟如下: :; 6)04. 02).9 . 0/(. 1)01. 02).3 . 0/(. 1)(xxyxfyfunction1.1.在在MatlabMatlab編輯窗口中編寫如下編輯窗口中編寫如下M M文件文件, ,并存盤并存盤2:2:在在MatlabMatlab命令窗口中鍵入下列命令命令窗口中鍵入下列命令: :)2 , 1, (fquadq結(jié)果結(jié)果: :3450.26q參考文獻(xiàn):張德豐.Matlab數(shù)值分析與應(yīng)用M,國防工業(yè)出版社,2007.注:tocyquadwtic),2

4、 , 1, 1(,;6-0.04)0.9)2-1./(x0.01)0.3)2-1./(xy11、被積函數(shù)可以用字符形式表示被積函數(shù)可以用字符形式表示在命令窗口依次鍵入下列命令在命令窗口依次鍵入下列命令 其中其中:tic,toc:tic,toc是秒表計(jì)時(shí)命令是秒表計(jì)時(shí)命令,tic,tic表示秒表計(jì)時(shí)開始表示秒表計(jì)時(shí)開始,toc,toc表表示秒表計(jì)時(shí)結(jié)束示秒表計(jì)時(shí)結(jié)束, ,運(yùn)行花費(fèi)時(shí)間輸出格式為運(yùn)行花費(fèi)時(shí)間輸出格式為“elapsed_time=elapsed_time=”, , 單位為秒。單位為秒。2 2、被積函數(shù)用內(nèi)嵌函數(shù)表示、被積函數(shù)用內(nèi)嵌函數(shù)表示); 6)04. 02)9 . 0/(. 1)

5、01. 02)3 . 0/(. 1(2; 2xxinlineyysyms( (二二) ) 、二元函數(shù)數(shù)值積分指令、二元函數(shù)數(shù)值積分指令1.dblquad1.dblquad指令指令dblquaddblquad指令是指令是MatlabMatlab系統(tǒng)提供的可求解矩形區(qū)域上二重系統(tǒng)提供的可求解矩形區(qū)域上二重 ( (閉型閉型) )積分積分格式格式: : dblquad(fun,xmin,xmax,ymin,ymax,tol,method)dblquad(fun,xmin,xmax,ymin,ymax,tol,method)其中其中:method:method是積分方法的選擇是積分方法的選擇, ,缺省方

6、法是缺省方法是quad,quad,可以可以選擇選擇quadlquadl等等例題例題: :211022)ln()(xdy:dxyxy計(jì)算積分可在命令窗口鍵入如下指令可在命令窗口鍵入如下指令: :)2 , 1 , 1 , 0),log(*).2.2.)(,(yxyxyxdblquad2、quad2d指令格式:quad2d(fun,ymin,ymax,xflowfun,xhighfun)說明:1.quad2d是matlab系統(tǒng)從R2009a版本起提供在一般區(qū)域上的二重積分指令;2.xlowfun,xhighfun分別表示函數(shù)形式的積分上下限;3.指令所采用的算法是將一般的積分區(qū)域映射到矩形區(qū)域,然后

7、利用自適應(yīng)Lobatton算法進(jìn)行計(jì)算.例題:計(jì)算dxyxyxx)sin()(22211222 可用如下指令:)2.1 )(,)(, 2 , 1.),2.sin(*).2.2.)(,(2xxxxyxyxyxdquad(三) 、三元及以上函數(shù)數(shù)值積分指令1.triplequad指令triplequad是matlab系統(tǒng)提供的可求解長方體區(qū)域上三重積分格式:triplequad(fun,xmin,xmax,ymin,ymax,zmin,zmax,tol,method)算法:根據(jù)method的選擇,缺省方法 quad50211022)ln()(:dxzyxyxdydz如 可用如下指令:)5 , 0

8、, 2 , 1 , 1 , 0),2.2.(*).log(),(yxzyxzyxtriplequad2.quadndg指令quadndg指令是積分工具箱INT中的指令,用于計(jì)算n重超長方體區(qū)域上的n元函數(shù)的數(shù)值積分格式: quadndg(fun,xlowdim,xhighdim,tol)說用:1.xmindim,xmaxdim分別表示積分下限和上限的n維常量;2.采用Gauss積分方法,在瑕點(diǎn)處采用Gauss-Chebyshev方法502110sinz)dx(xydydz:如計(jì)算可用指令)81,5 , 2 , 1 ,0 , 1 , 0),3(sin()2(*) 1 ()(exxxxquadnd

9、g二、一般數(shù)值積分的討論與研究 由于數(shù)值積分非常復(fù)雜,如振蕩積分,無界區(qū)域的積分,無界函數(shù)的積分,高維積分等,目前matlab系統(tǒng)還沒有直接提供一般區(qū)域上的三元及三元以上函數(shù)積分的指令,這需要人們從問題出發(fā),利用matlab系統(tǒng)現(xiàn)有的指令,探索求解一般區(qū)域的積分問題.以具體例題進(jìn)行討論:deyyx)(x2D2計(jì)算積分例1其中D由直線y=x,拋物線y=1+x2以及軸所圍成的區(qū)域.分析:1.由matlab可以作出積分區(qū)域的圖形,可見上述問題為廣義二重積分問題,積分區(qū)域?yàn)闊o界區(qū)域,可轉(zhuǎn)化為累次積分.dyeyxdxyxxx)(201222.該累次積分如何計(jì)算? 由于matlab系統(tǒng)二重積分沒有指持無界

10、區(qū)域的積分,但一元函數(shù)的quadgk指令支持無界區(qū)域的積分,為此可以利用兩個(gè) quadgk指令和arrayfun指令進(jìn)行組合處理二元函數(shù)無界區(qū)域 上的數(shù)值積分3.要求廣義二重積分是收斂的.積分指令如下:quadgk(x)arrayfun(x)quadgk(y)exp(-x-y).*(x.2+y.2),x,1+x.2),x),0,inf)運(yùn)行結(jié)果:0.6427例2.計(jì)算三重積分dvzxyxz)1ln()(222222160;1616; 44yxzxyxx 其中積分區(qū)域由xoy坐標(biāo)面與旋轉(zhuǎn)拋物面z=16-x2-y2所圍成的立體區(qū)域. 由題意可知,旋轉(zhuǎn)拋物面在xoy平面的投影為一圓,其方程滿足:x2

11、+y2=16,上述積分轉(zhuǎn)化為累次積分的積分限分別為:于是,累次積分的表示式為:44161616022222)1ln()(xxyxdzzxyxzdydx如何計(jì)算該三重積分?分析:利用matlab系統(tǒng)提供的一元函數(shù)及二元函數(shù)積分指令和arrayfun進(jìn)行組合求解.方法1:quad2d(x,y)arrayfun(x,y)quadgk(z)(x.*z+y.2).*log(1+x.2+z),0,16-x.2-y.2),x,y),-4,4,(x)-sqrt(16-x.2),(x)sqrt(16-x.2)結(jié)果:1.9302e+003, 運(yùn)行時(shí)間:4.925542秒.quadgk(z)arrayfun(z)quad2d(x,y)(x.*z+y.2).*log(1+x.2+z),-sqrt(16-z),sqrt(16-z),(x)-sqrt(16-z-x.2),(x)sqrt(16-z-x.2),z),0,16)方法2:結(jié)果:1.9302e+003,fcp 運(yùn)行時(shí)間:0.468731秒.方