定積分的近似計算

1、數(shù)學(xué)軟件與數(shù)學(xué)實驗專題 定積分的近似計算數(shù)學(xué)軟件與數(shù)學(xué)實驗數(shù)學(xué)軟件與數(shù)學(xué)實驗2q 問題背景和實驗?zāi)康膶ｎ}專題 定積分的近似計算定積分的近似計算u 定積分計算的基本公式是牛頓萊布尼茲公式。但當(dāng)被積函數(shù)的原函數(shù)不知道時，如何計算？這時就需要利用近似計算。特別是在許多實際應(yīng)用中，被積函數(shù)甚至沒有解析表達(dá)式，而是一條實驗記錄曲線，或一組離散的采樣值，此時只能用近似方法計算定積分。u 本實驗主要研究定積分的幾種近似計算算法：矩形法、梯形法和拋物線法；蒙特卡洛隨機投點法和蒙特卡洛樣本均值法。同時介紹 Matlab 計算定積分的相關(guān)函數(shù)。數(shù)學(xué)軟件與數(shù)學(xué)實驗數(shù)學(xué)軟件與數(shù)學(xué)實驗3q 矩形法u 定積分的定義：專題

2、專題 定積分的近似計算定積分的近似計算badxxf)(,1iiixx0 x1x2x1nxnx 1x2x1ix ix nx,)(iixf ,1iiixxxiixxmax0limxnixni 1數(shù)學(xué)軟件與數(shù)學(xué)實驗數(shù)學(xué)軟件與數(shù)學(xué)實驗4矩形矩形法法1( )(), nbiiaif x dxfx n 充分大，x 充分小u 定積分的近似：l 通常我們?nèi)xxx21nabh左點法右點法中點法l 點 可以任意選取，常見的取法有： 左端點 ，右端點 和中點 。,1iiixx1ixix2/ )(1iixx數(shù)學(xué)軟件與數(shù)學(xué)實驗數(shù)學(xué)軟件與數(shù)學(xué)實驗5()/ixhban, 1 2ixaihi, ,n步長節(jié)點u 右點法：11(

3、 )()()nnbiiaiif x dxfxhf xixu 中點法：1111( )()22(nnbiiiaiixxf x dxfxhfiixx111( )()()nnbiiaiif x dxfxhf x-1ixu 左點法：左點左點法、右點法和中點法法、右點法和中點法數(shù)學(xué)軟件與數(shù)學(xué)實驗數(shù)學(xué)軟件與數(shù)學(xué)實驗6解：矩形法舉例矩形法舉例h =1/100=0.01, xi = i*h, a=0, b=1, n=100 u 例：用不同的矩形法計算下面的定積分 ( 取 n=100 )， 并比較這三種方法的相對誤差。1021xdxl 左點法：niiniixhxfhxdx1121110211 )(1 0.7878

4、9399673078l 右點法：1201()1niidxhf xx 0.78289399673078 0.78540024673078niiixxfhxdx11102)2(1l 中點法：(i = 0,1,2,.,100)數(shù)學(xué)軟件與數(shù)學(xué)實驗數(shù)學(xué)軟件與數(shù)學(xué)實驗711020arctan1dxxx40.78789399673078/4/40.003178l 理論值：l 左點法相對誤差：0.78289399673078/4/40.0031880.78540024673078/4/4 -62.65310 u 誤差分析矩形法舉例矩形法舉例l 右點法相對誤差：l 中點法相對誤差：不同的方法有不同的計算精度有沒

5、有更好的近似計算定積分的方法有沒有更好的近似計算定積分的方法 ?數(shù)學(xué)軟件與數(shù)學(xué)實驗數(shù)學(xué)軟件與數(shù)學(xué)實驗8)(xfab1ixixxyobadxxfS)(1S2SiSnSniibaSdxxfS1)(定積分幾何意義定積分幾何意義數(shù)學(xué)軟件與數(shù)學(xué)實驗數(shù)學(xué)軟件與數(shù)學(xué)實驗9iSiiiixyyS21nixfyii, 2 , 1 ),(u 曲邊小梯形的面積可以由直邊小梯形的面積來近似u 整個曲邊梯形的面積：badxxfS)(iniiiniixyyS1112梯形法梯形法數(shù)學(xué)軟件與數(shù)學(xué)實驗數(shù)學(xué)軟件與數(shù)學(xué)實驗10u 如果我們 n 等分區(qū)間 a,b，即令：badxxfS)(niiiiniiiniiyyhxyyS11111

6、22nxxx21nabh則則22)(110nnbayyyyhdxxf梯形公式梯形公式梯形法梯形法梯形公式與中點公式有什么區(qū)別梯形公式與中點公式有什么區(qū)別 ?數(shù)學(xué)軟件與數(shù)學(xué)實驗數(shù)學(xué)軟件與數(shù)學(xué)實驗11解：101120122nndxyyhyyx 0.78539399673078u 例：用梯形法計算下面定積分 ( 取 n=100 )， 并計算相對誤差梯形法舉例梯形法舉例1021xdxa=0, b=1, n=100, f (x) = 1/( 1+x2 ) h =1/100=0.01, xi = i*h, yi = f (xi) l 相對誤差：0.78539399673078/4/4-65.30510 數(shù)

7、學(xué)軟件與數(shù)學(xué)實驗數(shù)學(xué)軟件與數(shù)學(xué)實驗12u 2n 等分區(qū)間 a,b ，得該直線用拋物線代替，計算精度是否會更好？11, , 0,1,22ibahxihinnu 計算每個節(jié)點上的函數(shù)值：nixfyii2 , 1 , 0 ),(拋物線法拋物線法u 在區(qū)間 x0, x2 上，用過以下三點),( ),( ),(222111000yxPyxPyxP的拋物線來近似原函數(shù) f (x) 。數(shù)學(xué)軟件與數(shù)學(xué)實驗數(shù)學(xué)軟件與數(shù)學(xué)實驗13u 設(shè)過以上三點的拋物線方程為：則在則在區(qū)間 x0, x2 上，有上，有y = x2 + x + = p1(x) 2020)()(1xxxxdxxpdxxf20) (2xxdxxx20

8、2 3 23xxxxx20012 (4)6xxyyy012 (4)6bayyyn拋物線法拋物線法數(shù)學(xué)軟件與數(shù)學(xué)實驗數(shù)學(xué)軟件與數(shù)學(xué)實驗14u 同理可得：)4(6)( )4(6)(2122243222242nnnxxxxyyynabdxxfyyynabdxxfnnu 相加即得：2221222121( )( ) (4) 6iinbxaxiniiiif x dxf x dxbayyyn拋物線法拋物線法數(shù)學(xué)軟件與數(shù)學(xué)實驗數(shù)學(xué)軟件與數(shù)學(xué)實驗15u 整理后可得：)(2 )(46 )(2242123120nnnbayyyyyyyynabdxxf或辛普森 (Simpson) 公式拋物線法公式拋物線法拋物線法數(shù)學(xué)

9、軟件與數(shù)學(xué)實驗數(shù)學(xué)軟件與數(shù)學(xué)實驗16)(461123120102nnyyyyynabxdx 0.78539816339745)(22242nyyyu 例：用拋物線法計算下面定積分 ( 取 n=100 )， 并計算相對誤差1021xdx解：a=0, b=1, n=100, yi = f (xi) = 1/( 1+xi2 ) 0.78539816339745/4/4-162.82710 l 相對誤差：拋物線法拋物線法數(shù)學(xué)軟件與數(shù)學(xué)實驗數(shù)學(xué)軟件與數(shù)學(xué)實驗17u 梯形法：trapztrapz(x,y)x 為分割點（節(jié)點）組成的向量，y 為被積函數(shù)在節(jié)點上的函數(shù)值組成的向量。22)(110nnbayyy

10、ynabdxxf,x10nxxx01y (),(),()nf xf xf xq Matlab 近似計算定積分的相關(guān)函數(shù)Matlab 計算定積分函數(shù)介紹計算定積分函數(shù)介紹數(shù)學(xué)軟件與數(shù)學(xué)實驗數(shù)學(xué)軟件與數(shù)學(xué)實驗18前面的做法u 例：用梯形法計算下面定積分 ( 取 n=100) 1021xdx解：a=0, b=1, n=100, yi = f (xi) = 1/( 1+xi2 ) x=0:1/100:1; y=1./(1+x.2); trapz(x, y)trapz函數(shù)1012120(/2/2)1nndxbayyyyyxntrapz(x,1./(1+x.2)trapz 舉例舉例數(shù)學(xué)軟件與數(shù)學(xué)實驗數(shù)學(xué)軟

11、件與數(shù)學(xué)實驗19quad(f,a,b,tol)f = f(x) 為被積函數(shù)，a,b 為積分區(qū)間，tol 為計算精度將自變量看成是將自變量看成是向量向量badxxf)(u 拋物線法：quadl 不用自己分割積分區(qū)間l 可以指定計算精度，若不指定，缺省精度是 10-6l 精度越高，函數(shù)運行的時間越長l 此處的函數(shù) f 是數(shù)值形式，應(yīng)該使用數(shù)組運算，即 點運算：.*，./ ，. ，. 注：拋物線法拋物線法數(shù)學(xué)軟件與數(shù)學(xué)實驗數(shù)學(xué)軟件與數(shù)學(xué)實驗201021xdx解：解： quad(1./(1+x.2),0,1) quad(1./(1+x.2),0,1,10e-10) quad(1./(1+x.2),0,

12、1,10e-16)函數(shù)表達(dá)式一定要用函數(shù)表達(dá)式一定要用 單引號單引號 括起來！括起來！涉及的運算一定要用涉及的運算一定要用 數(shù)組運算數(shù)組運算！u 例：用 quad 計算定積分：quad 舉例舉例數(shù)學(xué)軟件與數(shù)學(xué)實驗數(shù)學(xué)軟件與數(shù)學(xué)實驗21q 拋物線法計算二重積分： dblquaddblquad(f,a,b,c,d,tol)u tol 為計算精度，若不指定，則缺省精度為 10-6 badcdxdyyxf),(u f(x,y) 可以由 inline 定義，或通過一個函數(shù)句柄傳遞u a,b 是第一積分變量的積分區(qū)間，c,d 是第二積分變量 的積分區(qū)間按字母順序，大寫字母排在小寫字母的前面二重積分的計算二

13、重積分的計算數(shù)學(xué)軟件與數(shù)學(xué)實驗數(shù)學(xué)軟件與數(shù)學(xué)實驗2221201(43)Ixyydxdy f=inline(4*x*y+3*y2); I=dblquad(f, -1,1,0,2)u f(x,y) 中關(guān)于第一自變量的運算是數(shù)組運算， 即把 x 看成是向量，y 看成是標(biāo)量。 也可以全部采用數(shù)組運算例2：計算二重積分 20112)34(dxdyxxy dblquad(inline(4*x*y+3*x2),-1,1,0,2) dblquad(inline(4*x*y+3*x.2),-1,1,0,2)X例1：計算二重積分dblquad 舉例舉例數(shù)學(xué)軟件與數(shù)學(xué)實驗數(shù)學(xué)軟件與數(shù)學(xué)實驗23例：計算二重積分 20

14、112)34(dxdyxxy dblquad(x,y)4*x*y+3*x.2 , -1,1, 0, 2)指定 x、y 分別是第一第一和第二第二積分變量 dblquad(inline(4*x*y+3*x.2), -1,1, 0, 2)q 被積函數(shù) f (x,y) 的另一種定義方法：匿名函數(shù)dblquad 舉例舉例數(shù)學(xué)軟件與數(shù)學(xué)實驗數(shù)學(xué)軟件與數(shù)學(xué)實驗24int(f,a,b) 計算 f 關(guān)于默認(rèn)自變量 的定積分，積分區(qū)間為a,b。int(f) 計算 f 關(guān)于默認(rèn)自變量 的不定積分。int(f,v,a,b) 計算函數(shù) f 關(guān)于自變量 v 的定積分，積分區(qū)間為 a, bint(f,v) 計算函數(shù) f 關(guān)

15、于自變量 v 的不定積分badvvf)( )f v dvfindsym(f,1)q 符號積分： intint 符號積分符號積分?jǐn)?shù)學(xué)軟件與數(shù)學(xué)實驗數(shù)學(xué)軟件與數(shù)學(xué)實驗25 syms x y; f=y*sin(x); int(f,x) int(f,y) int(f) int(a+b)ans=-y*cos(x)ans=1/2*y2*sin(x)ans=-y*cos(x)ans=a*b+1/2*b2u 例：指出下面各條語句的輸出結(jié)果int 舉例舉例數(shù)學(xué)軟件與數(shù)學(xué)實驗數(shù)學(xué)軟件與數(shù)學(xué)實驗26u 例：用 int 函數(shù)計算定積分：解： syms x; f=1/(1+x2); int(f,x,0,1) f=sym

16、(1/(1+x2); int(f,x,0,1) int(1/(1+x2),x,0,1)或 int(1/(1+x2),0,1)1021xdx或或int 舉例舉例數(shù)學(xué)軟件與數(shù)學(xué)實驗數(shù)學(xué)軟件與數(shù)學(xué)實驗27double(a) 將將 a 轉(zhuǎn)化為雙精度型，若轉(zhuǎn)化為雙精度型，若 a 是字符，則取對應(yīng)的是字符，則取對應(yīng)的 ASCII 碼碼 a=3; double(a) double(a)例：ans = 3ans = 97其它相關(guān)函數(shù)其它相關(guān)函數(shù)數(shù)學(xué)軟件與數(shù)學(xué)實驗數(shù)學(xué)軟件與數(shù)學(xué)實驗28221dxex x=1:0.001:2; y=exp(x.(-2); trapz(x,y)l 梯形法：l 拋物線法： quad(

17、exp(x.(-2),1,2,10e-10)l 符號積分法： syms x int(exp(x(-2),x,1,2)例 1：用 Matlab 函數(shù)近似計算積分?jǐn)?shù)值實驗數(shù)值實驗數(shù)學(xué)軟件與數(shù)學(xué)實驗數(shù)學(xué)軟件與數(shù)學(xué)實驗2920112)(dyyxdxl 拋物線法： dblquad(inline(x+y2),0,2,-1,1)l 符號積分法： f=int(x+y2,y,-1,1); int(f,0,2)數(shù)值實驗數(shù)值實驗例 2：用 Matlab 函數(shù)近似計算二重積分?jǐn)?shù)學(xué)軟件與數(shù)學(xué)實驗數(shù)學(xué)軟件與數(shù)學(xué)實驗30q 用蒙特卡羅用蒙特卡羅 ( Monte Carlo ) 法近似計算法近似計算n=100000; %總模

18、擬次數(shù)總模擬次數(shù)m=0; %記錄有利實驗次數(shù)的變量初始化記錄有利實驗次數(shù)的變量初始化for i=1:n %開始模擬實驗開始模擬實驗 x=rand(1); %產(chǎn)生隨機點產(chǎn)生隨機點 y=rand(1); %產(chǎn)生隨機點產(chǎn)生隨機點 if ( 1/(x2 + 1) = y ) %判斷試驗點是否落在積分線下方判斷試驗點是否落在積分線下方 m=m+1; %有利實驗次數(shù)加有利實驗次數(shù)加1 endendfprintf(面積的近似值為：面積的近似值為：%fn, m/n)蒙特卡羅法近似計算定積分蒙特卡羅法近似計算定積分-投點法投點法1021xdx數(shù)學(xué)軟件與數(shù)學(xué)實驗數(shù)學(xué)軟件與數(shù)學(xué)實驗31蒙特卡羅法近似計算定積分蒙特卡

19、羅法近似計算定積分-樣本平均值樣本平均值法法設(shè)設(shè)g(x)是是( (a, b) )上的一個密度函數(shù)，改寫上的一個密度函數(shù)，改寫 dxxfba XgXfEdxxgxgxfbaniiixgxfnXgXfE1)()(1)()(n基本原理基本原理：對積分：對積分 其中，其中，X是服從是服從g(x)的隨機變量可見，積分可以表示的隨機變量可見，積分可以表示為為X的函數(shù)的期望。由矩法，若有的函數(shù)的期望。由矩法，若有n個來自個來自g( (x) )的觀測值的觀測值x1 1, , ,xn，則可給出，則可給出 的一個矩估計。的一個矩估計。數(shù)學(xué)軟件與數(shù)學(xué)實驗數(shù)學(xué)軟件與數(shù)學(xué)實驗32特別地，若特別地，若a，b有限，可取有限，可取 g(x) 為為 a,ba,b 上均勻分布此時上均勻分布此時, ,設(shè)設(shè)x1, xn是來自是來自U(a,b)的隨機數(shù)，則的隨機數(shù)，則 的一個估計為的一個估計為 niiniiixfnabxgxfn111 用上式來估計，和計算，隨機數(shù)