數值積分Matlab實現研究

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

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

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

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

5、01. 02)3 . 0/(. 1(2; 2xxinlineyysyms( (二二) ) 、二元函數數值積分指令、二元函數數值積分指令1.dblquad1.dblquad指令指令dblquaddblquad指令是指令是MatlabMatlab系統提供的可求解矩形區(qū)域上二重系統提供的可求解矩形區(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計算積分可在命令窗口鍵入如下指令可在命令窗口鍵入如下指令: :)2 , 1 , 1 , 0),log(*).2.2.)(,(yxyxyxdblquad2、quad2d指令格式:quad2d(fun,ymin,ymax,xflowfun,xhighfun)說明:1.quad2d是matlab系統從R2009a版本起提供在一般區(qū)域上的二重積分指令;2.xlowfun,xhighfun分別表示函數形式的積分上下限;3.指令所采用的算法是將一般的積分區(qū)域映射到矩形區(qū)域,然后

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

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

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

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

11、+y2=16,上述積分轉化為累次積分的積分限分別為:于是,累次積分的表示式為:44161616022222)1ln()(xxyxdzzxyxzdydx如何計算該三重積分?分析:利用matlab系統提供的一元函數及二元函數積分指令和arrayfun進行組合求解.方法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)結果:1.9302e+003, 運行時間: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:結果:1.9302e+003,fcp 運行時間:0.468731秒.方