第二章 插值法及其matlab實現 ppt課件

1、主講教師：牛曉穎河北大學質監(jiān)學院 描畫事物之間的數量關系：函數。 有兩種情況： 一是表格方式一組離散的數據來表示函數關系；另一種是函數雖然有明顯的表達式，但很復雜，不便于研討和運用。 從實踐需求出發(fā)：對于計算結果允許有一定的誤差，可以把函數關系用一個簡單的便于計算和處置的近似表達式來替代，從而使問題得到簡化。普通地，構造某種簡單函數替代原來函數。0 引言引言第二章第二章 插值插值Interpolation法法插值法就是一種根本方法插值法就是一種根本方法當準確函數當準確函數 y = f(x) 非常復雜或未知時，在一系列節(jié)非常復雜或未知時，在一系列節(jié)點點 x0 xn 處測得函數值處測得函數值 y0

2、 = f(x0), yn = f(xn)，由此構造一個簡單易算的近似函數由此構造一個簡單易算的近似函數 g(x) f(x)，滿足，滿足條件條件g(xi) = f(xi) (i = 0, n)。這里的。這里的 g(x) 稱為稱為f(x) 的插值函數。的插值函數。x0 x1x2x3x4xg(x) f(x) 根據實踐需求，可以用各種不同的函數來近似原來的函數。最常用的插值函數是最常用的插值函數是 ?多項式：多項式： 代數多項式最簡單，計算其值只需用到加、減乘運算，且積分和微分都很方便； 所以常用它來近似表示表格函數(或復雜函數),這樣的插值方法叫做代數插值法，簡稱插值法。1 拉格朗日多項式拉格朗日多

3、項式 niyxPiin,., 0,)(= = =求求 n 次多項式次多項式 使得使得nnnxaxaaxP = =10)(條件：無重合節(jié)點，即條件：無重合節(jié)點，即jixx ji n = 1知知 x0 , x1 ; y0 , y1 ，求，求xaaxP101)( = =使得使得111001)(,)(yxPyxP= = =可見可見 P1(x) 是過是過 ( x0 , y0 ) 和和 ( x1, y1 ) 兩點的直兩點的直線。線。)(1xP101xxxx- - -010 xxxx- - -= y0 + y11.1線性插值線性插值兩點式兩點式)()(0010101xxxxyyyxP- - - - = =點

4、斜式點斜式)(001010 xxxxxxy- - - - = =()ff1.2 二次插值二次插值n = 2知知 x0 , x1 , x2; y0 , y1 ,y2 , 求求22102)(xaxaaxP = =使得使得002,)(yxP112)(yxP= = =222)(yxP= =, 為求為求P2(x),將三點代入其表達式將三點代入其表達式,即可得到三個方程式即可得到三個方程式,從而聯立方程組解出系數從而聯立方程組解出系數a0, a1, a2即可即可:2020100 xaxaay = =2121101xaxaay = =2222102xaxaay = =方程組的解能否存方程組的解能否存在在?

5、假設存在解假設存在解,能能否獨一否獨一?!當當 x0 , x1 , x2互異時互異時,方程組的解存在且獨一方程組的解存在且獨一.注：顯然有注：顯然有, 求求n 次插值時次插值時, 由由n +1個點可有個點可有n +1個方程個方程, 聯立方程組即可求出插值多項式的聯立方程組即可求出插值多項式的n +1個系數個系數. 然而然而,方程組的求解也并不是一件容易的事。方程組的求解也并不是一件容易的事。1.2.1 待定系數法待定系數法 對于線性插值的兩種方式解進展適當的分析對于線性插值的兩種方式解進展適當的分析, , 從中尋從中尋求規(guī)律而得到啟發(fā)求規(guī)律而得到啟發(fā), ,就有了所謂的拉格朗日插值法就有了所謂的

6、拉格朗日插值法( (公式公式) )和牛頓插值和牛頓插值( (公式公式).). 我們先來看看如何得到二次拉格朗日插值公式。 首先, 線性插值的兩點式可看作是兩個特殊的一次式的一種線性組合.101xxxx- - -010 xxxx- - -)(1xP= y0 + y1 = = =10)(iiiyxl兩點式兩點式l0(x)l1(x)本質上本質上0lx和和1lx即是滿足函數表即是滿足函數表 的一次插值多項式的一次插值多項式 ,稱稱l0(x)和和l1(x)為以為以x0,x1為節(jié)點的根本插為節(jié)點的根本插值多項式，也稱為線性插值的插值基函數值多項式，也稱為線性插值的插值基函數 。 于是，線性插值即是用基函數

7、的線性組合來構造的于是，線性插值即是用基函數的線性組合來構造的. 1.2.2 基函數法基函數法稱為拉氏基函稱為拉氏基函數數 ，滿足，滿足 li(xj)=ij 顯然有顯然有l(wèi)0(x)+ l1(x)1.這里，這里， l0(x)和和l1(x)具有如下性質：具有如下性質：l0(x0)=1, l0(x1)=0, l1(x0)=0, l1(x1)=1, , 1, 0, 0, 0, 1, 0, 0, 0, 1221202211101201000=xlxlxlxlxlxlxlxlxl 由此啟發(fā)，我們希望二次插值也能由一些二次插由此啟發(fā)，我們希望二次插值也能由一些二次插值基函數來線性組合值基函數來線性組合:這時

8、，這時，l0(x), l1(x), l2(x)都是二次多項式，且應滿足都是二次多項式，且應滿足滿足滿足(2.1)式的式的 l i(x) 能否存在能否存在?假設存在，具有什么方式呢假設存在，具有什么方式呢?2.1 2211002yxlyxlyxlxp=同理可得同理可得 l1(x) 1(x x0)(x x2), l2(x) 2(x x0)(x x1),1(x1x0)(x1x2)12(x2x0)(x2x1)1此即二次拉格朗日插值公式此即二次拉格朗日插值公式, 其中其中, l0(x), l1(x), l2(x)是滿是滿足足(2.1)的特殊的特殊(根本根本)二次插值多項式二次插值多項式;稱為二次插值基函

9、數稱為二次插值基函數.P2(x)= y0+ y1+ y2(x -x0)(x -x2)(x1-x0)(x1-x2)(x -x1)(x -x2)(x0-x1)(x0-x2)(x -x0)(x -x1)(x2-x0)(x2-x1) 先思索先思索 l0(x)。因。因 l0(x)是以是以 x1, x2 為零點的二次多項為零點的二次多項式式,所以它可寫成所以它可寫成 l0(x) 0(x x1)(x x2), 其中其中0 是是待定系數。待定系數。 又由于又由于 l0( x0)=1，所以，所以0(x0 x1)(x0 x2)1，那么可有，那么可有0(x0 x1)(x0 x2)1 l0(x) 0(x x1)(x

10、x2), n 1希望找到希望找到li(x)，i = 0, , n 使得使得 li(xj)=ij ；然后；然后令令 = = =niiinyxlxP0)()(，那么顯然有，那么顯然有Pn(xi) = yi 。li(x)每個每個 li 有有 n 個根個根 x0 xi xn = =- -= =- - - -= =njj i jiniiixxCxxxxxxCxl00)().().()( - -= = =j i jiiiixxCxl)(11)(= = - - -= =njijjijixxxxxl0)()()( = = =niiinyxlxL0)()( 拉格朗日 多項式與與 有關，而與有關，而與 無關無關節(jié)

11、點節(jié)點f1.3 n 次插值次插值定理定理 (獨一性獨一性) 滿足滿足 的的 n 階插值階插值多項式是獨一存在的。多項式是獨一存在的。niyxPii,., 0,)(= = =證明：證明： ( 存在性可利用存在性可利用Vandermonde 行列式論證行列式論證)反證：假設不獨一，那么除了反證：假設不獨一，那么除了Ln(x) 外還有另一外還有另一 n 階多項式階多項式 Pn(x) 滿足滿足 Pn(xi) = yi 。調查調查 那么那么 Qn 的階的階數數, )()()(xLxPxQnnn- -= = n而而 Qn 有有 個不同的根個不同的根n + 1x0 xn注：假設不將多項式次數限制為注：假設不

12、將多項式次數限制為 n ，那么插值多項式不獨，那么插值多項式不獨一。一。例如例如 也是一個插值也是一個插值多項式，其中多項式，其中 可以是恣意多項式?？梢允琼б舛囗検?。= =- - = =niinxxxpxLxP0)()()()()(xp = = =niiinyxlxL0)()(設節(jié)點設節(jié)點)1( nf在在a , b內存在內存在, 調查截斷誤差調查截斷誤差)()()(xLxfxRnn- -= =, baCfn bxxxan 10，且，且 f 滿足條件滿足條件 ,Rolles Theorem: 假設假設 充分光滑，充分光滑， ，那么，那么存在存在 使得使得 。)(x 0)()(10= = =xx

13、 ),(10 xx 0)(= = 推行：假推行：假設設0)()()(210= = = =xxx ),(),(211100 xxxx 使得使得0)()(10= = = = ),(10 使得使得0)(= = 0)()(0= = = =nxx 存在存在),(ba 使得使得0)()(= = nRn(x) 至少有至少有 個根個根n+1 = =- -= =niinxxxKxR0)()()(恣意固定恣意固定 x xi (i = 0, , n), 調調查查 = =- - -= =niixtxKtRnt0)()()()( (t)有有 n+2 個不同的根個不同的根 x0 xn x),(, 0)()1(baxxn

14、= = !)1()()()1(-nxKRxnn 留意這里是對留意這里是對 t 求導求導= = - - - !)1)()()()1()1(nxKLfxnnxn !)1()()()1( = = nfxKxn = = - - = =niixnnxxnfxR0)1()(! ) 1()()( 1.4 插值余項插值余項 (Remainder) 注：注： 通常不能確定通常不能確定 x , 而是估計而是估計 , x(a,b) 將將 作為誤差估計上限。作為誤差估計上限。1)1()( nnMxf= = - - niinxxnM01|)!1(當當 f(x) 為任一個次數為任一個次數 n 的多項式的多項式時，時， ,

15、 可知可知 ，即插，即插值多項式對于次數值多項式對于次數 n 的多項式是準的多項式是準確的。確的。0)()1( xfn0)( xRn.)(應應用用的的高高階階導導數數存存在在時時才才能能余余項項表表達達式式只只有有在在xf,)()( )(61)(2,)( )(21)()(21)(1202102101021xxxxxxxxfxRnxxxxxxfxfxRn - - - - = = = - - - = = = = = ，時時，拋拋物物插插值值的的余余項項為為當當，時時，線線性性插插值值余余項項為為當當例例1 求經過求經過A(0,1)，B(1,2)，C(2,3)三個插值點的插值多項式三個插值點的插值多

16、項式.解：三個插值節(jié)點及對應的函數值為解：三個插值節(jié)點及對應的函數值為.322110221100= = = = = = =yxyxyx，；，；，13)12)(02()1)(0(2)21)(01()2)(0(1)20)(10()2)(1()()()()()()()(2120210121012002010212 = = - - - - - - - - - - - - - - -= =- - - - - - - - - - - - - - -= =xxxxxxxyxxxxxxxxyxxxxxxxxyxxxxxxxxxL由拋物插值公式得由拋物插值公式得例例2：知：知233sin,214sin,216s

17、in=分別利用分別利用 sin x 的的1次、次、2次次 Lagrange 插值計算插值計算 sin 50 并估計誤差。并估計誤差。 解：解：0 x1x2x185500 =n = 1分別利用分別利用x0, x1 以及以及 x1, x2 計算計算4,610 =xx利用利用216/4/6/214/6/4/)(1 - - - - - -= = xxxL這里這里)3,6(,sin)(,sin)()2( - -= = =xxxfxxf而而)4)(6(!2)()(,23sin21)2(1 - - -= = xxfxRxx00762. 0)185(01319. 01- - - - Rsin 50 = 0.7

18、660444)185(50sin10 L0.77614外推外推 (extrapolation ) 的實踐誤差的實踐誤差 0.010013,421 = = =xx利用利用sin 50 0.76008, 00660. 018500538. 01 R內插內插 (interpolation ) 的實踐誤差的實踐誤差 0.00596內插通常優(yōu)于外推。選擇要計算的內插通常優(yōu)于外推。選擇要計算的 x 所在的區(qū)間的端點，插值效果較好。所在的區(qū)間的端點，插值效果較好。n = 223)()(21)()(21)()()(4363463464363646342 - - - - - - - - - - - - - -

19、-= = xxxxxxxL)185(50sin20 L0.7654323cos21;)3)(4)(6(!3cos)(2 - - - - -= =xxxxxxR 00077. 018500044. 02 Rsin 50 = 0.76604442次插值的實踐誤差次插值的實踐誤差 0.00061高次插值通常優(yōu)于低高次插值通常優(yōu)于低次插值次插值,但絕對不是次但絕對不是次數越高就越好數越高就越好 例例3 思索下述的插值法問題：求二次多項式思索下述的插值法問題：求二次多項式P(x)，滿足，滿足P(x0) = y0, 其中其中 是已給的數據并給出使這一問題的解存在且獨一的條件是已給的數據并給出使這一問題的解

20、存在且獨一的條件.，2211)()(yxPyxP=21020yyyxx、，解：設解：設 那么那么 由知條件有由知條件有，cbxaxxP=2)(.2)(baxxP=11222200202ybaxycbxaxycbxax0012111222020 xxxxx0)()(22220201-xxxxx)0(220201-xxxxx即即 所以所以故原問題的獨一可解性就歸結為上述方程組的獨一可解性而后故原問題的獨一可解性就歸結為上述方程組的獨一可解性而后者獨一可解的充要條件為者獨一可解的充要條件為這就是這就是Px存在且獨一的條件。存在且獨一的條件。 回想回想 拉格朗日插值公式拉格朗日插值公式 .)10(.出

21、有牛頓插值公式出有牛頓插值公式形式，從而導形式，從而導項式變形為便于計算的項式變形為便于計算的這一缺點，可把插值多這一缺點，可把插值多為了克服為了克服的的實際計算中是很不方便實際計算中是很不方便式也要發(fā)生變化，這在式也要發(fā)生變化，這在均要隨之變化，整個公均要隨之變化，整個公，時，全部插值基函數時，全部插值基函數但是，當插值節(jié)點增加但是，當插值節(jié)點增加中非常方便中非常方便結構緊湊，在理論分析結構緊湊，在理論分析式，公式式，公式得到拉格朗日插值多項得到拉格朗日插值多項利用插值基函數很容易利用插值基函數很容易nklk= =1.5 拉格朗日插值公式的優(yōu)缺陷拉格朗日插值公式的優(yōu)缺陷n 1希望找到希望找到

22、li(x)，i = 0, , n 使得使得 li(xj)=ij ；然后；然后令令 = = =niiinyxlxP0)()(，那么顯然有，那么顯然有Pn(xi) = yi 。li(x)每個每個 li 有有 n 個根個根 x0 xi xn = =- -= =- - - -= =njj i jiniiixxCxxxxxxCxl00)().().()( - -= = =j i jiiiixxCxl)(11)(= = - - -= =njijjijixxxxxl0)()()( = = =niiinyxlxL0)()( 拉格朗日 多項式與與 有關，而與有關，而與 無關無關節(jié)點節(jié)點f拉格朗日插值公式拉格朗日

23、插值公式 Lagrange插值公式插值公式(利用插值基函數很容易得到利用插值基函數很容易得到): 含義直觀含義直觀,構造緊湊構造緊湊,在實際分析中非常方便在實際分析中非常方便; 計算機上實現也很容易計算機上實現也很容易. 也有一些缺陷：也有一些缺陷： 一是計算量大，這是顯然的；另外，還有一個更嚴重的一是計算量大，這是顯然的；另外，還有一個更嚴重的缺陷，當插值節(jié)點添加時，全部插值基函數均要隨之變化，缺陷，當插值節(jié)點添加時，全部插值基函數均要隨之變化，整個計算任務必需從頭開場：不僅原來的每一項都要改動，整個計算任務必需從頭開場：不僅原來的每一項都要改動，還要添加一項計算。還要添加一項計算。 為抑制

24、上述兩個缺陷為抑制上述兩個缺陷, 努力：把插值多項式變形為便于計算的方式。努力：把插值多項式變形為便于計算的方式。 希望：計算改動的過程中希望：計算改動的過程中,盡能夠能利用已有的計算結果盡能夠能利用已有的計算結果. 下面我們將看到下面我們將看到,這是能夠的。我們可以有具有這是能夠的。我們可以有具有“承襲性承襲性的所謂牛頓公式。的所謂牛頓公式。*1 拉格朗日插值的拉格朗日插值的Matlab實現實現常用的幾個函數常用的幾個函數(一一) POLY 函數函數調用格式：調用格式：Y = poly (V)(二二) POLYVAL 函數函數調用格式：調用格式：Y = polyval(p,x)前往矩陣特征多

25、項式的系數。前往n次多項式p在x處的值。輸入變量p=p0 p1 p2pn是一個長度為n+1的向量，其元素為按降陳列的多項式系數。如輸入程序如輸入程序Y = poly (3)；那么前往；那么前往1 -3，即其特征多項式為，即其特征多項式為x-3如對多項式如對多項式 ，計算在，計算在x=5,7,9的值。的值。 1232=xxxp(三三) POLY2SYM 函數函數調用格式一：調用格式一：poly2sym (C)調用格式二：調用格式二：f1=poly2sym(C,V) 或或 f2=poly2sym(C, sym (V) )輸入程序 p = 3 2 1; x=5,7,9; polyval(p,x)結果

26、為ans =86 162 262把多項式的系數向量轉化為符號多項式如輸入程序如輸入程序 c=1 -3;y=poly2sym(c) ；那么前往；那么前往y =x-3如輸入程序如輸入程序 c=1 -3;y=poly2sym(c,s) ；或者或者 c=1 -3;y=poly2sym(c,sym(s) )； 那么前往那么前往y =s-3 (四四) CONV 函數函數調用格式：調用格式：C =conv (A, B)前往兩個多項式乘積的多項式系數，即卷積和如輸入程序如輸入程序 A=1 -3;B=1 -2;C=conv(A,B)那么前往那么前往C = 1 -5 6A、B為兩個多項式的系數(五五) DECON

27、V 函數函數調用格式：調用格式：Q,R =deconv (B,A)conv函數的逆函數，前往兩個多項式相除的多項式系數及其他數，即反卷積和假設輸入程序假設輸入程序 C=1 -5 7;A=1 -3;B,R=deconv(C,A)那么輸出那么輸出B = 1 -2 R =0 0 1如輸入程序如輸入程序 C=1 -5 6;A=1 -3;B,R=deconv(C,A)那么輸出那么輸出B = 1 -2 R =0 0 0(六六) roots(poly(1:n)命令命令調用格式：調用格式：roots(poly(1:n) (七七) det(a*eye(size (A) - A)命令命令調用格式：調用格式：b=d

28、et(a*eye(size (A) - A)如輸入程序如輸入程序 C=1 -3 2;roots(C)輸出結果為輸出結果為ans = 2.0000 1.0000數值解數值解輸入輸入 poly(1:2)輸出輸出ans = 1 -3 2輸入輸入 roots(poly(1:2)輸出輸出ans = 2 11到到n的等差數列，的等差數列，其中其中n roots(poly(1:4)輸出輸出ans = 4.0000 3.0000 2.0000 1.0000假設常數a為矩陣A的特征值，那么前往0前往n到1的數值例例*1-1知函數知函數f(x)在在1,3上具有二階延續(xù)導數，上具有二階延續(xù)導數， ，且，且滿足條件滿

29、足條件f(1) =1，f(3)=2，求線性插值多項式和函數值，求線性插值多項式和函數值 f(1.5)，并估計其誤差。并估計其誤差。解：輸入程序解：輸入程序 X=1,3;Y=1,2; l01= poly(X(2)/( X(1)- X(2), l11= poly(X(1)/( X(2)- X(1), l0=poly2sym (l01),l1=poly2sym (l11), P = l01* Y(1)+ l11* Y(2), L=poly2sym (P),x=1.5; Y = polyval(P,x)運轉后輸出基函數運轉后輸出基函數l0和和l1及其插值多項式的系數向量及其插值多項式的系數向量P、插值

30、、插值多項式多項式L和插值和插值Y為為l0 = l1 = L = Y =-1/2*x+3/2 1/2*x-1/2 1/2*x+1/2 1.2500一、線性插值的一、線性插值的Matlab程序程序逗號會將本句程序的結果輸出到命令窗口，分號那么不會輸出輸入程序 M=5;R1=M*abs(x-X(1)* (x-X(2)/2運轉后輸出誤差限為 R1 = 1.8750例例*1-2求函數求函數f(x)=e-x在在0,1上線性插值多項式，并估計其誤差上線性插值多項式，并估計其誤差.解：解： 輸入程序輸入程序 X=0,1; Y =exp(-X) , l01= poly(X(2)/( X(1)- X(2), l

31、11= poly(X(1)/( X(2)-X(1), l0=poly2sym (l01),l1=poly2sym (l11), P = l01* Y(1)+ l11* Y(2), L=poly2sym (P),運轉后輸出基函數運轉后輸出基函數l0和和l1及其插值多項式的系數向量及其插值多項式的系數向量P和插值多項式和插值多項式L為為l0 = -x+1 l1 = x P =-0.6321 1.0000L =-1423408956596761/2251799885248*x+1 輸入程序 M=1;x=0:0.001:1; R1=M*max(abs(x-X(1).*(x-X(2)./2運轉后輸出誤差

32、限為 R1 = 0.1250.點乘和點除表示向量對應相乘和對應相除二、拋物線插值的二、拋物線插值的Matlab程序程序例例*1-3求將區(qū)間求將區(qū)間0,/2分成分成n等份等份n=1,2，用，用y=f(x)=cos(x)產生產生n+1個節(jié)點，分別作線性插值函數和拋物線插值函數，并個節(jié)點，分別作線性插值函數和拋物線插值函數，并用它們分別計算用它們分別計算cos (/6) (取四位有效數字取四位有效數字)，并估計其誤差，并估計其誤差.解：解：n=1輸入程序輸入程序 X=0,pi/2; Y =cos(X) ,l01= poly(X(2)/( X(1)- X(2), l11= poly(X(1)/( X(

33、2)- X(1), l0=poly2sym (l01),l1=poly2sym (l11), P = l01* Y(1)+ l11* Y(2), L=poly2sym (P),x=pi/6; Y = polyval(P,x)運轉后輸出基函數l0和l1及其插值多項式的系數向量P、插值多項式和插值為l0 =-5734161222659/9007199254740992*x+1l1 =5734161222659/9007199254740992*xP = -0.6366 1.0000L =-5734161222659/9007199254740992*x+1Y =0.6667輸入程序 M=1;x=p

34、i/6; R1=M*abs(x-X(1)*(x-X(2)/2運轉后輸出誤差限為R1 =0.2742.n=2輸入程序輸入程序 X=0:pi/4:pi/2; Y =cos(X) ,l01= conv (poly(X(2),poly(X(3)/( X(1)- X(2)* ( X(1)- X(3), l11= conv (poly(X(1), poly(X(3)/( X(2)- X(1)* ( X(2)- X(3),l21= conv (poly(X(1), poly(X(2)/( X(3)- X(1)* ( X(3)- X(2),l0=poly2sym (l01),l1=poly2sym (l11)

35、,l2=poly2sym (l21),P = l01* Y(1)+ l11* Y(2) + l21* Y(3), L=poly2sym (P),x=pi/6; Y = polyval(P,x)運轉后輸出基函數l01、l11和l21及其插值多項式的系數向量P、插值多項式L和插值Y為l0 =228155022448185/281474976710656*x2-2150310427208497/1125899906842624*x+1l1 =-228155022448185/140737488355328*x2+5734161222659/2251799885248*xl2 =22815502244

36、8185/281474976710656*x2-5734161222659/9007199254740992*xP = -0.3357 -0.1092 1.0000L=-6048395780875/18014398509481984*x2-7870612110600739/72057594037927936*x+1Y = 0.8508輸入程序 M=1;x=pi/6; R2=M*abs(x-X(1)*(x-X(2) *(x-X(3)/6運轉后輸出誤差限為R2 =0.0239.cos (/6)= 0.8660254三、三、n次拉格朗日插值的次拉格朗日插值的Matlab程序程序例例*1-4給出節(jié)點數

37、據給出節(jié)點數據f(-2.00)=17.00， f(0.00)=1.00 ，f(1.00)=2.00， f(2.00)=17.00 ，作三次拉格朗日插值多項式計，作三次拉格朗日插值多項式計算算f(0.6)，并估計其誤差。，并估計其誤差。解：輸入程序解：輸入程序 X=-2,0,1,2; Y =17,1,2,17;p1=poly(X(1); p2=poly(X(2);p3=poly(X(3); p4=poly(X(4); l01= conv ( conv (p2, p3), p4)/( X(1)- X(2)* ( X(1)- X(3) * ( X(1)- X(4), l11= conv ( conv

38、 (p1, p3), p4)/( X(2)- X(1)* ( X(2)- X(3) * ( X(2)- X(4),l21= conv ( conv (p1, p2), p4)/( X(3)- X(1)* ( X(3)- X(2) * ( X(3)- X(4),l31= conv ( conv (p1, p2), p3)/( X(4)- X(1)* ( X(4)- X(2) * ( X(4)- X(3),l0=poly2sym (l01),l1=poly2sym (l11),l2=poly2sym (l21), l3=poly2sym (l31),P = l01* Y(1)+ l11* Y(2)

39、 + l21* Y(3) + l31* Y(4),運轉后輸出基函數l0，l1，l2和l3及其插值多項式的系數向量P為l0 =-1/24*x3+1/8*x2-1/12*x，l1 =1/4*x3-1/4*x2-x+1l2 =-1/3*x3+4/3*x，l3 =1/8*x3+1/8*x2-1/4*xP = 1 4 -4 1輸入程序 L=poly2sym (P),x=0.6; Y = polyval(P,x)運轉后輸出插值多項式和插值為L = x3+4*x2-4*x+1 Y =0.2560.輸入程序 syms M; x=0.6; R3=M*abs(x-X(1)*(x-X(2) *(x-X(3) *(x

40、-X(4)/24運轉后輸出誤差限為R3 =91/2500*M 即 00. 2 ,00. 2,04. 043-fRsyms M為定義一個符號變量M四、拉格朗日多項式和基函數的四、拉格朗日多項式和基函數的Matlab程序程序求拉格朗日插值多項式和基函數的求拉格朗日插值多項式和基函數的MATLAB主程序如下主程序如下function C, L,L1,l=lagran1(X,Y)m=length(X); L=ones(m,m);for k=1: m V=1; for i=1:m if k=i V=conv(V,poly(X(i)/(X(k)-X(i); endendL1(k,:)=V; l(k,:)=

41、poly2sym (V)endC=Y*L1;L=Y*lreturn新建一個lagran1.m文件，將程序拷貝到文件中，并將文件保管到matlab的任務目錄下，即可直接調用函數。例例*1-5給出節(jié)點數據給出節(jié)點數據f(-2.15)=17.03， f(-1.00)=7.24 ，f(0.01)=1.05， f(1.02)=2.03 ， f(2.03)=17.06， f(3.25)=23.05，作五次拉格朗日插值多項式和基函數，作五次拉格朗日插值多項式和基函數，并寫出估計其誤差的公式。并寫出估計其誤差的公式。解：輸入程序解：輸入程序 X=-2.15 -1.00 0.01 1.02 2.03 3.25;

42、Y=17.03 7.24 1.05 2.03 17.06 23.05;C, L ,L1,l= lagran1(X,Y)運轉后輸出五次拉格朗日插值多項式L及其系數向量C，基函數l及其系數矩陣L1如下C = -0.2169 0.0648 2.1076 3.3960 -4.5745 1.0954L =1.0954-4.5745*x+3.3960*x2+2.1076*x3+0.0648*x4-0.2169*x5L1 = -0.0056 0.0299 -0.0323 -0.0292 0.0382 -0.0004 0.0331 -0.7 -0.0503 0.6305 -0.4852 0.0048 -0.0

43、693 0.2184 0.3961 -1.2116 -0.3166 1.0033 0.0687 -0.1469 -0.5398 0.6528 0.9673 -0.0097 -0.0317 0.0358 0.2530 -0.0426 -0.2257 0.0023 0.0049 0.0004 -0.0266 0.0001 0.0220 -0.0002l = -0.0056*x5+0.0299*x4-0.0323*x3-0.0292*x2+0.0382*x-0.0004 0.0331*x5-0.7*x4-0.0503*x3+0.6305*x2-0.4852*x+0.0048 -0.0693*x5+0

44、.2184*x4+0.3961*x3-1.2116*x2-0.3166*x+1.0033 0.0687*x5-0.1469*x4-0.5398*x3+0.6528*x2+0.9673*x-0.0097 -0.0317*x5+0.0358*x4+0.2530*x3-0.0426*x2-0.2257*x+0.0023 0.0049*x5+0.0004 *x4-0.0266*x3+0.0001*x2+0.0220*x-0.0002估計其誤差的公式為 25. 3 ,15. 225. 303. 202. 101. 000. 115. 2! 665-=xxxxxxfxR五、拉格朗日插值及其誤差估計的五、拉

45、格朗日插值及其誤差估計的Matlab程序程序function y,R=lagranzi(X,Y,x,M)n=length(X); m=length(x);for i=1:m z=x(i);s=0.0; for k=1:n p=1.0; q1=1.0; c1=1.0;for j=1:n if j=kp=p*(z-X(j)/(X(k)-X(j); end q1=abs(q1*(z-X(j);c1=c1*j; end s=p*Y(k)+s; end y(i)=s;endR=M*q1/c1;return例例*1-6知知sin30 =0.5，sin45 =0.7071，sin60 =0.8660，用拉格

46、朗日及其誤差估計的用拉格朗日及其誤差估計的MATLAB函數求函數求sin40 的的近似值，并估計其誤差。近似值，并估計其誤差。解：輸入程序解：輸入程序 x=2*pi/9; M=1; X=pi/6 ,pi/4, pi/3;Y=0.5,0.7071,0.8660; y,R=lagranzi(X,Y,x,M)運轉后輸出插值運轉后輸出插值y及其誤差限及其誤差限R為為y = R =0.6434 8.8610e-004.sin (40 )= 0.6427876)()(0010101xxxxyyyxP- - - - = =)(001010 xxxxxxy- - - - = =()fffx0,x1 二次牛頓插

47、值多項式二次牛頓插值多項式 我們再看線性插值的點斜式: )(00 xxy- - = =fx0,x1常數常數(差商差商) 由此啟發(fā)，我們希望二次插值也能類似地有有規(guī)律的由此啟發(fā)，我們希望二次插值也能類似地有有規(guī)律的組合表達式組合表達式:P2(x)=0 + 1(x-x0) + 2(x-x0)(x-x1)利用利用P2(x0)=y0有有: 0 = y0 ,利用利用P2(x1)=y1有有: 1 = 0101xxxx- - -()ff= fx0,x1 ,利用利用P2(x2)=y2有有: 2 = fx0,x1 (x2-x0)(x2-x1) (x2-x0)(x2-x1)0 xx2 - -()ff (x2-x0

48、)-fx0,x2fx0,x1 x2 - x1 =-= fx0,x1,x2 ;P2(x)=f(x0) + (x-x0) + (x-x0)(x-x1) fx0,x1 fx0,x1,x2 fx0,x2 x=x0時0注注: 1. 現實上現實上,從上述可看出二次牛頓插值公式是用待從上述可看出二次牛頓插值公式是用待定系數法求得的定系數法求得的; 2. 它也可看作是三個特殊函數的一種線性組合它也可看作是三個特殊函數的一種線性組合:P2(x)=f(x0) + (x-x0) + (x-x0)(x-x1) fx0,x1 fx0,x1,x2 fx0,x1 , fx0,x1,x2 f(x0), 1 , (x-x0)

49、, (x-x0)(x-x1)即函數 的線性組合,組合系數為 本質上還是基函數法本質上還是基函數法. 更普通地，更普通地，n+1個節(jié)點的插值多項式，我們希望由上個節(jié)點的插值多項式，我們希望由上述類似的一組特殊函數：述類似的一組特殊函數：來線性組合為：來線性組合為： 1 , (x-x0) , (x-x0)(x-x1)，(x-x0)(x-x1)(x-xn).(.)()()(10102010- - - - - - - - - = =nnnxxxxaxxxxaxxaaxN那么其組合系數是什么樣的呢？怎樣求呢？那么其組合系數是什么樣的呢？怎樣求呢？我們同樣可用待定系數法我們同樣可用待定系數法. 容易發(fā)現容

50、易發(fā)現,計算計算a0, a1, a2 , an 是很有規(guī)律的是很有規(guī)律的.一、均差及其性質一、均差及其性質2 牛頓插值牛頓插值當當x=x0時，時，Pn(x0)=a0=f0.當當x=x1時，時，Pn(x1)=a0+a1(x1-x0)=f1, 推得推得a1=f1-f0 x1-x0當當x=x2時，時，Pn(x2)=a0+a1(x2-x0)+a2(x2-x0)(x2-x1)= f2,推得推得f2-f0 x2-x0- -f1-f0 x1-x0a2=x2-x1 依次遞推可得到依次遞推可得到a3, , an. 為寫出系數為寫出系數 ak的普通表達式的普通表達式,先引進如下均差定義先引進如下均差定義. 定義定

51、義2 稱稱 為函數為函數f(x)關于點關于點x0,xk的的一階均差一階均差.稱稱 為為f(x) 的二階均差的二階均差.普通地普通地, 稱稱 為為 f(x) 的的k 階均差階均差(差商差商). fx0,xk =f(xk)-f(x0)xk-x0 fx0,x1,xk=fx0,xk- fx0,x1- fx0,x1xk-x1 fx0,x1,xk=fx0, xk-2,xk- fx0,x1, ,xk-1- fx0,x1, ,xk-1xk-xk-1均差有如下的根本性質均差有如下的根本性質: 1 k 階均差可表示為函數值階均差可表示為函數值f(x0), f(x1), f(xk)的線性組合的線性組合,即即 fx0

52、,x1,xk=f(xj)(xj-xj+1) (xj-xk)(xj-xj-1)(xj-x0) kj=0這個性質可用歸納法證明這個性質可用歸納法證明. 這個性質也闡明均差與節(jié)點的陳列這個性質也闡明均差與節(jié)點的陳列次序無關次序無關,稱為均差的對稱性稱為均差的對稱性,即即 fx0,x1,xk= fx1,x0,x2,xk= = fx1, , xk ,x0 fx0, x1,xk =f(n)()n!,ba 3 假設假設f(x)在在a,b上存在上存在n階導數階導數, 且節(jié)點且節(jié)點x0,x1,xn a,b,那么那么n階均差與導數關系如下階均差與導數關系如下: 這個公式可直接用羅爾定理證明這個公式可直接用羅爾定理

53、證明.，0!)()(0)()(= =- -= =nnnxxfnfq .!)()(0nfxxfnn = =， 所以所以，使使記記為為個個零零點點內內至至少少有有在在理理，可可知知零零點點；反反復復應應用用羅羅爾爾定定個個內內至至少少有有在在個個零零點點，故故的的兩兩個個零零點點間間至至少少有有一一在在，個個零零點點，根根據據羅羅爾爾定定理理上上有有在在處處均均為為零零，所所以以，在在，證證明明：設設, , 1 ,)(,)()()(1,)()( ),()()( )(0100babaxqnbaxqxqxqnbaxqxxxqxxxxxxxxxfxqnnnn - - - -= = fx0,x1,xk=f

54、x1, xk-1,xk- fx0,x1, ,xk-1- fx0,x1, ,xk-1xk-x02 由性質由性質1可得可得:,)()()(000 xxfxxxfxf- - = =,)(,101100 xxxfxxxxfxxf- - = =,.,)(,.,.,0010nnnnxxxfxxxxfxxxf- - = =- -).(.)()()(10102010- - - - - - - - - = =nnnxxxxaxxxxaxxaaxN12 n1依次將后一式代入前一式依次將后一式代入前一式.)(,)(,)()(102100100 - - - - - = =xxxxxxxfxxxxfxfxf).(,.,

55、100- - - - nnxxxxxxf)().(,.,100nnnxxxxxxxxxf- - - - - -Nn(x)Rn(x)ai = f x0, , xi 二、牛頓插值公式二、牛頓插值公式)().(,.,100nnnxxxxxxxxxf- - - -= =- -Rn(x).)(,)(,)(102100100 - - - - - = =xxxxxxxfxxxxfxf).(,.,100- - - - nnxxxxxxfNn(x)n+1(x)10(0nkxxfakk，= = = 多項式多項式Nn(x)顯然滿足插值條件顯然滿足插值條件,即即Nn(xj)=f(xj),(j=1, n),且次數不超越

56、且次數不超越n,由獨一性定理它就是前述的由獨一性定理它就是前述的Ln(x),其系數為其系數為 Nn(x)稱為牛頓均差插值多項式稱為牛頓均差插值多項式,它比拉格朗日插值多項式它比拉格朗日插值多項式計算量省計算量省,且便于程序設計且便于程序設計.注：注： 由獨一性可知由獨一性可知 Nn(x) Ln(x)， 只是算法不同，故只是算法不同，故其他項也一樣，即其他項也一樣，即)(! ) 1()()(,.,1) 1(10 xnfxxxxfnxnnn=),(,!)(,.,maxmin)(0 xxkfxxfkk= 實踐計算過程實踐計算過程為為f (x0)f (x1)f (x2)f (xn1)f (xn)f x

57、0, x1f x1, x2 f xn1, xnf x0, x1 , x2 f xn2, xn1, xnf x0, , xn f (xn+1) f xn, xn+1 f xn1, xn, xn+1 f x1, , xn+1 f x0, , xn+1均差計算可列均差表如下：均差計算可列均差表如下：, 2, 1 , 1, 0 ,11 = = =- - -= =- - iikixxxxfxxfxxfikkikiki 例例1 根據如下函數值表建立不超越根據如下函數值表建立不超越3次的拉格朗日插值多次的拉格朗日插值多項式及牛頓插值多項式項式及牛頓插值多項式Nn(x),并驗證插值多項式的獨一性并驗證插值多項

58、式的獨一性. 解解: (1)拉格朗日插值多項式拉格朗日插值多項式Ln(x).插值基函數插值基函數xk0124f (xk)19233拉格朗日插值多項式為拉格朗日插值多項式為:121445411 )(3)(23)(9)()()(233210303-=xxxxlxlxlxlyxlxLiii,12181241)24)(14)(04()2)(1)(0()(,4541)42)(12)(02()4)(1)(0()(,38231)41)(21)(01()4)(2)(0()(, 1478781)40)(20)(10()4)(2)(1()(233232231230 xxxxxxxlxxxxxxxlxxxxxxxl

59、xxxxxxxl-=-=-=-=-=-=-=-=xkf (xk) 一階均差一階均差二階均差二階均差三階均差三階均差0119822314343- -10- -8(2) 牛頓插值多項式牛頓插值多項式Nn(x).建立如下差商表建立如下差商表牛頓插值多項式為牛頓插值多項式為:121445411 )2)(1)(0(411) 1)(0(3)0(81)(233-=-=xxxxxxxxxxN411-(3) 獨一性驗證獨一性驗證.經過比較牛頓插值多項式和拉格朗日插值多項式經過比較牛頓插值多項式和拉格朗日插值多項式,知知: Nn(x) = Ln(x)這一現實與插值多項式的獨一性一致這一現實與插值多項式的獨一性一致

60、. 例例2留意：當標題中沒有指明用哪一種方法建立插值多項式時，原留意：當標題中沒有指明用哪一種方法建立插值多項式時，原那么上拉格朗日和牛頓插值法都可行。近似計算時，牛頓插值那么上拉格朗日和牛頓插值法都可行。近似計算時，牛頓插值法計算量小，當添加節(jié)點時只需添加一項，因此比較方便。相法計算量小，當添加節(jié)點時只需添加一項，因此比較方便。相對之下，拉格朗日插值法沒有上述優(yōu)點，但它在實際證明上因對之下，拉格朗日插值法沒有上述優(yōu)點，但它在實際證明上因其插值基函數的許多特點而得到廣泛運用。其插值基函數的許多特點而得到廣泛運用。知函數知函數y=f(x)的數據如下表的數據如下表i0123xi0123yi=f(x