2019年數(shù)值積分的matlab實(shí)現(xiàn)

1、f(x9h)?f(x?h)?oo?(x)f o2h?3y?4y?y?2io?)(fx y?4y?3y?Xmm“?x)f(三點(diǎn)公式。其誤差均為102數(shù)值微分的MATLAB實(shí)現(xiàn)實(shí)驗(yàn)10數(shù)值積分實(shí)驗(yàn)?zāi)康模?. 了解數(shù)值積分的基本原理：2. 熟練掌握數(shù)值積分的MATLAB實(shí)現(xiàn)：3. 會(huì)用數(shù)值積分方法解決-些實(shí)際問(wèn)題。實(shí)驗(yàn)內(nèi)容：積分是數(shù)學(xué)中的個(gè)基本概念，在實(shí)際問(wèn)題中也有很廣泛的應(yīng)用。同微分樣，在微積分中， 它也是通過(guò)極限定義的，由于實(shí)際問(wèn)題中遇到的函數(shù)般都以列農(nóng)形式給岀，所以常常不能用來(lái) 宜接進(jìn)行積分。此外有些函數(shù)雖然有解析式，但其原函數(shù)不是初等函數(shù)，所SlnX 1?Xd這時(shí)我們般考慮用數(shù)值方法計(jì)算其如

2、不定積分。以仍然得不到積分的精確值， Xo近似值，稱(chēng)為數(shù)值積分。10. 1數(shù)值微分簡(jiǎn)介x)f(x?y在設(shè)函數(shù)可導(dǎo)，則其導(dǎo)數(shù)為-)x?f(x(?h)f*?lim?(Xf) )(10. 1 homy?f(x)以列農(nóng)形式給出(見(jiàn)表10-1),則其精確值無(wú)法求得，但可由 下式求得如果函數(shù)其近似值)(fk?h)f(x?xf)( 10. 2)(h 表 10-1Xyhy?f(x)在(10.2)式右端項(xiàng)的分了稱(chēng)為函數(shù)般的，步長(zhǎng)越小，所得結(jié)果越精確C= xx的差分， 所以右端項(xiàng)又稱(chēng)為差商。的差分，分母稱(chēng)為自變量在數(shù)值微分即用差商近似代替微商。常用的差 商公式為：(10.3 ) (10.4)oh2.(10.5)

3、n21)h(,稱(chēng)為統(tǒng)稱(chēng)MATLAB提供了 個(gè)指令求解階向前差分，其使用格式為: dx=diff(x)?nx?,XXX,x?X?l?n,這樣展于兩點(diǎn)的數(shù)值導(dǎo)dx維數(shù)組，為其中X是維數(shù)組】述1數(shù)可通過(guò)指令diff(x)h實(shí)現(xiàn)。對(duì)于三點(diǎn)公式，讀者可參考例1的M函數(shù)文件diff3.mox7f(xx)yf的值由下衣給用三點(diǎn)公式計(jì)算處的導(dǎo)數(shù)值，在1. 0,1. 2,1. 41例X)(xf1.01.31. 40. 25000. 226S0. 20660. 1S900.1736解:建立三點(diǎn)公式的M函數(shù)文件diff3.m如下:Rinction f=df(x,y)n=length(x) ;h=x(2)-x( 1)

4、;f(lH-3y(l)÷4*y(2)-y(3)(2*h);forj=2:n-lf(j)=(y(ji)-y(j-i)y(2*h);endf(nXy(n-2)-4y(n4)÷3*y(n)(2*h);在MATLAB指令窗中輸入指令：x=1Q1.1,1.2,L3,14; y=0.2500,0.2268,0.2066,0.1890,0.1736;diff3(%y)y?f(x)所以，-0.0014。，運(yùn)行得各點(diǎn)的導(dǎo)數(shù)值為：-0.2470, -0.2170-0. 1890,-0. 1650x71. 0, 1. 2, 1. 4 處的導(dǎo)數(shù)值分別為-0 2470> -0. 1890 和-

5、0. 0014 在。對(duì)于高階導(dǎo)數(shù)，MATLAB提供了幾個(gè)指令借助于樣條函數(shù)進(jìn)行求導(dǎo)，詳細(xì)使用步驟如下：stepl：對(duì)給定數(shù)據(jù)點(diǎn)(x, y),利用指令PP=SPIine(x, y),獲得三次樣條函數(shù)數(shù)據(jù)pp,供后面 PPVaI等指令使用。其中，PP是個(gè)分段多項(xiàng)式所對(duì)應(yīng)的行向量，它包含此多項(xiàng)式的階數(shù)、段數(shù)、 節(jié)點(diǎn)的橫坐標(biāo)值和各段多項(xiàng)式的系數(shù)。step2:對(duì)于上而所求的數(shù)據(jù)向 PP»利用指令breaks, coefs, m, n=unmkpp(pp)進(jìn)行處理，生 成幾個(gè)有序的分段多項(xiàng)式PP。step3:對(duì)各個(gè)分段多項(xiàng)式PP的系數(shù)，利用函數(shù)PPVaI生成其相應(yīng)導(dǎo)數(shù)分段多項(xiàng)式的系數(shù)，再利 用指

6、令mkpp生成相應(yīng)的導(dǎo)數(shù)分段多項(xiàng)式step4:將待求點(diǎn)代入此導(dǎo)數(shù)多項(xiàng)式，即得樣條導(dǎo)數(shù)值。上述過(guò)程可建立M函數(shù)文件ppd. m實(shí)現(xiàn)如下：ficton dy=d(p)breaks.coefs,m=mmk();for=l :mCOefSm(ij)=olyder(coefs(,:);enddy=nk(breaks,coefsm);于是，如果已知節(jié)點(diǎn)處的值x,y,可用下面指令計(jì)算XX處的導(dǎo)數(shù)dyy：PP=SPIine(X,y).dr=ppd(pp):dy3r=ppval(dy.xx);y?SlnX的數(shù)據(jù)點(diǎn)，利用三點(diǎn)公式和三次樣條插值分別求導(dǎo)，并與基于正弦函數(shù)例2解析所求得 的導(dǎo)數(shù)進(jìn)行比較。解:編寫(xiě)M腳本

7、文件bi jiao, m如下：h=0.1*i=0*iir=sm(x); dyl=d3(x,y); PP=SPlme(X,y)：dy=ppd(pp):dy2=ppval(dy.x); Z=COS(x);error 1 =norm(dy 1 -z).error2=norm(dy2-z) )rb,r-,.z.lot(x,dyl,k x,dy2,運(yùn)行得結(jié)果為:errorl =0. 0666, error2 =0. 0025,生成圖形見(jiàn)圖10IO0.5 -0 -0.5 -1o三點(diǎn)公式、三次樣條插值與解析求詁比較圖圖10.1顯然利用三次樣條插值求導(dǎo)所得謀差比三點(diǎn)公式求導(dǎo)小很多，同時(shí)由圖2. 15可知利用三

8、次樣條 插值求導(dǎo)所得曲線(xiàn)與解析求導(dǎo)曲線(xiàn)基本重合,而三點(diǎn)公式在極值點(diǎn)附近和兩個(gè)端點(diǎn)附近謀差較人, 其它點(diǎn)吻合的較好。10.3應(yīng)用示例：湖水溫度變化問(wèn)題問(wèn)題：湖水在夏天會(huì)出現(xiàn)分層現(xiàn)象，其特點(diǎn)是接近湖面的水的溫度較高，越往下水的溫度越低。 這種現(xiàn)象會(huì)影響水的對(duì)流和混合過(guò)程，使得下層水域缺氧，導(dǎo)致水生魚(yú)類(lèi)死亡。對(duì)某個(gè)湖的水溫 進(jìn)行觀(guān)測(cè)得數(shù)據(jù)見(jiàn)表10-2 o表10-2某湖的水溫觀(guān)測(cè)數(shù)據(jù)深度()02.34.99.113.71S. 322.927.2溫度(C)22. S22. S22. S20.613.911.711. 111.1試找出湖水溫度變化最人的深度。1. 問(wèn)題的分析湖水的溫度可視為關(guān)于深度的函數(shù)，

9、于是湖水溫度的變化問(wèn)題便轉(zhuǎn)化為溫度函數(shù)的導(dǎo)數(shù)問(wèn)題，顯 然導(dǎo)函數(shù)的最人絕對(duì)值所對(duì)應(yīng)的深度即為溫度變化最人的深度。對(duì)于給定的數(shù)據(jù)，但考慮到所給 從而得到溫度變化最大的深度，可以利用數(shù)值微分計(jì)算各深度的溫度變化值，的數(shù)據(jù)較少，由此計(jì)算的深度不夠精確，所以采用插值的方法計(jì)算加密深度數(shù)據(jù)的導(dǎo)數(shù)值，以得 到更準(zhǔn)確的結(jié)果。2. 模型的建立及求解TT7T(h)T(h)可(°C),相應(yīng)的溫度為，且有記湖水的深度為，并假定函數(shù)(m)導(dǎo)。T(h)的插值導(dǎo)函數(shù)：然后將給定對(duì)給定的數(shù)據(jù)進(jìn)行三次樣條插值，并對(duì)其求導(dǎo)，得到的深度數(shù)據(jù) 加密,捜索加密數(shù)據(jù)的導(dǎo)數(shù)值的絕對(duì)值，找出其最人值及其相應(yīng)的深度，相應(yīng)的MATLA

10、B指令如 下：h=0 2.3 4.9 9.1 13.7 18.3 22.9 27.2;T=22.8 22.8 22.8 20.6 13.9 11.7 11.1 11.1;hh=O:0.1:27.2;PP=SPIme(h.T):dT=ppd(pp):dTT=ppval(dT.hli);dTTmax.=max(abs(dTT),hh()').grd Onplot(hh,dTThh().dTT(),b,r.h= 11. 4,最人值為運(yùn)行得導(dǎo)函數(shù)絕對(duì)值的最人值點(diǎn)為：1.6139,即湖水 在深度為11.4m時(shí)溫度變化最人，如圖10.2所示(黑點(diǎn)為溫度變化最人的點(diǎn)。0.5圖10. 2湖水沮度變化曲

11、線(xiàn)圖10.4數(shù)值積分簡(jiǎn)介考慮定積分b?f(x)dx(10. 6) df(x)是以列衣形式給出，則其求解思想同數(shù)值微分類(lèi)似，即用逼近多項(xiàng)如果彼積函數(shù)b?P(x)dx)P(x)(xf, 得(式10.6,然后計(jì)算積分近似地代替被積函數(shù))式的近似值；皿如果彼積函數(shù)的原函數(shù)不是 初等函數(shù)，則將積分區(qū)間進(jìn)行細(xì)分，對(duì)每個(gè)小區(qū)間，用個(gè)近似f(x),然后積分得(10.6)式的 近似值。這兩種類(lèi)型最終都可歸結(jié)為函函數(shù)代替被積函數(shù)Xf(X)f(x的某種線(xiàn)性組合，即下面數(shù) 值求積公式：在節(jié)點(diǎn)上的函數(shù)值數(shù)Iinb?)(XX(X)d?fAIf ) 10.7或(kka"knb?PRAf(x)?f(x)dx)10.

12、8(IiaO匕?fR 為截?cái)嘀\差。此謀差可用代數(shù)精度衡量，代數(shù)精度越高，謀差越?。悍粗`其中差越人。m)(xf的是個(gè)次數(shù)不超過(guò)代數(shù)精度是用來(lái)衡量數(shù)值積分公式近似程度的辦法，如果lm?)f(x)式不能精是 個(gè)(代數(shù)多項(xiàng)式，(10. 7)式等號(hào)成立：而當(dāng)10.7次多項(xiàng)式時(shí)，m。確成立，則稱(chēng)(10. 7) 式的代數(shù)精度為選取不同的近似函數(shù)，可產(chǎn)生不同的數(shù)值求積公式，常見(jiàn)的有：梯形公式、辛普 森公式和高斯公式。MATLAB實(shí)現(xiàn)10. 5數(shù)值積分的MATLAB捉供了下面幾個(gè)函數(shù)計(jì)算積分，其使用格式分別為：)nO,lf(klh?»?,x) 1 O trapz(X)采用梯形公式計(jì)算積分(k法計(jì)算積

13、分)quad('fun', a, b, tol)采用自適應(yīng)SimPSOn (2法計(jì)算積分)QUadI (, fun', a, bf tol)釆用自適應(yīng)GaUSS-LObattO (3 a, b為積分的上、下限。其中fun為被積函數(shù)：tol是可選項(xiàng)，農(nóng)示絕對(duì)謀差,i2?x?XdL并與其法計(jì)算GaUSS-LObattO 例1分別利用梯形公式、SilnPSOn公式和o精確值比較。解：先對(duì)積分作符號(hào)運(yùn)算，然后將其計(jì)算結(jié)果轉(zhuǎn)換為數(shù)值型，再將其與這三種方 法求得的數(shù)值解比較，其MATLAB指令為：SymSXXZO=Suile(int(,sqrt(l+xxA2)0.1)z=doubl

14、e(z0);Z=VPQ(Z,8)X=OrOOl :1; y=sqrt(l+x 人2);zl=traz(y)*0.0Ezl=rpa(zL8).el=z-zlerrl=,a(errl,8) z2=quad('sqrt(l+xA2)0.1);z2=pa(z2,8),err2=z-z2;err2=rpa(err2.8) z3=quadlCsqft(l+x.A2y,0,l);z3=pa(z3,8Xen3=z-z3;err3=pa(err3,8)l(2?ln(2?l)?l. 1477936,三種公式計(jì)算得數(shù)值積分值分別為運(yùn)行得精確值為一 21.1477995,1. 1477935和1. 14779

15、36,其相應(yīng)誤差分別為-59e5, . le-6和0.,由三者誤差可見(jiàn), GaUSS-LObattO法計(jì)算最為持確，SimPSOn公式次之，梯形公式最差，但它也能持確到小數(shù)點(diǎn)后 5位數(shù)。例2人造地球衛(wèi)星軌道可視為平面上的橢圓。我國(guó)第顆人造地球衛(wèi)星近地點(diǎn)距地球農(nóng)面439km, 遠(yuǎn)地點(diǎn)距地球農(nóng)面2384km,地球半徑為6371km,求該衛(wèi)星的軌道長(zhǎng)度。?),a,2(0t?t?bsm?,cos?Xatyb分別是長(zhǎng)、衛(wèi)星軌道橢圓的參數(shù)方程為解:短半ab=6371+439=68IOo 軸，則根據(jù)所給數(shù)據(jù)知=6371+2384=8755,由對(duì)弧長(zhǎng)的曲線(xiàn)積分知識(shí)知，橢圓的長(zhǎng)度為?12222?tdcostL?

16、4)(asmt?b22o上積分稱(chēng)為橢圓積分，它無(wú)法用解析方法計(jì)算，可用計(jì)算其數(shù)值解，編寫(xiě)H函數(shù)文件如下：Rinction y=y(t)a=8755:b=6810;)'=4*sqrt(aA2*sin(t).A2-bA2*cos(t).A2);在MATLAB指令窗中輸入以下指令：l=quad('y'.0.pi2)運(yùn)行得結(jié)果為：1=4. 9090e+004即衛(wèi)星軌道長(zhǎng)度為4909OknlO對(duì)于用列農(nóng)形式給出的函數(shù)，上述方法不再適用，可利用指令SPIine構(gòu)造三次樣條插值函數(shù)， 再計(jì)算積分，具體步驟可參考例2。例3在橋梁的端每隔段時(shí)間記錄Imin有幾輛車(chē)過(guò)橋，得到數(shù)據(jù)見(jiàn)<

17、 10-3：表10-3過(guò)橋車(chē)輛數(shù)據(jù)時(shí)間0:002:004:005:006:007:008:00輛車(chē)輛數(shù)/29025S25時(shí)間9:0010:3011:3012:3014:0016:0017:00輛車(chē)輛數(shù)/12 O10127928時(shí)間18:0019:0020:0021:0022:0023:0021:00輛車(chē)輛數(shù)/10911893試估計(jì)天通過(guò)橋梁的車(chē)流量。24?tx(t)dt)(tx,解：記記錄時(shí)刻為時(shí)，相應(yīng)的車(chē)輛數(shù)為，則所求車(chē)流量即為計(jì)算積分0MATLAB 則在指令窗中輸入下面指令:x=02,4,5,6,7,8,9,103,113,123,14,16,17,18,19,20,21,22,2324；

18、 y=22025,8,25,12,5,10.12,7.9,2&22,10911,893;% 利用三次樣條插值計(jì)算積分 PP=SPIme(X.y)isl=quadl(fun0,24,.)s2=trapz(x,y)%利用梯形公式計(jì)算枳分其中H函數(shù)文件fun.m為：fiction XFf=ftn(x.)Vrf=PPVai(PP.x);運(yùn)行得三次樣條插值計(jì)算所得車(chē)流量為212輛，梯形公式計(jì)算所得車(chē)流量為216輛。10.6數(shù)值積分的建模示例：煤炭?jī)?chǔ)量計(jì)算問(wèn)題問(wèn)題:某煤礦為估計(jì)其煤炭的儲(chǔ)量在該礦區(qū)內(nèi)進(jìn)行勘探,得到數(shù)據(jù)如衣10-4試估l?x",l?y?5) 煤炭的儲(chǔ)量。算出該礦區(qū)(表IoT

19、某煤礦勘探數(shù)據(jù)表編號(hào)193467S910X坐標(biāo)(km1111122222y坐標(biāo)(km1234 O1234511 ) (m煤炭師度13. 7225. 80S. 4725. 2799 3°15.4721.3311. 4924. S326. 19編號(hào)11121314151617IS1920X ) km坐標(biāo)(3333344441y坐標(biāo)(km12341231511 (In)煤炭伸度23. 2826. 4S29. 1412. 0414. 5819.9523. 7315. 351S.0116. 291.問(wèn)題的分析問(wèn)題給出了很多點(diǎn)對(duì)應(yīng)的煤炭厚度，顯然整個(gè)煤礦可以看作是個(gè)巨人的曲頂柱體(見(jiàn)圖10.3,

20、 此圖經(jīng)過(guò)插值得到)，而煤炭的儲(chǔ)雖即為此立體的體積。要計(jì)算此立體的體積，可以利用插值得 到曲頂柱體的頂而函數(shù)，再對(duì)其積分；也可以將此曲頂柱體分割成若干個(gè)細(xì)的曲頂柱體，用數(shù)值 方法計(jì)算這些細(xì)曲頂柱體的體積，再對(duì)其求和即得原曲頂柱體的體積。1000 1000圖10.3煤炭M度曲而圖2.模型的建立及求解Zli.則它是坐標(biāo)坐標(biāo)建立空間坐標(biāo)系。記煤炭的厚度為以煤炭的厚度為三維空間中的 ?W)y(xZEyX,為的二元函數(shù)，即，則由二重積分的知識(shí)可知，此煤礦的煤炭?jī)?chǔ)量?(x，y)dxdyW?(10.9)D?|D?(x,y)l?x?4,ly?5其中。7)yx,(H給出了 些離散點(diǎn)上的函數(shù)值，無(wú)法直接計(jì)算上述二

21、重 積分，所以由于函數(shù)下面采用數(shù)值積分的方法計(jì)算其值。由數(shù)值積分的知識(shí)知，計(jì)算定積分有復(fù)合梯形公式為 nnhb?f(x)(b)?2df(x)xf(a)f(10. 10)k2akn?)n 1 ,x(k01ikhxa»?O為步長(zhǎng)，其中為節(jié)點(diǎn)，且有吐由(10. 9)式得Mb? ? ? ? 7xdydx9g(x)d? WyX(J (10.11)?®a?ydy)(xg(x)?,)式可得其中，則由(10. IOehIQ'b?xg(x?2)(a?)?g(b)g(x)dx?gW(10. 12)j2aj?l 而hm“!d?y?,)?2y)(g(a)?(a,c)?a(a,d(a,y)d

22、y?一I叱krid?y?2y)(b,(b,c)?(b,dg(b)?)bCy)dy? _k2cik?hm?id?y? ?)(x,)?y(x,dg(x)?)(x,y)dy?2x(,cu±kkk2ci?k 所以有hhnri?Y?)yb)(a,y)?(b,c)?(b0)W?2Q(a,c?(aqMiwi加？?y?2)(xc(xj?(Xq)2?頃ij( 10. 13) hli?(b,dc)?5)(a,d?)?5(b,?5(a,c)?5 43?區(qū)丫)(人(1)勺(b,)?24(x,c)?2?y(a,)呻Jg 考慮到給定的數(shù)據(jù)較少，由此產(chǎn)生的誤差較人，所以利用插值后的數(shù)據(jù)計(jì)算(10. 13)式，相應(yīng)

23、的MATLAB計(jì)算指令 如下：x=l 111122222333334444 4n000;y=l 234512345123451234 5*1000:z=-13.72 25.80 8.47 25.27 22.32 15.47 21.33 14.49 24.83 26.19 23.28 26.48 29.14 12.04 14.5819.95 23.73 1535 18.01 16.29;hx=10:hy=10;CX=IoO(xW000;Cy=IoOOhy5000;X, Y= meShgnd(CXXy);II=Iength(Cx);m=length(cy);Z=gndd 也(xy,zXY,PV);%插值Surf(X YZ)% 繪制圖 10.3W=-hx*hy*(-5*Z(l,l)-5*Z(l,n)-5*Z(md)-5*Z(m.n)+2*(suni(Z(l,:)*Z(mj)+sum(Z(:,l)+Z(:>n)+4*sum(sum(Z) ”4w：s.ml,即煤礦的煤炭?jī)?chǔ)雖約為2. 5242X10運(yùn)行得X=N 5242實(shí)驗(yàn)任務(wù)：D?D(t)如下農(nóng)所示：1. 一個(gè)物體的運(yùn)動(dòng)距離tS.09.010.011.012.0OtD1