定積分的近似計(jì)算

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

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

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

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

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

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

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

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

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

10、ynabdxxf,x10nxxx01y (),(),()nf xf xf xq Matlab 近似計(jì)算定積分的相關(guān)函數(shù)Matlab 計(jì)算定積分函數(shù)介紹計(jì)算定積分函數(shù)介紹數(shù)學(xué)軟件與數(shù)學(xué)實(shí)驗(yàn)數(shù)學(xué)軟件與數(shù)學(xué)實(shí)驗(yàn)18前面的做法u 例：用梯形法計(jì)算下面定積分 ( 取 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í)驗(yàn)數(shù)學(xué)軟

11、件與數(shù)學(xué)實(shí)驗(yàn)19quad(f,a,b,tol)f = f(x) 為被積函數(shù)，a,b 為積分區(qū)間，tol 為計(jì)算精度將自變量看成是將自變量看成是向量向量badxxf)(u 拋物線法：quadl 不用自己分割積分區(qū)間l 可以指定計(jì)算精度，若不指定，缺省精度是 10-6l 精度越高，函數(shù)運(yùn)行的時(shí)間越長(zhǎng)l 此處的函數(shù) f 是數(shù)值形式，應(yīng)該使用數(shù)組運(yùn)算，即 點(diǎn)運(yùn)算：.*，./ ，. ，. 注：拋物線法拋物線法數(shù)學(xué)軟件與數(shù)學(xué)實(shí)驗(yàn)數(shù)學(xué)軟件與數(shù)學(xué)實(shí)驗(yàn)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á)式一定要用 單引號(hào)單引號(hào) 括起來(lái)！括起來(lái)！涉及的運(yùn)算一定要用涉及的運(yùn)算一定要用 數(shù)組運(yùn)算數(shù)組運(yùn)算！u 例：用 quad 計(jì)算定積分：quad 舉例舉例數(shù)學(xué)軟件與數(shù)學(xué)實(shí)驗(yàn)數(shù)學(xué)軟件與數(shù)學(xué)實(shí)驗(yàn)21q 拋物線法計(jì)算二重積分： dblquaddblquad(f,a,b,c,d,tol)u tol 為計(jì)算精度，若不指定，則缺省精度為 10-6 badcdxdyyxf),(u f(x,y) 可以由 inline 定義，或通過(guò)一個(gè)函數(shù)句柄傳遞u a,b 是第一積分變量的積分區(qū)間，c,d 是第二積分變量 的積分區(qū)間按字母順序，大寫字母排在小寫字母的前面二重積分的計(jì)算二

13、重積分的計(jì)算數(shù)學(xué)軟件與數(shù)學(xué)實(shí)驗(yàn)數(shù)學(xué)軟件與數(shù)學(xué)實(shí)驗(yàn)2221201(43)Ixyydxdy f=inline(4*x*y+3*y2); I=dblquad(f, -1,1,0,2)u f(x,y) 中關(guān)于第一自變量的運(yùn)算是數(shù)組運(yùn)算， 即把 x 看成是向量，y 看成是標(biāo)量。 也可以全部采用數(shù)組運(yùn)算例2：計(jì)算二重積分 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：計(jì)算二重積分dblquad 舉例舉例數(shù)學(xué)軟件與數(shù)學(xué)實(shí)驗(yàn)數(shù)學(xué)軟件與數(shù)學(xué)實(shí)驗(yàn)23例：計(jì)算二重積分 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í)驗(yàn)數(shù)學(xué)軟件與數(shù)學(xué)實(shí)驗(yàn)24int(f,a,b) 計(jì)算 f 關(guān)于默認(rèn)自變量 的定積分，積分區(qū)間為a,b。int(f) 計(jì)算 f 關(guān)于默認(rèn)自變量 的不定積分。int(f,v,a,b) 計(jì)算函數(shù) f 關(guān)于自變量 v 的定積分，積分區(qū)間為 a, bint(f,v) 計(jì)算函數(shù) f 關(guān)

15、于自變量 v 的不定積分badvvf)( )f v dvfindsym(f,1)q 符號(hào)積分： intint 符號(hào)積分符號(hào)積分?jǐn)?shù)學(xué)軟件與數(shù)學(xué)實(shí)驗(yàn)數(shù)學(xué)軟件與數(shù)學(xué)實(shí)驗(yàn)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 例：指出下面各條語(yǔ)句的輸出結(jié)果int 舉例舉例數(shù)學(xué)軟件與數(shù)學(xué)實(shí)驗(yàn)數(shù)學(xué)軟件與數(shù)學(xué)實(shí)驗(yàn)26u 例：用 int 函數(shù)計(jì)算定積分：解： 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í)驗(yàn)數(shù)學(xué)軟件與數(shù)學(xué)實(shí)驗(yàn)27double(a) 將將 a 轉(zhuǎn)化為雙精度型，若轉(zhuǎn)化為雙精度型，若 a 是字符，則取對(duì)應(yīng)的是字符，則取對(duì)應(yīng)的 ASCII 碼碼 a=3; double(a) double(a)例：ans = 3ans = 97其它相關(guān)函數(shù)其它相關(guān)函數(shù)數(shù)學(xué)軟件與數(shù)學(xué)實(shí)驗(yàn)數(shù)學(xué)軟件與數(shù)學(xué)實(shí)驗(yàn)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 符號(hào)積分法： syms x int(exp(x(-2),x,1,2)例 1：用 Matlab 函數(shù)近似計(jì)算積分?jǐn)?shù)值實(shí)驗(yàn)數(shù)值實(shí)驗(yàn)數(shù)學(xué)軟件與數(shù)學(xué)實(shí)驗(yàn)數(shù)學(xué)軟件與數(shù)學(xué)實(shí)驗(yàn)2920112)(dyyxdxl 拋物線法： dblquad(inline(x+y2),0,2,-1,1)l 符號(hào)積分法： f=int(x+y2,y,-1,1); int(f,0,2)數(shù)值實(shí)驗(yàn)數(shù)值實(shí)驗(yàn)例 2：用 Matlab 函數(shù)近似計(jì)算二重積分?jǐn)?shù)學(xué)軟件與數(shù)學(xué)實(shí)驗(yàn)數(shù)學(xué)軟件與數(shù)學(xué)實(shí)驗(yàn)30q 用蒙特卡羅用蒙特卡羅 ( Monte Carlo ) 法近似計(jì)算法近似計(jì)算n=100000; %總模

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

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