數(shù)值分析第8次課

1、1問(wèn)題問(wèn)題: 給定給定n個(gè)數(shù)據(jù)點(diǎn)個(gè)數(shù)據(jù)點(diǎn) (xi , yi ) (i=1, 2, , n)ix2x1xnxiy1y2yny求直線求直線 y=a0+a1x 使得使得 niiixaay1102)(達(dá)到最小達(dá)到最小. 最小二乘一次多項(xiàng)式擬合最小二乘一次多項(xiàng)式擬合1 數(shù)據(jù)擬合的最小二乘法數(shù)據(jù)擬合的最小二乘法2 , , , , 21201011112xxxyxxyxxxxxyyxxxyyxxxxniiyiniixyniixniixxSnSSSnSaSnSSSSSaSSaaSSSnySyxSxSxS 設(shè)設(shè) 由上式求得由上式求得a0, a1, 代入代入 y=a0+a1x 得到得到最小二乘擬最小二乘擬合合(直

2、線直線)一次多項(xiàng)式一次多項(xiàng)式.3問(wèn)題問(wèn)題: 給定給定n個(gè)數(shù)據(jù)點(diǎn)個(gè)數(shù)據(jù)點(diǎn) (xi , yi ) (i=1, 2, , n)ix2x1xnxiy1y2yny求求,2210 xaxaay 使得使得 niiiixaxaay122102)(達(dá)到最小達(dá)到最小. 最小二乘二次多項(xiàng)式擬合最小二乘二次多項(xiàng)式擬合4 niiiniiiniiniiniiniiniiniiniiniiniiyxyxyaaaxxxxxxxxn121121014131213121121 用用 Cholesky分解法求此對(duì)稱正定陣分解法求此對(duì)稱正定陣 用用 MATLAB 函數(shù)函數(shù) z = Ar 由上式求得由上式求得a0, a1, a2, 得

3、到得到最小二乘擬合二次多項(xiàng)式最小二乘擬合二次多項(xiàng)式正則方程組正則方程組5xiyi例例 給定一組實(shí)驗(yàn)數(shù)據(jù)如下給定一組實(shí)驗(yàn)數(shù)據(jù)如下求求x, y的函數(shù)關(guān)系的函數(shù)關(guān)系.1 2 3 4 6 7 82 3 6 7 5 3 2u (1) 作散點(diǎn)分布圖作散點(diǎn)分布圖點(diǎn)的分布近似為拋物線點(diǎn)的分布近似為拋物線6 (2)確定近似表達(dá)式確定近似表達(dá)式設(shè)擬合曲線為二次多項(xiàng)式設(shè)擬合曲線為二次多項(xiàng)式22102)(xaxaaxPy (3) 建立正則方程組建立正則方程組, 7 n,3171 iix,179712 iix,1171713 iix,8147714 iix,2871 iiy,12171 iiiyx,635712 iii

4、xy7故正則方程組為故正則方程組為 6351212881471171179117117931179317210aaa (4) 求解正則方程組得求解正則方程組得,3864. 0,4318. 3,3182. 1310 aaa故所求擬合曲線為故所求擬合曲線為.3864. 04318. 33182. 1)(22xxxPy 8xiyi例例 給定一組實(shí)驗(yàn)數(shù)據(jù)如下給定一組實(shí)驗(yàn)數(shù)據(jù)如下求求x, y的函數(shù)關(guān)系的函數(shù)關(guān)系.1 2 3 4 6 7 82 3 6 7 5 3 2 Matlab解法解法: polyfit(1, 2, 3, 4, 6, 7, 8, 2, 3, 6, 7, 5, 3, 2, 2) ans=

5、-0.3864 3.4318 -1.31829例例 測(cè)得一發(fā)射源的發(fā)射強(qiáng)度測(cè)得一發(fā)射源的發(fā)射強(qiáng)度 I 與時(shí)間與時(shí)間 t 的一組數(shù)據(jù)如下的一組數(shù)據(jù)如下tiIi0.2 0.3 0.4 0.5 0.6 0.7 0.83.16 2.38 1.75 1.34 1.00 0.74 0.56試用最小二乘法確定試用最小二乘法確定 I 與與 t 的函數(shù)關(guān)系的函數(shù)關(guān)系.tIn (1)作散點(diǎn)分布圖作散點(diǎn)分布圖可以考慮用指數(shù)函數(shù)近似可以考慮用指數(shù)函數(shù)近似10 列數(shù)據(jù)表列數(shù)據(jù)表tiIi0.2 0.3 0.4 0.5 0.6 0.7 0.83.16 2.38 1.75 1.34 1.00 0.74 0.56lnIi1.1

6、506 0.8671 0.5596 0.2927 0 0.3011 0.5798求求lnI與與t的最小二乘直線的最小二乘直線. 將上表數(shù)據(jù)代入正則方程組將上表數(shù)據(jù)代入正則方程組得得 1858. 09891. 103. 25 . 35 . 3710aa其解為其解為.89. 2,73. 110 aa故所求擬合曲線為故所求擬合曲線為.64. 589. 289. 273. 1lnttIeeeI 11 列數(shù)據(jù)表列數(shù)據(jù)表tiIi0.2 0.3 0.4 0.5 0.6 0.7 0.83.16 2.38 1.75 1.34 1.00 0.74 0.56lnIi1.1506 0.8671 0.5596 0.29

7、27 0 0.3011 0.5798求求lnI與與t的最小二乘直線的最小二乘直線. 將上表數(shù)據(jù)代入正則方程組將上表數(shù)據(jù)代入正則方程組得得 1858. 09891. 103. 25 . 35 . 3710aa其解為其解為.89. 2,73. 110 aa故所求擬合曲線為故所求擬合曲線為.64. 589. 289. 273. 1lnttIeeeI Matlab解法解法:polyfit(0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 1.1506, 0.8671, 0.5596, 0.2927, 0, -0.3011, -0.5798, 1) ans= -2.8883 1.72

8、8312 求數(shù)據(jù)組的最小二乘擬合函數(shù)的步驟求數(shù)據(jù)組的最小二乘擬合函數(shù)的步驟 (1) 由給定數(shù)據(jù)確定近似函數(shù)的表達(dá)式由給定數(shù)據(jù)確定近似函數(shù)的表達(dá)式, 一般可一般可通過(guò)描點(diǎn)觀察或經(jīng)驗(yàn)估計(jì)得到通過(guò)描點(diǎn)觀察或經(jīng)驗(yàn)估計(jì)得到. (2) 按最小二乘原則確定表達(dá)式中的參數(shù)按最小二乘原則確定表達(dá)式中的參數(shù), 即由即由殘差平方和最小導(dǎo)出正則方程組殘差平方和最小導(dǎo)出正則方程組, 求解得參數(shù)求解得參數(shù).13第七章第七章 數(shù)值微分與數(shù)值積分?jǐn)?shù)值微分與數(shù)值積分1 數(shù)值微分?jǐn)?shù)值微分2 NewtonCotes求積公式求積公式3 復(fù)化求積公式復(fù)化求積公式4 Romberg求積公式求積公式14 利用離散點(diǎn)上函數(shù)的信息求函數(shù)導(dǎo)數(shù)近

9、似值利用離散點(diǎn)上函數(shù)的信息求函數(shù)導(dǎo)數(shù)近似值的方法的方法, 稱為稱為數(shù)值微分?jǐn)?shù)值微分. 差商型數(shù)值微分公式差商型數(shù)值微分公式1 數(shù)值微分?jǐn)?shù)值微分15由導(dǎo)數(shù)定義由導(dǎo)數(shù)定義hxfhxfxfh)()(lim)(0 當(dāng)當(dāng)h很小時(shí)很小時(shí), 可用可用差商差商近似導(dǎo)數(shù)近似導(dǎo)數(shù).16 差商型求導(dǎo)公式差商型求導(dǎo)公式 (3)中心差商公式中心差商公式0,)()()( hhxfhxfxf,)()()(hhxfxfxf .2)()()(hhxfhxfxf (1) 向前差商公式向前差商公式 (2) 向后差商公式向后差商公式17 幾何意義幾何意義hx xhx ABChxfhxfkBC)()( hhxfxfkAB)()( hh

10、xfhxfkAC2)()( B點(diǎn)切線斜率點(diǎn)切線斜率)(xf 從幾何直觀看從幾何直觀看: 中心差商效果最好中心差商效果最好18 截?cái)嗾`差截?cái)嗾`差)(2)( )()()( 1hOhhxfhxfhxfxf )(2)( )()()( 2hOhhxfhhxfxfxf )(12)()(2)()()( 223)3(3)3(hOhhxfhxfhhxfhxfxf 其中其中1,0321 由由Taylor公式可得公式可得19 二階導(dǎo)數(shù)的中心差商公式二階導(dǎo)數(shù)的中心差商公式2)()(2)()( hhxfxfhxfxf 截?cái)嗾`差截?cái)嗾`差)(12)()(2)()( )4(22 fhhhxfxfhxfxf 20近似計(jì)算近似計(jì)

11、算 badxxfI)(數(shù)值積分?jǐn)?shù)值積分21 依據(jù)微積分基本定理依據(jù)微積分基本定理, 只要找到被積函數(shù)只要找到被積函數(shù) f (x)的的原函數(shù)原函數(shù) F (x), F (x)=f (x), 便有便有)()()(aFbFdxxfba 為什么還要對(duì)積分進(jìn)行近似計(jì)算為什么還要對(duì)積分進(jìn)行近似計(jì)算 大量的被積函數(shù)找不到用初等函數(shù)表示的原函數(shù)大量的被積函數(shù)找不到用初等函數(shù)表示的原函數(shù) 實(shí)驗(yàn)測(cè)量或數(shù)值計(jì)算給出的通常是一張離散函數(shù)實(shí)驗(yàn)測(cè)量或數(shù)值計(jì)算給出的通常是一張離散函數(shù)表表, 即被積函數(shù)的表達(dá)式未知即被積函數(shù)的表達(dá)式未知.數(shù)值積分?jǐn)?shù)值積分222 Newton-Cotes 公式公式基本思想基本思想: 利用利用插值

12、多項(xiàng)式插值多項(xiàng)式).()(xfxLn banbadxxLdxxffI.)()()(其中其中Ln(x)是是n階階Lagrange插值多項(xiàng)式，用插值多項(xiàng)式，用Ln (x)的的積分近似積分近似 f (x)的積分，即的積分，即插值型求積公式插值型求積公式23 bankjjjkjkdxxxxxA,0)()(由由 決定決定,與與 無(wú)關(guān)無(wú)關(guān).節(jié)點(diǎn)節(jié)點(diǎn) f (x) 在在a, b上取上取 a x0 x1 0, 使得使得|,|max0knk 則稱該求積公式是則稱該求積公式是穩(wěn)定穩(wěn)定的的. 求積公式的穩(wěn)定性求積公式的穩(wěn)定性48 若求積公式是穩(wěn)定的若求積公式是穩(wěn)定的, 則則 f (x)的觀察值的較小的的觀察值的較小的

13、誤差引起的求積結(jié)果的誤差也是較小的誤差引起的求積結(jié)果的誤差也是較小的. 求積公式求積公式?jīng)]有把沒(méi)有把 f (x)的誤差的誤差“放大放大”很多很多.49), 1 , 0()(nkfxfkk nkkkknnfxfAfIfI0)()()(證明證明因此復(fù)化梯形公式是數(shù)值穩(wěn)定的因此復(fù)化梯形公式是數(shù)值穩(wěn)定的.當(dāng)當(dāng)).(2)1(2abhhnh 定理定理 復(fù)化梯形公式是復(fù)化梯形公式是數(shù)值穩(wěn)定數(shù)值穩(wěn)定的的. nkkkkfxfA0)( nkkA0 50 x0 x2xf (x)x4hhxn 2hxnmnnabh2, .hx3x1xn 1 復(fù)化復(fù)化Simpson公式公式分片二次多項(xiàng)式近似分片二次多項(xiàng)式近似51 將積分

14、區(qū)間將積分區(qū)間a, b劃分為劃分為n=2m等分等分, 步長(zhǎng)步長(zhǎng) h=( b a )/n, 分點(diǎn)分點(diǎn) xk= a+kh ( k=0, 1, , n). 在每個(gè)在每個(gè)小區(qū)間小區(qū)間 x2k 2 , x 2k ( k=1, , m)上用上用Simpson公式：公式： kkxxkkkkkxfxfxfxxdxxf222)()(4)(6)(21222222 )()(4)(321222kkkxfxfxfh 復(fù)化復(fù)化Simpson公式公式k=1, , m52 111122)(4)(2)()(3mkmkkkxfxfbfafh= Sn( f ) mkxxbadxxfdxxffIkk1222)()()( )()(4)

15、(3212221kkkmkxfxfxfh 復(fù)化復(fù)化Simpson公式公式53)()()(fSfIfRnn 當(dāng)當(dāng) f (x)在在a, b上具有四階連續(xù)導(dǎo)數(shù)時(shí)上具有四階連續(xù)導(dǎo)數(shù)時(shí), ),(),()(1)4()4(bamffmkk ),(222kkkxx 故得故得),(),(180)()(90)()4(4)4(5bafhabfmhfRn 復(fù)化復(fù)化Simpson公式的截?cái)嗾`差公式的截?cái)嗾`差 mkkfh1)4(5),(90 54 由復(fù)化由復(fù)化Simpson公式的截?cái)嗾`差知公式的截?cái)嗾`差知, 誤差階為誤差階為 h4, 收斂性是顯然的收斂性是顯然的, 事實(shí)上事實(shí)上, 只要只要 f (x) Ca, b則則可得

16、到可得到收斂性收斂性, 即即.)()(limdxxffSbann 由于求積系數(shù)均為正由于求積系數(shù)均為正, 與復(fù)化梯形公式一樣的與復(fù)化梯形公式一樣的證法可得復(fù)化證法可得復(fù)化 Simpson公式是公式是數(shù)值穩(wěn)定數(shù)值穩(wěn)定的的.55例例 對(duì)于函數(shù)對(duì)于函數(shù),sin)(xxxf 給出給出n=8的函數(shù)表的函數(shù)表, 試用試用復(fù)化梯形公式及復(fù)化復(fù)化梯形公式及復(fù)化Simpson公式計(jì)算積分公式計(jì)算積分 10.sindxxxI解解ix)(ixf0.00.1250.250.3750.50.6250.750.8751.01.00.99739780.98961580.97672670.95885100.93615560.

17、90885160.84147090.8771925應(yīng)用復(fù)化梯形公式求得應(yīng)用復(fù)化梯形公式求得T8=0.9456909應(yīng)用復(fù)化應(yīng)用復(fù)化Simpson公式求得公式求得S8=0.9460832準(zhǔn)確值準(zhǔn)確值 I=0.9460831兩者運(yùn)算量基本相同兩者運(yùn)算量基本相同56trapz: 復(fù)化梯形公式求積分復(fù)化梯形公式求積分.用法用法: trapz(X, Y), 其中其中X, Y為相同維數(shù)的向量為相同維數(shù)的向量.例例: X=0.125:0.125:1.0;Y=sin(X)./X;X=0,X;Y=1,Y;trapz(X,Y)ans = 0.94569086358270Matlab函數(shù)函數(shù)57例例 若用復(fù)化求積公

18、式計(jì)算積分若用復(fù)化求積公式計(jì)算積分dxeIx 10的近似值的近似值, 若要求計(jì)算結(jié)果有若要求計(jì)算結(jié)果有4位有效數(shù)字位有效數(shù)字, n應(yīng)取多大應(yīng)取多大?解解, 1110 dxeIex.105 . 04 1 , 0, 1| )(|)( xexfxk復(fù)化梯形公式的誤差復(fù)化梯形公式的誤差)( )(12|2 fabhRT .83.40 n若用復(fù)化梯形公式求積分若用復(fù)化梯形公式求積分, n取取41能達(dá)到精度要求能達(dá)到精度要求.2121n 41021 58故應(yīng)取故應(yīng)取n=4. 該例表明該例表明, 為達(dá)到相同的精度為達(dá)到相同的精度, 用復(fù)化用復(fù)化Simpson公式所需的計(jì)算量比復(fù)化梯形公式要少公式所需的計(jì)算量比

19、復(fù)化梯形公式要少, 這也說(shuō)明這也說(shuō)明了復(fù)化了復(fù)化Simpson公式的精度高公式的精度高.復(fù)化復(fù)化Simpson公式的誤差公式的誤差)()(180|)4(4 fabhRS .25. 3 n41801n 41021 59 復(fù)化復(fù)化梯形公式的逐次分半算法梯形公式的逐次分半算法將區(qū)間將區(qū)間a, b分成分成n=2m等分等分, 記記,2mmabh .)(2)()(21212 mkmmkhafbfafhTm), 2 , 1 , 0( m稱稱 為為梯形值序列梯形值序列.2mT60 mhTTmm2122所有新增加節(jié)點(diǎn)的函數(shù)值之和所有新增加節(jié)點(diǎn)的函數(shù)值之和其中其中.2mmabh 復(fù)化復(fù)化梯形公式的逐次分半算法梯形

20、公式的逐次分半算法61以以n=8, m=3為例為例. 記記 fk= f (xk)x0 x2x4x6x3x1x5x7x8 )(2167654321808fffffffffabT )(21664280fffffab 75318ffffab 24T 所有新增加節(jié)點(diǎn)的函數(shù)值之和所有新增加節(jié)點(diǎn)的函數(shù)值之和. 3h62 復(fù)化梯形公式余項(xiàng)的后驗(yàn)估計(jì)復(fù)化梯形公式余項(xiàng)的后驗(yàn)估計(jì));,(),(12112bafhabTIn );,(),(2122222bafhabTIn f ( 1 ), f ( 2 ) 分別是分別是 f (x) 在在a, b上的上的n個(gè)點(diǎn)與個(gè)點(diǎn)與 2n 個(gè)點(diǎn)處的算術(shù)平均值個(gè)點(diǎn)處的算術(shù)平均值 (每個(gè)

21、小區(qū)間上取一個(gè)點(diǎn)每個(gè)小區(qū)間上取一個(gè)點(diǎn)). 當(dāng)當(dāng)n較大時(shí)較大時(shí), 有有.)( 1)()(21dxxfabffba 63因此因此, 若事先給定誤差限若事先給定誤差限 , 則當(dāng)則當(dāng).3|2 nnTT時(shí)時(shí), 就可停止計(jì)算就可停止計(jì)算, 并認(rèn)為并認(rèn)為 T2n是滿足精度要求的近是滿足精度要求的近似值似值.;412 nnTITI);()(21 ff ).(3122nnnTTTI 64 復(fù)化復(fù)化Simpson公式的逐次分半算法公式的逐次分半算法將區(qū)間將區(qū)間a, b分成分成 n=2m 等分等分, 記記,2mmabh evenodd2)(2)(4)()(3kmkmmkhafkhafbfafhSm, 2 , 1 m

22、稱稱 為為Simpson序列序列.2mS65;1612 nnSISI);,(),(18011)4(4bafhabSIn );,(),(218022)4(42bafhabSIn );()(2)4(1)4( ff ).(15122nnnSSSI 因此因此, 若事先給定誤差限若事先給定誤差限 , 則當(dāng)則當(dāng).15|2 nnSS時(shí)可停止計(jì)算時(shí)可停止計(jì)算, 取取 S2n為滿足精度要求的近似值為滿足精度要求的近似值. 復(fù)化復(fù)化Simpson公式余項(xiàng)的后驗(yàn)估計(jì)公式余項(xiàng)的后驗(yàn)估計(jì)66nnnnnTTTTTI3134)(31222 4 Romberg求積公式求積公式啟示啟示: 是否用是否用 復(fù)化梯形公式余項(xiàng)的后驗(yàn)估

23、計(jì)表明復(fù)化梯形公式余項(xiàng)的后驗(yàn)估計(jì)表明nnTT31342 逼近逼近 I ( f ) 比用比用 T2n要好要好. 事實(shí)上有事實(shí)上有.313422nnnSTT 即梯形值序列的巧妙線性組合得到即梯形值序列的巧妙線性組合得到Simpson序列序列!67以以n=4為例加以說(shuō)明為例加以說(shuō)明. 記記 fk= f (xk), )(2167654321808fffffffffabT )(28642804fffffabT 8abkaxk )(422847654321808fffffffffabT .3)4(488TTS )( 2)( 48642753180fffffffffab 484TT 68).(15122nn

24、nSSSI .1516)(151222nnnnnSSSSSI 逼近逼近 I ( f ) 比用比用 S2n要好要好.回答回答: 是的是的, 記記,15)16(22nnnSSC 15)16(2nnSS 則它恰為復(fù)化則它恰為復(fù)化Cotes公式公式; 且有如下誤差估計(jì)式且有如下誤差估計(jì)式).(| )(|62hOfICn 復(fù)化復(fù)化Simpson公式余項(xiàng)的后驗(yàn)估計(jì)表明公式余項(xiàng)的后驗(yàn)估計(jì)表明 問(wèn)題問(wèn)題: 是否用是否用69.631636422nnnRCC 類似地可以得到類似地可以得到2nR其中其中被稱為被稱為Romberg序列序列.).(| )(|82hOfIRn 截?cái)嗾`差截?cái)嗾`差:701T2T4T8T16T

25、32T2S4S8S16S32S4C8C16C32C8R16R32R3422nnnTTS 151622nnnSSC 636422nnnCCR 停機(jī)準(zhǔn)則停機(jī)準(zhǔn)則:梯形值序列梯形值序列Simpson序列序列Cotes序列序列Romberg序列序列 Romberg求積公式求積公式71例例 計(jì)算計(jì)算.sin10dxxxI 2)1()0(1ffT =0.9207355)5 . 0(212112fTT =0.9397933 )43()41(412124ffTT=0.9445135=0.9456909 )87()85()83()81(812148ffffTT解解先求梯形值序列先求梯形值序列72nnTnSnCn

26、R24180.92073550.93979330.94451350.94569090.94614590.95608690.94608330.94608300.94608310.94608313422nnnTTS 151622nnnSSC .636422nnnCCR 利用只有兩三位有效數(shù)字利用只有兩三位有效數(shù)字的的T1, ,T8 經(jīng)過(guò)三次外推得經(jīng)過(guò)三次外推得到到7位有效數(shù)字位有效數(shù)字. 可見(jiàn)加速的可見(jiàn)加速的效果十分顯著效果十分顯著.用用Romberg算法計(jì)算如下算法計(jì)算如下73 理論依據(jù)理論依據(jù): 復(fù)化梯形公式的余項(xiàng)展開(kāi)復(fù)化梯形公式的余項(xiàng)展開(kāi).記記),(hTTn 定理定理 設(shè)設(shè),)(baCxf 則則 kkhhhIhT24221)( 其中系數(shù)其中系數(shù) k ( k=0, 1, )是與是與 h 無(wú)關(guān)的常數(shù)無(wú)關(guān)的常數(shù). T (h) 逼近逼近 I 的速度是的速度是 O ( h2 )階階.74,3)()2(4)(262411 kkhhhIhThThT 當(dāng)區(qū)間當(dāng)區(qū)間a, b 2n等分時(shí)等分時(shí), 則有則有),2(2hTTn 在定理中以在定理中以 h/2 代替代替 h 得得,2164224221 k