插值法和數(shù)據(jù)擬合 (1)

1、 2 2 3 3問(wèn)題的提出問(wèn)題的提出 010011( )1,( ),(0,1, ),( )( ),(0,1, ).( )1,( )( )( )niiiinnyf xnx xxyf xinyP xP xyinyf xnxyx yxyxa bR xf xP x插值問(wèn)題的數(shù)學(xué)提法：已知函數(shù)在個(gè)點(diǎn)上的函數(shù)值求一個(gè)多項(xiàng)式，使其滿(mǎn)足即要求該多項(xiàng)式的函數(shù)曲線(xiàn)要經(jīng)過(guò)上已知的個(gè)點(diǎn)同時(shí)在其他上要估計(jì)誤差 4 4( )( )( )Rxf xP x 同時(shí)在其它同時(shí)在其它上要估計(jì)誤差上要估計(jì)誤差。當(dāng)當(dāng)插值問(wèn)題插值問(wèn)題n 1 ( )( )( )Rxf xP x n 1 個(gè)點(diǎn)個(gè)點(diǎn)同時(shí)在其它同時(shí)在其它上要估計(jì)誤差上要估計(jì)誤差

2、。當(dāng)當(dāng) 5 5一次插值一次插值 100111(),nPxxyxy當(dāng)時(shí) ， 求 一 次 多 項(xiàng) 式要 求 通 過(guò)兩 點(diǎn) 6 6 20011222( ),nP xxyx yxy當(dāng)時(shí)，求二次多項(xiàng)式要求通過(guò)三點(diǎn) 7 7 設(shè)函數(shù)設(shè)函數(shù) y=f(x)在區(qū)間在區(qū)間a,b連續(xù)連續(xù), , 給定給定n+1個(gè)點(diǎn)個(gè)點(diǎn)已知已知 f(xk)=yk (k=0,1,n) ,在函數(shù)類(lèi)中尋找一函數(shù)在函數(shù)類(lèi)中尋找一函數(shù)g(x) 作為作為 f(x) 的近似表達(dá)式的近似表達(dá)式, ,使?jié)M足使?jié)M足這時(shí)稱(chēng)這時(shí)稱(chēng) y=f(x)為為被插值函數(shù)被插值函數(shù), , g(x) 稱(chēng)為稱(chēng)為插值函數(shù)插值函數(shù), xk 稱(chēng)稱(chēng)為為插值節(jié)點(diǎn)插值節(jié)點(diǎn), ,式式(4.2

3、)(4.2)稱(chēng)為稱(chēng)為插值條件插值條件, ,尋求插值函數(shù)尋求插值函數(shù)g(x) 的方法稱(chēng)為的方法稱(chēng)為插值方法插值方法. .一、一、插值問(wèn)題插值問(wèn)題a x0 x1 xnb (4.1)g(xk)=f(xk)=yk , k=0,1, ,n (4.2) 8 8 nnnnnnnnnnyxaxaxaayxaxaxaayxaxaxaa22101121211000202010(4.4) 設(shè)設(shè) pn(x)=a0+ a1x+an xn ，當(dāng)滿(mǎn)足如下的插值條件當(dāng)滿(mǎn)足如下的插值條件時(shí)，即時(shí)，即得到關(guān)于系數(shù)得到關(guān)于系數(shù) a0、 a1、 、an 的線(xiàn)性方程組的線(xiàn)性方程組pn(xk) = f(xk) = yk , k = 0,

4、1, ,n (4.3) 9 9其系數(shù)行列式為其系數(shù)行列式為Vandermonde(范德蒙范德蒙)行列式行列式 njiijnnnnnnnxxxxxxxxxxxxxxV021211020010)(111),(只要只要 ij 就有就有xi xj , ,因而因而 V(x0,x1,xn) 0 ,于是方程組于是方程組(4.4)(4.4)有唯一解有唯一解, ,也就是說(shuō)也就是說(shuō), ,當(dāng)節(jié)點(diǎn)不重合時(shí)當(dāng)節(jié)點(diǎn)不重合時(shí), , n+1+1個(gè)節(jié)點(diǎn)的個(gè)節(jié)點(diǎn)的插值條件能唯一確定一個(gè)插值條件能唯一確定一個(gè)n次插值多項(xiàng)式次插值多項(xiàng)式, ,從而有從而有: : 定理定理4.1 給定給定n+1 個(gè)互異節(jié)點(diǎn)個(gè)互異節(jié)點(diǎn) x0 , x1 ,

5、 , xn 上的函數(shù)上的函數(shù)值值y0 , y1 , , yn ,則滿(mǎn)足插值條件則滿(mǎn)足插值條件 pn(xk )=f(xk )= yk , k=0, 1, ,n的的n 次插值多項(xiàng)式次插值多項(xiàng)式 pn(x) 存在且唯一。存在且唯一。 1010這樣只要求解方程組這樣只要求解方程組1. 拉格朗日插值多項(xiàng)式拉格朗日插值多項(xiàng)式2. 牛頓插值多項(xiàng)式牛頓插值多項(xiàng)式 3. 分段線(xiàn)性插值多項(xiàng)式分段線(xiàn)性插值多項(xiàng)式 4. 三次樣條插值多項(xiàng)式三次樣條插值多項(xiàng)式便可確定插值多項(xiàng)式便可確定插值多項(xiàng)式 pn(x) ，但相對(duì)來(lái)講計(jì)算較復(fù)雜。，但相對(duì)來(lái)講計(jì)算較復(fù)雜。而插值多項(xiàng)式的唯一性保證了無(wú)論用什么方法獲得滿(mǎn)足而插值多項(xiàng)式的唯一

6、性保證了無(wú)論用什么方法獲得滿(mǎn)足插值條件的多項(xiàng)式都是同一個(gè)多項(xiàng)式插值條件的多項(xiàng)式都是同一個(gè)多項(xiàng)式 pn(x) ，因此可以，因此可以采用其它更簡(jiǎn)便的方法來(lái)確定多項(xiàng)式。下面就介紹幾種采用其它更簡(jiǎn)便的方法來(lái)確定多項(xiàng)式。下面就介紹幾種常用的方法：常用的方法： nnnnnnnnnnyxaxaxaayxaxaxaayxaxaxaa22101121211000202010(4.4) 1111的函數(shù)值的函數(shù)值 nyyy,10已知已知 y=f(x) 在在n+1 個(gè)點(diǎn)個(gè)點(diǎn)bxxxan 10nkyxfxpkkkn, 1 , 0,)()( 構(gòu)造構(gòu)造n次多項(xiàng)式次多項(xiàng)式 pn(x) ,使得使得從而得到從而得到 f(x) 的

7、近似計(jì)算式的近似計(jì)算式 ,),()(baxxpxfn 1212求解求解 p1(x)=a1 x+a0ixiy0 x1x0y1y已知已知使得使得 f(x) p1(x), x x0 , x1.10100101yxxxxyxxxx 根據(jù)點(diǎn)斜式得到根據(jù)點(diǎn)斜式得到如果令如果令1010)(xxxxxl 0101)(xxxxxl 則稱(chēng)則稱(chēng) l0(x) , l1(x)為一次插值多項(xiàng)式的為一次插值多項(xiàng)式的基函數(shù)基函數(shù)。這時(shí)。這時(shí)：并稱(chēng)其為一次并稱(chēng)其為一次Lagrange插值多項(xiàng)式。插值多項(xiàng)式。 0)(, 1)(1000 xlxl1)(, 0)(1101 xlxlf(x) p1(x)=y0l0(x)+y1 l1(x

8、)1010010( )()yyP xyxxxx 1313例例 已知已知 的兩個(gè)采樣點(diǎn)（的兩個(gè)采樣點(diǎn)（100,10）和（）和（121,11）,求求 。115xy x 100 121y=f(x) 10 11)()(0010101xxxxyyyxp 101001011)(yxxxxyxxxxxp 1115 121115 100( )1151011 10.714285100 121121 100p xg10 1 , lg20 1.3010lg12( )lg( )lg(10)1(20)1.3010102011.3010201101( )(20)( )(10)10 201020

9、1010f xx f xx ffxxyyxxlxxl xx例：已知，利用插值一次多項(xiàng)式求的近似值。解：，設(shè)，則插值基本多項(xiàng)式為：， 151510 01 1111.3010( )( )( )(20)(10)101011.3010(12)(12 20)(12 10) 1.06021010lg12 lg10 lg20 lg12 1.0792).P xy l xyl xxxP于是，拉格朗日型一次插值多項(xiàng)式為：故即由和兩個(gè)值的線(xiàn)性插值得到，且具有兩位有效數(shù)字（精確解 1616求解求解 p2(x)=a2x2+a1 x+a0使得使得 f(x) p2(x), x x0 , x2.關(guān)于二次多項(xiàng)式的構(gòu)造采用如下方

10、法：令關(guān)于二次多項(xiàng)式的構(gòu)造采用如下方法：令已知已知ixiy0 x1x0y1y2x2y并由插值條件并由插值條件得到得到)(20100 xxxxyA )(21011xxxxyB )(12022xxxxyC p2(x)=A(x-x1)(x-x2)+B(x-x0)(x-x2)+C(x-x0)(x-x1)p2(x0)=y0 , p2(x1)=y1 , p2(x2)=y2 1717于是得到于是得到則有則有 jijixlijij, 0, 1)( f(x) p2(x)=y0l0(x)+y1 l1(x)+ y2 l2(x)()()(2010210 xxxxxxxxxl 如果令如果令)()()(2101201xx

11、xxxxxxxl )()()(1202102xxxxxxxxxl 并稱(chēng)其為二次并稱(chēng)其為二次Lagrange 插值多項(xiàng)式。插值多項(xiàng)式。 則稱(chēng)則稱(chēng) l0(x) , l1(x)，l(x)為二次插值多項(xiàng)式的為二次插值多項(xiàng)式的基函數(shù)基函數(shù)。這時(shí)。這時(shí)：0201122012010210122021()()()()()()( )()()()()()()x xx xx xx xx x x xp xyyyxx xxxxxxxxxx 1818例例已知已知 , , ,10100 11121 115y解解: : 這里這里x0 0=100=100，y0 0=10=10，x1 1=121=121，y1 1=11, =11

12、, x2 2=144=144，y2 2=12=12，利用拋物線(xiàn)插值公式，利用拋物線(xiàn)插值公式 利用拋物線(xiàn)插值求利用拋物線(xiàn)插值求144122115(121)(144)(100)(144)(115)1011(100 121)(100 144)(121 100)(121 144)(100)(121)1210.72275(144 100)(144 121)xxxxxLxx 19192例 ：已知ixlgiiyxixlgiiyx10152011.17611.3010lg12利用此三值的二次插值多項(xiàng)式求的近似值012012101520(15)(20)1( )1520(10 15)(1020)50(10)(20

13、)1( )1020(15 10)(1520)25(10)(15)1( )1015(20 10)(20 15)50 xxxxxlxxxxxl xxxxxlxxx 解：設(shè)，則 2020ixlgiiyx20 01 12 221( )( )( )( )2015501.17611.301010201015255011.1761(12)1220 12 1512 10 122050251.301012 10 12 151.076650lg12lg121.07923P xy lxy l xy lxxxxxxxP故所以：利用三個(gè)點(diǎn)進(jìn)行拋物線(xiàn)插值得到的的值，與精確值相比，具有 位有效數(shù)字。 2121已知已知n+1

14、組離散數(shù)據(jù)組離散數(shù)據(jù)ixiy0 x1x0y1ynxny按照二次按照二次Lagrange插值多項(xiàng)式的構(gòu)造方法，令：插值多項(xiàng)式的構(gòu)造方法，令：將插值條件將插值條件 pn( x0 )= y0 代入，得到：代入，得到：)()(0201000nxxxxxxyA 同理，由插值條件同理，由插值條件 pn( x1 )= y1 ，得到：得到：)()(1210111nxxxxxxyA 012102011( )()()()()()()()()()nnnnnpxA xxxxxxA xxxxxxA xxxxxx 2222由插值條件由插值條件 pn( xn )= yn ,得得)()(110 nnnnnnxxxxxxyA將

15、將)()(0201000nxxxxxxyA )()(1210111nxxxxxxyA )()(,110 nnnnnnxxxxxxyA代入下式：代入下式：得到：得到：012102011( )()()()()()()()()()nnnnnp xA xxxxxxA xxxxxxA xxxxxx1202010102010121011011()()()()()()( )()()()()()()()()()()()()nnnnnnnnnnnxxxxxxxxxxxxp xyyxxxxxxxxxxxxxxxxxxyxxxxxx 2323也就是也就是011011()()()()()()()()jjnjjjjjj

16、jnxxxxxxxxyxxxxxxxx1202010102010121()()()()()()( )()()()()()()nnnnnx x x xx xx xx xx xp xyyxx xxxxxxxxxx011011()()()()()()nnnnnnxxxxxxyxxxxxx 0110011()()()()()()()()njjnjjjjjjjjnxxxxxxxxyxxxxxxxx 從而就得到了在區(qū)間從而就得到了在區(qū)間 a, ,b 上的關(guān)于函數(shù)上的關(guān)于函數(shù) f(x) 的的n次多項(xiàng)式近似計(jì)算式：次多項(xiàng)式近似計(jì)算式： ( )( ), , nfxpxxa b 242401234012340(2

17、,0)(4,3)(6,5)(8,4)(10,1)24681003541(4)(6)(8)(10)( )(24)(26)(28)(2 10)1(4)(6)(8)(10)384xxxxxyyyyyxxxxlxxxxx例3：求過(guò)點(diǎn)的拉格朗日型插值多項(xiàng)式。解：用4次插值多項(xiàng)式對(duì)5個(gè)點(diǎn)插值：，； 2525123(2 ) (6 ) (8 ) (1 0 )()( 42 ) ( 46 ) ( 48 ) ( 41 0 )1(2 ) (6 ) (8 ) (1 0 )9 6(2 ) (4 ) (8 ) (1 0 )()( 62 ) ( 64 ) ( 68 ) ( 61 0 )1(2 ) (4 ) (8 ) (1 0

18、 )6 4(2 ) (4 ) (6 ) (1 0 )()( 82 ) ( 84 ) ( 86 ) ( 81 0 )1xxxxlxxxxxxxxxlxxxxxxxxxlx，4(2 ) (4 ) (6 ) (1 0 )9 6(2 ) (4 ) (6 ) (8 )()(1 02 ) (1 04 ) (1 06 ) (1 08 )1(2 ) (4 ) (6 ) (8 )3 8 4xxxxxxxxlxxxxx，； 262640 01 12 23 34 4( )( )( )( )( )( )3(2)(6)(8)(10)965(2)(4)(8)(10)644(2)(4)(6)(10)961(2)(4)(6)

19、(8)384P xy lxy l xy lxy l xy lxxxxxxxxxxxxxxxxx 所以 2727拉格朗日插值多項(xiàng)式的截?cái)嗾`差拉格朗日插值多項(xiàng)式的截?cái)嗾`差 ,( )( ),( )( )( )( )( )( )( )0nnnniniinia bP xf xR xR xf xP xxxR xf xP x我們?cè)谏嫌枚囗?xiàng)式來(lái)近似替代函數(shù)其截?cái)嗾`差記作，。當(dāng) 在插值結(jié)點(diǎn) 上時(shí)，估計(jì)截?cái)嗾`差。 2828( )(1)01(1)011( )( ) , ( )( , )( ) , 1( )( )( )(1)!( , ),( )()()()nnnnnnnnnyf xnyfxa byfxa baxxxb

20、P xnxa bR xfxna bxxxxxxx 定理 ：設(shè)函數(shù)的 階導(dǎo)數(shù)在上連續(xù)，在上存在；插值結(jié)點(diǎn)為，是次拉格朗日插值多項(xiàng)式；則對(duì)任意有：其中 2929對(duì)于線(xiàn)性插值，其余項(xiàng)為對(duì)于二次插值，其余項(xiàng)為 303012212lg121lg10lg20lg12P 121.0602,lg121.0792e1.0792 1.06020.0190( )lg1( ),10,20ln101|( )|0.043ln10 101|( )(12 10)(1220)| 8 0.0430.3442f xxfxxxff 分析例，例 中計(jì)算的截?cái)嗾`差在例中，用和計(jì)算，估計(jì)誤差：當(dāng)時(shí) 3131243242lg10 lg15l

21、g20lg12.P 121.0766,e1.0792 1.07660.00262( )8.686 10ln101|(12)| |( )(12 10)(12 15)(1220)|3!18.686 102 3 80.006956fxxRf 在例 中，用，和計(jì)算估計(jì)誤差：故 3232作業(yè)1.已知函數(shù)表求拋物插值多項(xiàng)式及其余項(xiàng)。2已知插值基函數(shù)，證明：當(dāng)m7時(shí)時(shí),舍入誤差亦會(huì)增大舍入誤差亦會(huì)增大.)()!()()()()()(xwnfxLxfxRnnnn1 考考慮慮 6464,01n+1niiiiii+1i : a,bn+1a= x xxx =by = f(x )(i=1,2,.,n),x ,1 j

22、x = y ,(i=1,2,.,n)(2)x ,x,(i=0,1,2,.,n-1)xxx定義設(shè)在區(qū)間上給定個(gè)節(jié)點(diǎn)及節(jié)點(diǎn)上的函數(shù)值作一個(gè)插值函數(shù)（ ） 使其滿(mǎn)足（） （ ）在每個(gè)小區(qū)間上，（ ）都是線(xiàn)性插值函數(shù)。稱(chēng)函數(shù)（ ）為a,b上關(guān)于數(shù)據(jù)（)(0,1,2,., )iyin的分段線(xiàn)性插值函數(shù)。分段分段線(xiàn)性線(xiàn)性插值插值 6565由Lagrange線(xiàn)性插值公式容易寫(xiě)出 的分段表達(dá)式也可以通過(guò)構(gòu)造基函數(shù)的方法來(lái)求 。首先構(gòu)造一組基函數(shù) ,每個(gè) 都滿(mǎn)足2） 在每個(gè)小區(qū)間 上都是線(xiàn)性函數(shù)。這組函數(shù)稱(chēng)為分段線(xiàn)性插值函數(shù)。其表達(dá)式也可以寫(xiě)為) x(111)=(0,1,2,., )i+1iiiiiii+1ii

23、x- xxxxyyxxxinx - xxx() x( )(0,1,2,., )il x in( )il x01) ()( ,0,1,2,., )1ijjil xi jnji()ilx1 ,(0,1,2,., )iix xin10100101( ),0,xxxxlxxx xxx x 6666111111011,( ) ,(1,2,.,1)0,)(,iiiiiiiiiiiiinxxxxxxxxxl xxx xinxxxx xxx11101,( )0,nnnnnnnxxxxxxxlxxx x 6767類(lèi)似于lagarange插值多項(xiàng)式的構(gòu)造，函數(shù)就是所求的分段插值函數(shù)。例已知函數(shù)y=f(x)=1/(

24、1+x2),在區(qū)間0,5上取等距節(jié)點(diǎn)xi=0+i(i=0,1,5)。求分段線(xiàn)性插值函數(shù)，并由此計(jì)算f(4.5)的近似值。解 節(jié)點(diǎn)處的函數(shù)值見(jiàn)下表0()()niiixy lxxi012345yi1.000000.500000.200000.100000.058820.03846 6868解解 分段插值基函數(shù)為分段插值基函數(shù)為5 , 1 (0 1 , 01)(0 xxxxl 1, 15 , 00 1,) 1(, 11)(iixiixixiixixxli)4 , 005 , 44)(xxxxln 6969所以，分段插值函數(shù)為所以，分段插值函數(shù)為)(10000. 0)(20000. 0)(50000.

25、 0)()(3210 xlxlxlxlx)(03846. 0)(05882. 054xlxl)5 . 4(03846. 0)5 . 4(05882. 0)5 . 4()5 . 4(54llf5 . 003846. 05 . 005882. 004864. 004706. 0)5 . 4(f比較，結(jié)果是比較精確的。比較，結(jié)果是比較精確的。 7070分段線(xiàn)性函數(shù)的優(yōu)缺點(diǎn)：1）簡(jiǎn)便易行，而且當(dāng)節(jié)點(diǎn)加密時(shí)，分段線(xiàn)性插值誤差變小，收斂性可以得到保證。2）由于每個(gè)小區(qū)間的插值函數(shù)只依賴(lài)于本段的節(jié)點(diǎn)值，因而每個(gè)節(jié)點(diǎn)只影響相鄰的一兩個(gè)小區(qū)間，不會(huì)擴(kuò)大誤差，從而保證插值結(jié)果的穩(wěn)定性。3）由于分段插值函數(shù)僅在a,

26、b上連續(xù)。一般的，在節(jié)點(diǎn)處插值函數(shù)不可微，因此不滿(mǎn)足工程上光滑度的問(wèn)題。 7171例 對(duì)sinx進(jìn)行線(xiàn)性插值示例解在編輯窗口輸入clear;x=0:10;y=sin(x);xi=0:0.25:10;yi=interp1(x,y,xi);plot(x,y,o,xi,yi) 7272 4.5 Hermite插值法插值法Lagrange插值雖然構(gòu)造比較簡(jiǎn)單，但插值曲線(xiàn)只是在節(jié)點(diǎn)處與原函數(shù)吻合，若還要求在節(jié)點(diǎn)處兩者相切，即導(dǎo)數(shù)值相等，使之與被插函數(shù)的”密切”程度更好，這就要用到帶導(dǎo)數(shù)的插值.0101( ),nnf xax xxbfff設(shè)在節(jié)點(diǎn)處的函數(shù)值為值函數(shù)上的具有一階導(dǎo)數(shù)的插的在區(qū)間為設(shè),)()(

27、baxfxP處必須滿(mǎn)足在節(jié)點(diǎn)即nxxxxP,)(10()()iiiP xf xf()()iiiP xfxf)(,)()1(一階光滑度上具有一階導(dǎo)數(shù)在若要求baxP-(1)0,1,in 7373處必須滿(mǎn)足在節(jié)點(diǎn)即nxxxxP,)(10()()iiiP xf xf()()iiiP xfxf)(,)(,)2(階光滑度階導(dǎo)數(shù)上具有在若要求同樣mmbaxP個(gè)待定的系數(shù)可以解出22n次的多項(xiàng)式可以是最高次數(shù)為因此12)(nxP次多項(xiàng)式作為插值函數(shù)兩個(gè)節(jié)點(diǎn)就可以用3112()()iiiP xfxf()()()()()mmmiiiPxfxfni, 1 ,0-(2)22n 共個(gè)方程 7474定義1. 稱(chēng)滿(mǎn)足(1

28、)或(2)式的插值問(wèn)題為Hermite插值，稱(chēng)滿(mǎn)足(1)或(2)式的插值多項(xiàng)式P(x)為Hermite插值多項(xiàng)式，記為 , 為多項(xiàng)式次數(shù).,( )(),kkHermiteHxkRungek一般次插值多項(xiàng)式的次數(shù) 如果太高會(huì)影響收斂性和穩(wěn)定性現(xiàn)象因此 不宜太大.一、兩點(diǎn)三次Hermite插值先考慮只有兩個(gè)節(jié)點(diǎn)的插值問(wèn)題0101( ),f xxxff設(shè)在節(jié)點(diǎn)處的函數(shù)值為0101,xxff在節(jié)點(diǎn)處的的一階導(dǎo)數(shù)值為作為插值函數(shù)多項(xiàng)式次兩個(gè)節(jié)點(diǎn)最高可以用)(33xHHermite( )kHxk 7575應(yīng)滿(mǎn)足插值條件)(3xH300()Hxf311()Hxf300()Hxf311()Hxf示應(yīng)用四個(gè)插值

29、基函數(shù)表)(3xH3 ,2 ,1 ,0),()(3ixhxHi的插值基函數(shù)為設(shè))()()()()(332211003xhaxhaxhaxhaxH希望插值系數(shù)與Lagrange插值一樣簡(jiǎn)單重新假設(shè)300110011( )( )( )( )( )Hxfxfxfxfx 7676其中1)(00 x0)(00 x1)(00 x300110011( )( )( )( )( )Hxfxfxfxfx 0)(10 x0)(01x1)(11x0)(10 x0)(01x0)(11x0)(00 x0)(10 x0)(01 x0)(11 x0)(10 x0)(01 x1)(11 x可知即可假設(shè)的二重零點(diǎn)是,)(01xx

30、)()()(210baxxxx1)(00 x0)(00 x由 7777可得310)(2xxa3100210)(2)(1xxxxxb)()()(210baxxxx21)(xx 310)(2xxx3100210)(2)(1xxxxx21021)()(xxxx102xxx10021xxx01021xxxx2101xxxx)()(21(201xlxlLagrange插值基函數(shù) 7878)(1x)()(21(210 xlxl類(lèi)似可得)(0 x)()(200 xlxx)(1x)()(211xlxx10121xxxx2010 xxxx0 xx 2101xxxx2010 xxxx1xx )(0 x01021x

31、xxx2101xxxx)()(21(201xlxl即將以上結(jié)果代入 7979300110011( )( )( )( )( )Hxfxfxfxfx得兩個(gè)節(jié)點(diǎn)的三次Hermite插值公式300110011( )( )( )( )( )Hxfxfxfxfx2101(12 ( )( )flxlx2000()( )fxxlx2111()( )f xxlx2010(12 ( )( )fl xlx110112xxfxx2010 xxxx00fxx2101xxxx2010 xxxx11fxx001012xxfxx2101xxxx 8080二、兩點(diǎn)三次Hermite插值的余項(xiàng)兩點(diǎn)三次Hermite插值的誤差為)

32、()()(33xHxfxR0)()()(33iiixHxfxR0)()()(33iiixHxfxR1 , 0i因此可設(shè)的二重零點(diǎn)均為,)(,310 xRxx21203)()()(xxxxxKxR待定其中)(xK 818121203)()()()()(xtxtxKtHtft構(gòu)造輔助函數(shù)0)()()()()(21203xxxxxKxHxfxiiiii1 , 0i0)()()()()(21203xxxxxKxHxfx均是二重零點(diǎn)個(gè)零點(diǎn)至少有因此5)(t連續(xù)使用4次Rolle定理，可得，,10 xx至少存在一點(diǎn)使得0)()4( 82820)(! 4)()()4()4(xKf即! 4)()()4(fxK

33、所以,兩點(diǎn)三次Hermite插值的余項(xiàng)為2120)4(3)()(! 4)()(xxxxfxR10 xx以上分析都能成立嗎？上述余項(xiàng)公式成立上存在時(shí)在當(dāng),)(10)4(xxxf 838311( )( ),0,1,( )( )nniiiiiiHfinHfxxxxyyHermite滿(mǎn)足插值條件:Hermite值問(wèn)題的解插值法定理1 存在且唯一.(22)221210( )( )( )() , , (2(2 !)nnnnjja bf xHfxnRxxx Hermite插值多項(xiàng)式的余項(xiàng)為Hermite插2:值法定理 8484例.1)2(,0)1(21)(3)2(,2)1(21)(ffxfffxf處的導(dǎo)數(shù)值

34、為，在節(jié)點(diǎn)處的函數(shù)值為，在節(jié)點(diǎn)已知.7 . 1 , 5 . 1)(,)(處的函數(shù)值在及的兩點(diǎn)三次插值多項(xiàng)式求xxfxf解:2, 110 xx012,3ff010,1ff 300110011( )( )( )( )( )Hxfxfxfxfx110112xxfxx2010 xxxx00fxx2101xxxx2010 xxxx11fxx001012xxfxx2101xxxx 8585)2(213x21x21x2 x)1(212x22x)(3xH91713323xxx)5 . 1(f)5 . 1(3H625. 2)7 . 1(f)7 . 1(3H931. 2作為多項(xiàng)式插值作為多項(xiàng)式插值,三次已是較高的

35、次數(shù)，次數(shù)再高就有三次已是較高的次數(shù)，次數(shù)再高就有可能發(fā)生可能發(fā)生Runge現(xiàn)象現(xiàn)象,因此，對(duì)有因此，對(duì)有n+1節(jié)點(diǎn)的插值問(wèn)題，節(jié)點(diǎn)的插值問(wèn)題，我們可以使用分段兩點(diǎn)三次我們可以使用分段兩點(diǎn)三次Hermite插值插值 8686例 按下表求hermite插值多項(xiàng)式解 將表中數(shù)據(jù)代入hermite插值多項(xiàng)式 xi01yi01yi3922222320110H ()(12)()0(12)()11001011010(0)()3(1)()90110(32)3(1)9(1)10123xxxxxxxxxx xx xxxxxx 8787Heemite插值的matlab實(shí)現(xiàn)function f = Hermite(

36、x,y,y1,x0)syms t;f = 0.0;if(length(x) = length(y) if(length(y) = length(y1) n = length(x); else disp(y和y的導(dǎo)數(shù)的維數(shù)不相等); return; endelse disp(x和y的維數(shù)不相等！ ); return;end for i=1:n h = 1.0; a = 0.0; for j=1:n if( j = i) h = h*(t-x(j)2/(x(i)-x(j)2); a = a + 1/(x(i)-x(j); end end f = f + h*(x(i)-t)*(2*a*y(i)-y

37、1(i)+y(i);endif (nargin=4) f=subs(f,t,x0);end 8888例 hermite插值應(yīng)用實(shí)例，根據(jù)下表計(jì)算hermite插值多項(xiàng)式，并計(jì)算當(dāng)x=1.56時(shí)的y值。解：在matlab編輯窗口輸入x=1.0:0.2:1.8;y=1 1.0954 1.1832 1.2649 1.3416;y1=0.5 0.4564 0.4226 0.3953 0.3727;f=Hermite(x,y,y1)運(yùn)行結(jié)果的x1.01.21.41.61.8y11.4095 1.18321.24691.3416y 0.50000.4564 0.42260.39530.3727 8989f

38、 =390625/576*(t-6/5)2*(t-7/5)2*(t-8/5)2*(t-9/5)2*(-61/3+64/3*t)+390625/36*(t-1)2*(t-7/5)2*(t-8/5)2*(t-9/5)2*(-114418258618928321/10995116277760000+337232823972231/35184372088832*t)+390625/16*(t-1)2*(t-6/5)2*(t-8/5)2*(t-9/5)2*(6660373488918492043/11258999068426240000+7612884810106783/1801439850948198

39、4*t)+390625/36*(t-1)2*(t-6/5)2*(t-7/5)2*(t-9/5)2*(384779664999124623/21990232555520000-2855713758717179/281474976710656*t)+390625/576*(t-1)2*(t-6/5)2*(t-7/5)2*(t-8/5)2*(8968626627620006931/175921860444160000-7762319875242775/281474976710656*t) 9090下面計(jì)算x=1.56時(shí)的值f=Hermite(x,y,y1,1.56)運(yùn)行結(jié)果f = 1.2490 例4

40、.16 p142頁(yè)解t=0.1 0.5 1 1.5 2 2.5 3;y=0.95 0.84 0.86 1.06 1.5 0.72 1.9;y1=1 1.5 2 2.5 3 3.5 4;yy=hermite(t,y,y1,1.8)t1=0.1:0.01:3;yy1=hermite(t,y,y1,t1);plot(t,y,o,t,y1,*,t1,yy1) 9191運(yùn)行結(jié)果yy = 1.3298 92924.6三次樣條插值1.三次差值和三次樣條插值的matlab命令對(duì)于三次差值和三次樣條插值，matlab中都設(shè)有專(zhuān)用命令三次插值命令的調(diào)用格式如下：yk=pchip(x,y,xk)三次樣條插值命令的調(diào)

41、用格式如下：yk=spline(x,y,xk)參數(shù)x,y, xk, yk的意義及要求與interp1中的完全相同，插值效果與interp1中的參數(shù)”method”分別選用pchip和spline等價(jià)。 9393例 在y=f(x)的函數(shù)關(guān)系中，當(dāng)x=-6.2 -5 -3 -1.7 0 1 3 5 6.2 7.5 8.2時(shí)，y=1 0 -1 0 1 0 -1 0 0.5 -1 0。用“線(xiàn)性插值”、“最近插值”分別求出當(dāng)xi=-6.2:0.2:8.2時(shí)所對(duì)應(yīng)的函數(shù)值，并畫(huà)出曲線(xiàn)圖。解：x=-6.2 -5 -3 -1.7 0 1 3 5 6.2 7.5 8.2;y=1 0 -1 0 1 0 -1 0

42、0.5 -1 0 ;xi=-6.2:0.2:8.2;y1=interp1(x,y,xi,linear);y2=interp1(x,y,xi,nearest);plot(x,y,ro,xi,y1,xi,y2,-);legend(,); 9494 9595例 在y=f(x)的函數(shù)關(guān)系中，當(dāng)x=0 2 3 4 5 6時(shí)，y=-4 1 -2 -4 -4 -5 。用分段三次插值和三次樣條插值分別求出當(dāng)xi=-0.3：0.3：6.8時(shí)所對(duì)應(yīng)的函數(shù)值，并畫(huà)出曲線(xiàn)。解x=0 2 3 4 5 6;y=-4 1 -2 -4 -4 -5 ;xi=-0.3:0.3:6.8;y1=interp1(x,y,xi,pchi

43、p) y1=pchip(x,y,xi)y2=interp1(x,y,xi,spline) y2=spline(x,y,xi)plot(x,y,ro,xi,y1,r,xi,y2,-,linewidth,2);legend(樣本點(diǎn),分段三次插值,三次樣條插值);grid; 9696 9797例對(duì)離散的分布在y=exsinx函數(shù)曲線(xiàn)上的數(shù)據(jù)點(diǎn)進(jìn)行樣條插值計(jì)算。解 x=0 2 4 5 8 12 12.8 17.2 19.9 20;y=exp(x).*sin(x);xx=0:0.5:20;yy=spline(x,y,xx)plot(x,y,o,xx,yy) 9898得到的結(jié)果得到的結(jié)果yy共有共有41個(gè)

44、數(shù)據(jù)點(diǎn)個(gè)數(shù)據(jù)點(diǎn) 9999( )( ) , S xf xa b三三次次樣樣則則條條稱(chēng)稱(chēng)為為在在上上的的插插值值函函數(shù)數(shù)q 三次樣條插值函數(shù)三次樣條插值函數(shù)(2)( ),( ),( ) , S xS xSxa b 都都在在區(qū)區(qū)間間上上連連續(xù)續(xù)1(1)( ),iiS xxx 在在每每個(gè)個(gè)小小區(qū)區(qū)間間上上都都是是三三次次多多項(xiàng)項(xiàng)式式( ) , S xa b則則稱(chēng)稱(chēng)為為區(qū)區(qū)間間上上的的三三次次樣樣條條函函數(shù)數(shù)(3)( )S x而而三三次次樣樣條條函函數(shù)數(shù)滿(mǎn)滿(mǎn)足足(),0,1,iif xy in(),0,1,iiS xy in01, , naxxxba b 為為區(qū)區(qū)間間的的一一個(gè)個(gè)分分割割( ) , :S

45、 xa b如如果果函函數(shù)數(shù)在在區(qū)區(qū)間間上上滿(mǎn)滿(mǎn)足足條條件件 100100q 三次樣條插值原理三次樣條插值原理在在n n個(gè)小區(qū)間構(gòu)造個(gè)小區(qū)間構(gòu)造S(xS(x) )，共有，共有n n個(gè)三次多項(xiàng)式，需確定個(gè)三次多項(xiàng)式，需確定4n4n個(gè)參數(shù)個(gè)參數(shù)在所有節(jié)點(diǎn)上在所有節(jié)點(diǎn)上n+1個(gè)方程個(gè)方程在除端點(diǎn)外的節(jié)點(diǎn)上在除端點(diǎn)外的節(jié)點(diǎn)上3(n-1)個(gè)方程個(gè)方程這樣就得到這樣就得到4n - 2個(gè)方程，為保證待定參數(shù)的唯一性，還差兩個(gè)個(gè)方程，為保證待定參數(shù)的唯一性，還差兩個(gè)方程。為此，常用的方法是對(duì)邊界節(jié)點(diǎn)除函數(shù)值外附加要求，方程。為此，常用的方法是對(duì)邊界節(jié)點(diǎn)除函數(shù)值外附加要求，這就是所謂的邊界條件。根據(jù)實(shí)際問(wèn)題的不同

46、，三次樣條插值這就是所謂的邊界條件。根據(jù)實(shí)際問(wèn)題的不同，三次樣條插值常用到下列常用到下列三類(lèi)邊界條件三類(lèi)邊界條件。( ),0,1,iiS xy in1, 1,)()(limniyxSxSiixxi1, 1),()(limnixSxSixxi1, 1,)()(lim niMxSxSiixxi 101101 m邊界條件：邊界條件： 即給定端點(diǎn)處的一階導(dǎo)數(shù)值。即給定端點(diǎn)處的一階導(dǎo)數(shù)值。 M邊界條件：邊界條件： 即給定端點(diǎn)處的二階導(dǎo)數(shù)值。即給定端點(diǎn)處的二階導(dǎo)數(shù)值。 周期邊界條件周期邊界條件： 當(dāng)當(dāng)y = g(x)是以是以 b-a= x0 - xn 為周期的周期函數(shù)時(shí)，要求為周期的周期函數(shù)時(shí)，要求S(x

47、)也是周期函數(shù)，故端點(diǎn)處要滿(mǎn)足也是周期函數(shù)，故端點(diǎn)處要滿(mǎn)足S(a + 0) = S(b - 0), S(a + 0) = S(b - 0)，此條件稱(chēng)為周期條件此條件稱(chēng)為周期條件。00() , ()nnS xyS xy00() , ()nnSxySxy 102102三次樣條插值的matlab實(shí)現(xiàn)實(shí)例例 在區(qū)間的兩個(gè)端點(diǎn)處給出的條件，稱(chēng)為邊界條件。給定兩端點(diǎn)處的二階導(dǎo)數(shù)值S(a)=y0,S(b)=yn,然后根據(jù)下面的數(shù)據(jù)點(diǎn)求出其三次樣條插值多項(xiàng)式，并計(jì)算當(dāng)x=3.5時(shí)y的值。解：在matlab命令窗口輸入x=1:8;y=0.84 0.91 0.14 -0.76 -0.96 -0.28 0.66,0

48、.99f,f0=ThrSample1(x,y,0.54,-0.15,3.5)x12345678y0.840.910.14-0.76-0.96-0.280.660.99y0.5403-0.1455 103103程序運(yùn)行結(jié)果如下f =(7/25*t-7/10)*(t-4)2+(-171/25+38/25*t)*(t-3)2+(-71823/72775*t+215469/72775)*(4-t)2-(-37788/14555+9447/14555*t)*(t-3)2f0 = -0.3522 104104程序運(yùn)行時(shí)調(diào)用的三次樣條函數(shù)如下：function f,f0 = ThrSample1 (x,y,

49、y_1, y_N,x0)syms t;f = 0.0;f0 = 0.0;if(length(x) = length(y) n = length(x);else disp(x和y的維數(shù)不相等！); return;end %維數(shù)檢查for i=1:n if(x(i)=x0) index = i; break; endend %找到x0所在區(qū)間for i=1:n if(x(i)=x0) index = i; break; endend %找到x0所在區(qū)間 A = diag(2*ones(1,n); %求解m的系數(shù)矩陣u = zeros(n-2,1);lamda = zeros(n-1,1);c =

50、zeros(n,1);for i=2:n-1 u(i-1) = (x(i)-x(i-1)/(x(i+1)-x(i-1); lamda(i) = (x(i+1)-x(i)/(x(i+1)-x(i-1); c(i) = 3*lamda(i)*(y(i)-y(i-1)/(x(i)-x(i-1)+ . 3*u(i-1)*(y(i+1)-y(i)/(x(i+1)-x(i); A(i, i+1) = u(i-1); A(i, i-1) = lamda(i); %形成系數(shù)矩陣及向量cendc(1) = 2*y_1;c(n) = 2*y_N;m = followup(A,c); %用追趕法求解方程組h = x

51、(index+1) - x(index); %x0所在區(qū)間長(zhǎng)度f(wàn) = y(index)*(2*(t-x(index)+h)*(t-x(index+1)2/h/h/h + . y(index+1)*(2*(x(index+1)-t)+h)*(t-x(index)2/h/h/h + . m(index)*(t-x(index)*(x(index+1)-t)2/h/h - . m(index+1)*(x(index+1)-t)*(t-x(index)2/h/h; %x0所在區(qū)間的插值函數(shù)f0 = subs(f,t,x0); %x0處的插值 105105其中的followup函數(shù)如下：function

52、 x=followup(A,b)n=rank(A);for i=1:n if A(i,i)=0 disp (error:對(duì)角線(xiàn)有元素0！) return; endendd=ones(n,1);a=ones(n-1,1);c=ones(n-1);for i=1:n-1 a(i,1)=A(i+1,i); c(i,1)=A(i,i+1); d(i,1)=A(i,i);endd(n,1)=A(n,n);for i=2:n d(i,1)=d(i,1)-(a(i-1,1)/d(i-1,1)*c(i-1,1); b(i,1)=b(i,1)-(a(i-1,1)/d(i-1,1)*b(i-1,1);endx(n

53、,1)=b(n,1)/d(n,1);for i=n-1:-1:1 x(i,1)=(b(i,1)-c(i,1)*x(i+1,1)/d(i,1);end 106106例 在區(qū)間的兩個(gè)端點(diǎn)處給出的條件，稱(chēng)為邊界條件。給定兩端點(diǎn)處的二階導(dǎo)數(shù)值S(a)=y0,S(b)=yn,然后根據(jù)下面的數(shù)據(jù)點(diǎn)求出其三次樣條插值多項(xiàng)式，并計(jì)算當(dāng)x=3.5時(shí)y的值。解：x=1:8;y=0.84 0.91 0.14 -0.76 -0.96 -0.28 0.66,0.99f,f0=ThrSample2(x,y,0.84,-0.99,3.5)x12345678y0.840.910.14-0.76-0.96-0.280.660.

54、99y-0.8415-0.9894 107107運(yùn)行結(jié)果如下f =(7/25*t-7/10)*(t-4)2+(-171/25+38/25*t)*(t-3)2+(-23849/23288*t+71547/23288)*(4-t)2-(-5730726096775419/2251799813685248+5730726096775419/9007199254740992*t)*(t-3)2f0 = -0.3585 108108程序運(yùn)行調(diào)用的三次樣條函數(shù)如下function f,f0 = ThrSample2 (x,y,y2_1, y2_N,x0)syms t;f = 0.0;f0 = 0.0;if

55、(length(x) = length(y) n = length(x);else disp(x和y的維數(shù)不相等！); return;end %維數(shù)檢查for i=1:n if(x(i)=x0) index = i; break; endend %找到x0所在區(qū)間for i=1:n if(x(i)=x0) index = i; break; endend %找到x0所在區(qū)間 A = diag(2*ones(1,n); %求解m的系數(shù)矩陣A(1,2) = 1;A(n,n-1) = 1;u = zeros(n-2,1);lamda = zeros(n-1,1);c = zeros(n,1);for

56、 i=2:n-1 u(i-1) = (x(i)-x(i-1)/(x(i+1)-x(i-1); lamda(i) = (x(i+1)-x(i)/(x(i+1)-x(i-1); c(i) = 3*lamda(i)*(y(i)-y(i-1)/(x(i)-x(i-1)+ . 3*u(i-1)*(y(i+1)-y(i)/(x(i+1)-x(i); A(i, i+1) = u(i-1); A(i, i-1) = lamda(i); %形成系數(shù)矩陣及向量cendc(1) = 3*(y(2)-y(1)/(x(2)-x(1)-(x(2)-x(1)*y2_1/2;c(n) = 3*(y(n)-y(n-1)/(x(

57、n)-x(n-1)-(x(n)-x(n-1)*y2_N/2;m = followup(A,c); %用追趕法求解方程組h = x(index+1) - x(index); %x0所在區(qū)間長(zhǎng)度f(wàn) = y(index)*(2*(t-x(index)+h)*(t-x(index+1)2/h/h/h + . y(index+1)*(2*(x(index+1)-t)+h)*(t-x(index)2/h/h/h + . m(index)*(t-x(index)*(x(index+1)-t)2/h/h - . m(index+1)*(x(index+1)-t)*(t-x(index)2/h/h; %x0所在

58、區(qū)間的插值函數(shù)f0 = subs(f,t,x0); %x0處的插值 109109例4-25根據(jù)表中的數(shù)據(jù)，用程序求三次樣條插值函數(shù)解 用程序進(jìn)行計(jì)算x=0,1,2,3;y=0,0.5,2,1.5;dx0=0.2;dxn=-1; s=csfit(x,y,dx0,dxn)回車(chē)得到s = 0.4800 -0.1800 0.2000 0 -1.0400 1.2600 1.2800 0.5000 0.6800 -1.8600 0.6800 2.0000 x0123f(x)00.521.5f(x)0.2-1 1101103213223230.480.180.2 ,0,11.041.261.280.5,1,

59、20.681.860.682,2,3SxxxxSxxxxSxxxx 111111程序調(diào)用時(shí)的函數(shù)如下：function S=csfit(x,y,dx0,dxn)% x,y分別為n個(gè)節(jié)點(diǎn)的橫坐標(biāo)所組成的向量及縱坐標(biāo)組成的向量% dx0,dn分別為S的導(dǎo)數(shù)在x0,xn處的值n=length(x)-1;h=diff(x);d=diff(y)./h;a=h(2:n-1);b=2*(h(1:n-1)+h(2:n);c=h(2:n);u=6*diff(d);b(1)=b(1)-h(1)/2;u(1)=u(1)-3*(d(1)-dx0);b(n-1)=b(n-1)-h(n)/2;u(n-1)=u(n-1)-

60、3*(dxn-d(n);for k=2:n-1 temp=a(k-1)/b(k-1); b(k)=b(k)-temp*c(k-1); u(k)=u(k)-temp*u(k-1);endm(n)=u(n-1)/b(n-1);for k=n-2:-1:1 m(k+1)=(u(k)-c(k)*m(k+2)/b(k);endm(1)=3*(d(1)-dx0)/h(1)-m(2)/2;m(n+1)=3*(dxn-d(n)/h(n)-m(n)/2;for k=0:n-1 S(k+1,1)=(m(k+2)-m(k+1)/(6*h(k+1); S(k+1,2)=m(k+1)/2; S(k+1,3)=d(k+1