Newton插值法的程序設(shè)計與應(yīng)用

1、引言在我們生活中，許多實際問題都都可以利用函數(shù)來表示某種內(nèi)在規(guī)律的數(shù)量關(guān) 系，但是在很多應(yīng)用領(lǐng)域，函數(shù) f(x)有時不能直接寫出表達(dá)式或者過于復(fù)雜不易計 算，而只能給出函數(shù)f(x)在若干個點上的函數(shù)值或?qū)?shù)值，通常也是造一張函數(shù)表， 當(dāng)遇到要求表中未列出的變量的函數(shù)值時， 就必須做數(shù)值逼近.例如給定了函數(shù)f(x) 在幾個特定點上的函數(shù)值，為了研究函數(shù)的變化規(guī)律，就要根據(jù)這個表，將其公式化 尋求某一函數(shù)y(x)去逼近f(x),在給定點上等于函數(shù)值，在其他點上約等于函數(shù)值， 這樣y(x)既能反映f(x)的函數(shù)特性，又便于計算.稱 y(x)的插值函數(shù)，f(x)為被插 信函數(shù).我們可以求一個待定函數(shù)來

2、近似反映函數(shù)的特性，使得待定函數(shù)在給定點上等于 函數(shù)值，在其它點上的函數(shù)的值作為函數(shù)的近似值，這種方法稱為插值法.利用差值函數(shù)很容易得到Lagrange插值多項式，公式結(jié)構(gòu)緊湊，在理論分析中甚為方便.由于 公式中的每一項與所有的插值結(jié)點有關(guān).因此，如果需要增加一個插值結(jié)點，則Lagrange插值公式中的每一項都需要改變計算量大，為了克服這一缺點，于是產(chǎn)生 了 Newton插值法.有的要求更高需要在某些點處的導(dǎo)數(shù)也相同于是前兩種方法達(dá)不 到精度于是產(chǎn)生了 Hermite插值法.本題主要討論牛頓插值多項式.主要是對牛頓插值法的理論總結(jié)，利用C+取序和Matlab語言程序編寫，以及對其應(yīng)用討論.第一

3、章插值方法求得近似函數(shù)的方法通常有兩種：一種是插值法，另一種是曲線擬合法.本章主 要討論插值法用多項式來逼近列表函數(shù)問題，即具有唯一插值函數(shù)的多項式插值,對其中的多項式插值主要討論n次多項式插值的方法，即給定 n+1各點處的函數(shù)值后， 怎樣構(gòu)造一個n次插值多項式的方法.§ 1.1 值問題的由來在許多實際問題及科學(xué)研究中，因素之間往往存在函數(shù)關(guān)系.然而這種關(guān)系經(jīng)常 很難有明顯的解析表達(dá)式，通常只是由觀察與測試得到一些離散數(shù)值 .有時由于表達(dá) 式過于復(fù)雜，不僅使用不便，而且不易于進(jìn)行計算與理論分析 .解決此類問題的方法 有兩種：一種是插值法，另一種是擬合法.插值法是一種古老的數(shù)學(xué)方法，

4、他來自生產(chǎn)實踐，早在一千多年前，我國科學(xué)家 在研究歷法上就應(yīng)用了線性插值與二次插值，但它的基本理論卻是在微積分產(chǎn)生之后才逐漸完善的，具應(yīng)用也逐步增多.特別是在計算機軟件中，許三多庫函數(shù).等的計算 實際上歸結(jié)于它的逼近函數(shù)的計算.在工程的實際問題中，由于航空、造船、精密機 械加工等實際問題的需要，插值方法在實踐上和理論上顯得更為重要，并得到了進(jìn)一 步的發(fā)展和廣泛的應(yīng)用.因此我們希望能夠得到一個“簡單函數(shù)”逼近被計算函數(shù)的 函數(shù)值.于是就得到了一種方法.叫做插值逼近或者插值法.§ 1.2 個基本概念§ 1.2.1 值函數(shù)的概念及相關(guān)性質(zhì)定義1.1 設(shè)連續(xù)函數(shù)y - f(x)在區(qū)

5、間a,b上有定義，已知在n+1個互異的點 Xo,X1，Xn上取值分別為y°, y1,yn (設(shè)a Mx。x 4 «b ).若在函數(shù)類中存在 以簡單函數(shù)p(x),使得p(xi)= yi (i =0,1，n ),則稱p(x)為f (x)的插值函數(shù).稱Xo,X1，Xn為插值節(jié)點,稱a,b為插值區(qū)間.定義1.2設(shè)如下表XiXoX1 .XnV、VoV1 .yn求一個 p(x),使其滿足：(1)p(x) w IPn(x) P(Xi) = f (Xi) = yi , i =0,1，n則稱p(x)為f (x)的n次代數(shù)插值多項式，該問題稱為n次代數(shù)插值問題.定理1.1n次代數(shù)插值問題的解存

6、在且唯§ 1.2.2 Lagrange插值函數(shù)的概念及相關(guān)性質(zhì)定義1.3 若n次多項式1k(x) (k=0,1n)在n+1個互異節(jié)點x0,x1，xn上滿 足%, j = k兀(七)=跖Jk,j=0,1;,n)P j # k則稱1k(x) ( k =0,1，n )為節(jié)點x0,x1，xn上的n次插值基函數(shù).n次差值多項式Ln(x)為nLn(x)=" yJk(x)k =0稱之為拉格朗日多項式.定理1.2 設(shè)fwCn由a,b,則余項為:Rn(x)=fn1() (n 1)!Wn 1(x)其中 Wn i(X)=(x -Xo)(X -Xi) (x -xn).§ 1.2.3 商的

7、概念及相關(guān)性質(zhì)定義1.4稱函數(shù)f (x)于點Xi,Xj ( x #Xj)上的平均變化率為f(x)在Xi,Xj處的一階差商，記作fx,Xj,即fX,Xjf(Xj)-f(x)Xj -XiXi 二 Xj .一階差商的差商為f(x)在點xixXk處的二階差商.記作flxixXk,即fXj,Xk 一 fXi,Xjfx,Xj,Xk:Xk -Xi般地，在求出f (x)的m -1階差商后就可以構(gòu)造f(x)的m階差商.稱f Xi0,Xi1, . ,XimfXi1,% - fXi0,xim . xi0為f (x)在為0?小xim處的m階差商.特別地，規(guī)定零階差商fxj = f (xi).fXo,Xi,X2二性質(zhì)1

8、.1任意調(diào)換插值節(jié)點的次序其值不變即fXi,X0,X2 = fXo,X2,Xi性質(zhì)1.2差商是一個數(shù)，可表示為某個函數(shù)的線性組合.性質(zhì)1.3 當(dāng)f(X)的m階導(dǎo)數(shù)f(m)(X)在包含節(jié)點Xi0,Xi1，Xm的某個區(qū)間上存在時，在Xi0, 4,Xim之間至少存在一點之使一，f(m)()f Xi 0 , X 1, Xim =m!注1.1由性質(zhì)2可以推出，n次多項式的一階差商是n-1多項式.表1-1差商表Xkf(Xk)一階差商二階差商三階差商四階差商X0f(Xk)X1f(Xk)fx。, xX2f (Xk)fX1,X2f Xo,x,x2X3f(Xk)f X2,X3f X1, X2,X3fx,x,%,x

9、X4f (Xk)f X3,X4f X2, X3,X4f X1,X2, X3f%,x,x,x,x§ 1.2.3 Newton插值函數(shù)與差分 的概念及相關(guān)性質(zhì)定義1.5我們定義牛頓插值公式為Nn(X)= f (X。)+ f Xo, x(x - X。)+ f Xo,X1,x2(x- Xo)(x- X1)(11)f X0,X1,Xn(X - X0)(X - X1) (X -Xni).定義1.6 設(shè)函數(shù)y = f (X)在等距節(jié)點Xk =x。+kh (k =0,1，n)上的函數(shù)值yk=f(x。(k=0,1，n),其中h為常數(shù)，稱為步長.稱函數(shù)在每個小區(qū)間上的增量y-yk為函數(shù)f(x)在Xk處以

10、h為步長的一階向前差分，記作 ",即-7k =yk 1 -yk稱一階差分的差分Ayk書-Ayk為二階差分，記作Hyk,即2yk =AykHtAyk = ykt2yk + yk由此可得m-1階差分的差分myk=”'yk1 一即“為f(x)在Xk點處的m階差分.稱yk-yg為函數(shù)f(x)在Xk處以h為步長的一階向后 差分，記作''yk.稱Vmbk Vmyj為函數(shù)f(x)在Xk處以h為步長的m階向后差分，記作Vmyk.即，r m, m 1I m' yk=yk-ykj表1-2向前差分表XiV、y、2y3yKvXoy1Xiy1y0X2y2K y。X3y3y22w

11、3y0X4VaK y2Nv。表1-3向后差分表X、y、y、V2ybybyXoy1X1y1VyoX2y2V2y0X3y3內(nèi)2可礦y0X4Va內(nèi)2五丁y0性質(zhì)1.4各階差分可以表示成函數(shù)值的線性組合.即n nyj = " (-1)kC：yj nA k =0 n(2) “j =" (7)kC：yj/k =0性質(zhì)1.5差商與差分有如下關(guān)系fXo,Xi,Xnnfo n!hnfon!hn定義1.7 設(shè)n十1個等距插值節(jié)點的順序為Xo,Xi，2二,4,則由差商表示的牛頓插值公式為Nn(X) = fXo fX),X1(X-Xo) fX),X1,X2(X-Xo)(X-X1) fXo,X1,

12、,Xn(X-Xo) (X- XnA)又設(shè)插值節(jié)點x = x0th (0 ： t ： n)由性質(zhì)1.5差商與差分關(guān)系fX0,Xi； ,Xk:可得N n(X)= fo.2f. . t(t -1) (t-n 1) .nf!_ f o '- Ton!n=zk ZZ0t(t 1) (t k 1)(1.2)k!稱為牛頓前插公式.其余項公式為fnZ )n1 _Rn(x)= Rn(Xo th)= -(-t(t -1) (t-n)h ,Xo,Xn(n 1)!定義1.8 設(shè)n+1個等距插值節(jié)點的順序為4,4,X1,X0,則由差商表示的牛頓 插值公式為Nn(X)=fXn fXn,Xn(XXn)fXn,XnL

13、Xn/(XXn)(XXn)(XXn/),f Xn , Xn 1, Xo( X -Xn ) (X -X1 )又設(shè)插值節(jié)點x = Xn -th (o : t :二 n)由性質(zhì)1.5差商與差分關(guān)系fXn,Xn4, Xn-k=fnk!hk可得、 t(t -1)、2n t(t -1) (t -n 1) nNn(x) = fn Yfn ' J fn -( -1)fn(1.3)2!n!= 1(7)5(tk1)：kfn小k!稱為牛頓后插公式.其余項公式為fn1( )n1 n1 .R(X) =Rn(Xn-th)t(t-1) (t-n)(-1) h ,Xo,Xn(n 1)!第二章Newton插值法的程序設(shè)

14、計與應(yīng)用利用插值基函數(shù)很容易得到 Lagrange插值多項式，公式結(jié)構(gòu)緊湊，在理論分析 中也很方便，但當(dāng)插值節(jié)點增減時全部插值基函數(shù)L(x) i =0,1，n都要隨之變化，整個公式也將發(fā)生變化，這在實際計算中是很不方便的，于是我們導(dǎo)出插值多項式的 又一種表示形式一一Newton插值公式.遠(yuǎn)遠(yuǎn)減少了工作量.但由于通過相同結(jié)點的差 值多項式是惟一的，因此，拉格朗日差值公式與牛頓公式本質(zhì)上是一個多項式，只不過形式不同而已.為了更方便的應(yīng)用Newton插法，本章主要討論Newton插值法的程序設(shè)計，并對其 應(yīng)用進(jìn)行討論.§ 2.1 Newton插值法的算法與程序設(shè)計本節(jié)主要討論n次多項式插值

15、的方法，即給定n+1各點處的函數(shù)值后，怎樣構(gòu)造 一個n次插值多項式的方法.為了更方便的計算并與計算機結(jié)合，我們便作出了相應(yīng) 的程序.并利用Matlab進(jìn)行繪圖. 2.1.1 2.1.1 Newton插值法的應(yīng)用步驟步驟 首先我們按照表1-1,求得各點的差商.然后利用Newton前差或后差公式 即(公式1.3)把數(shù)值帶入.即可以求得n次多項式.它在計算機上的應(yīng)用步驟如下：步驟1輸入所要求的牛頓多項式的次數(shù)n ,并依次輸入n+1個節(jié)點(Xi,yi). for : i=0 to n+1scanf("%f',&xi) scanf("%f",&y0i

16、);步驟2計算各階差商 for : i=1 to n+1for : j=i to n+1 if(i>1)yij=(y1j-y11)/(xj-xi);elseyij=(yi-1j-yi-1j-1)/(xj-xj-1);步驟3代入牛頓插值公式，可計算得出結(jié)果.for : i=1 to n+1temp=temp*(xx-xi-1);newton=newton+yii*temp;由于篇幅所限，流程圖、C語言程序和Matlab程序詳見附錄(附錄A、B、C). 2.1.2 利用Newton插值法程序的實例為了更方便的應(yīng)用牛頓插值法，我們進(jìn)行了與計算機的結(jié)合，下面我們將展示幾 個例子.例2.1已知y

17、= sin x的一組數(shù)據(jù)為0冗兀冗JIx6432sinx01V2V31222(1)構(gòu)造牛頓插值函數(shù)并作圖分析.(2)并分別利用程序估計sin-, sin啊的估計值.1215分析 首先我們可以通過程序求出差商表：xkf (xk)一階差商二階差商三階差商四階差商006120.954929268714三 20.791089631-2.0860762L3而 20.607024424-0.35153865-0.13648909匹 210.25587263-0.44710035-0.091254700.028797106帶入定義1.5中可求得如下牛頓插值多項式如下冗N2(x) =0.95492926x -

18、0.2086076x(x -)(2.1)JIJIN3(x) =0.95492926x -0.2086076x(x - ) -0.13648909x(x -一)(x -一) 664(2.2)71JE31N4(x) = 0.95492926x -0.2086076x(x ) - 0.13648909x(x-)(x -)(2.3)664JIJJ0.028797106x(x -)(x -)(x -)643第二步 我們分別利用C+?序和Matlab分別計算sin工和sin粵的值.步驟如下： 1215利用C+E序：首先我們輸入所要求的牛頓多項式的次數(shù)n,然后/&入n+1個節(jié)點的值.即可以得出sin

19、工和sin巴的值為0.2586和0.0804 ; 1215利用Matlab程序：首先我們在Matlab中新建一個文件，文件內(nèi)容即為附錄C中Matlab程序代碼，然后輸入以下步驟：x=0,pi/6,pi/4,pi/3,pi/2;y=0,0.5,sqrt(2)./2,sqrt(3)./2,1;x0=pi/12;Newton(x,y,x0ans =0.2586x=0,pi/6,pi/4,pi/3,pi/2;y=0,0.5,sqrt(2)./2,sqrt(3)./2,1;x0=14.*pi/15;Newton(x,y,x0) ans =0.0804 即結(jié)果為0.2586和0.0804 利用Matlab

20、分別對求得的牛頓差值多項式作圖(代碼詳見附錄D)圖形如下.-1%IIM2X ifcpi/12.O.250e)閭i-如g)XI1-D511.52圖2-1例2.2設(shè)f (x) =ln(1 +x)的函數(shù)表如下:x0.250.300.360.390.450.2231440.2623640.3074850.3293040.371564試計算 ln(1.275) , f(x) =ln2.3分析：同上題步驟我們先求差商表，進(jìn)而代入公式可得N4(x) =0.22314 0.7844(x-0.25) 0.294393939(x - 0.25)(x - 0.30)+ 0.1411736357(x-0.25)(x-

21、0.30)(x-0.36)+0.057720032( 2.4 )(x -0.25)(x -0.30)(x-0.36)(x -0.39)N3(x) = 0.22314 + 0.7844(x 0.25) -0.294393939(x-0.25)(x-0.30)250.1411736357(x -0.25)(x -0.30)(x -0.36)M(x) =0.22314 +0.7844(x0.25) -0.294393939 (x - 0.25)(x - 0.30)(2.6)利用C+?序我們可以得到計算結(jié)果：W(0.275) =0.242946,電(0.275) = 0.242945,N2(0.275

22、) =0.242938 N4(1.3) = 0.928823 , N3(1.3) = 0.877002 , N2(0.275) = 0.737653我們利用Matlab進(jìn)行繪圖.1.J上上1一0 .cP :.遍NFSJ/i明出處100,511,5圖2-2從上述兩個例子我們可以看出，多項式(2.1 )在區(qū)間0,1周圍與原函數(shù)逼近的 較好.離這個區(qū)間越遠(yuǎn)與原函數(shù)的誤差越大在 x =1.5處時，該圖像就已經(jīng)開始背離圖 像了.所以該多項式只能在一個小白區(qū)間里可以逼近原函數(shù)，不適合作為原函數(shù)的逼近函數(shù).也可以看出多項式(2.2)在0,；區(qū)間的周圍逼近的較好，但是x = 2處時， 該圖像就離原圖像誤差較大

23、.多項式(2.3)在區(qū)間0,2.5都逼近的挺好，從圖中我 們看出在遠(yuǎn)離這個區(qū)間的圖像誤差相對較大，但是在這三個多項式中是逼近最好的，同樣在例2.2中也是如此.于是我們我們可以得出節(jié)點越多，函數(shù)逼近的相對較好.在節(jié)點附近逼近的越好， 越遠(yuǎn)離節(jié)點誤差越大，所以公式適用于計算節(jié)點附近的值于是為了減小誤差，在下一節(jié)的等距節(jié)點下的插值公式根據(jù)所求的點的函數(shù)值的不同分別采取了前插和后插公 式.§ 2.2 距節(jié)點下的Newton插值算法與程序設(shè)計前面我們講述了一般節(jié)點下的牛頓插值公式，為了計算方便于是有了對等距節(jié)點下的牛頓多項式的研究，本節(jié)將對等距節(jié)點下的插值多項式進(jìn)行總結(jié)討論.§2.2

24、.1 等距節(jié)點下的Newton插值法的程序算法步驟步驟 我們按照表1-2或1-3 ,求差分.然后利用Newton前插公式(公式1.2)或 Newton后插公式(公式1.3)并把數(shù)值帶入.即可以求得n次多項式.它在計算機上的應(yīng)用步驟如下：步驟1輸入所要求的牛頓多項式的次數(shù)n，步長h，并依次輸入n+1個節(jié)點(Xi, yi).for : i=0 to n+1scanf("%f',&xi) scanf("%f',&y0i);步驟2求得各界差分for : i=1 to n+1 for : j=i to n+1yij=(yi-1j+1-yi-1j-1)/

25、(xj-xj);/求向前差分for : i=1 to n+1 for : j=i to n+1yij=(yi-1j-yi-1j-1)/(xj-xj-i)；/求向后差分步驟3代入牛頓插值公式，可計算得出結(jié)果for(i=1;i<n+1;i+) temp=temp*(t-i+1)/i);newton=newton+yii*temp;printf(" 求得的結(jié)果為：N(%.4f)=%9fn",xx,newton);/求得運用前插公式的值for(i=n;i<0;i-)temp=temp*(t-i+1)/i);newton=newton+yn-in-i*temp;print

26、f("求得的結(jié)果為：N(%.4f)=%9fn",xx,newton);/求得運用后插公式的值由于篇幅所限C語言程序(詳見附錄D E).§222等距節(jié)點下Newton插值的實例例2.4 已知tan x的值列表如下:x1.31.311.321.33tanx3.60213.74713.90334.0723近似計算 tan 1.325, tan1.305.采用Newton向后插公式.為此，做差分表xifiV2fiV3fi1.33.60210.14500.01120.00161.313.74710.15620.01281.323.90330.16901.334.0723從而

27、，有p3(x) =4.0723 0.1690t 00128t(t 1) 0.0016 t(t 1)(t 2).2!3!1.325 -1.33將t = -0.5代入上式,得0.01tan1.325 p3(1.325) =3.9869.將t = 1.305 -1.33 =_2.5帶入后插多項式中可以得到0.01tan1.305 5(1.305) =3.0797現(xiàn)在我們利用C+E序首先輸入所求插值的次數(shù)3和步長0.01.然后輸入各個節(jié)點，并輸入所要求的點1.325既可以求出該點的函數(shù)值.即f (1.325) =3.9869. f (1.26) ： p3(1.26) =3.0797卜面我們可得圖圖2-

28、4例2.4設(shè)函數(shù)在各節(jié)點的取值如下X00.20.40.60.81.0f兇1.00.8187310.6703200.5488120.4493290.367879試?yán)貌逯倒角骹 (0.3)的值.采用Newton向前插公式，同上題我們先做差分表，然后相應(yīng)帶入到差分公式 . 中求得后插公式.利用C語言程序步驟如下：首先輸入所求插值的次數(shù) 5和步長0.2.然后輸入各個節(jié)點，并輸入所要求的點 0.3既可以求出該點的函數(shù)值.即f(0.3) =0.7880.由以上例子我們看到例1用了牛頓后插公式，例2用了牛頓前插公式，我們該怎 樣選取.這個經(jīng)過驗證得出，如果所要求的點較靠近節(jié)點 %,則采用前插公式；如果

29、靠近Xn ,則采用牛頓后插公式.附錄A:牛頓插值法的流程圖Ni<nYj<n?i>1輸入i<n+1輸入n/分別輸入 n+1/個輸入各節(jié)點YNNYY/輸出n的值yi-ij -yi-i,j-iy Xj - Xjyi-1j - yi-1,j-1 yj 二Xj - jt=t*(x- Xi J )N=N+y iii=1;t=1;N= y°o開始i=i+1i=i+1j=j+1j=ii=0附錄B:牛頓插值法的程序#include<stdio.h>void main()float x11,y1111,xx,temp,newton;int i,j,n；printf(&

30、quot;Newton插值:n請輸入要運算的值:x=");scanf("%f",&xx);printf("請輸入插值的次數(shù)(n<11):n=");scanf("%d",&n);printf("請輸入 %dffl 值:n",n+1);for(i=0;i<n+1;i+) printf("x%d=",i);scanf("%f",&xi);printf("y%d=",i);scanf("%f",&a

31、mp;y0i);for(i=1;i<n+1;i+)for(j=i;j<n+1;j+) if(i>1)yij=(yi-1j-yi-1j-1)/(xj-xj-i);elseyij=(yi-1j-yi-1j-1)/(xj-xj-1);printf("%fn",yii);temp=1;newton=y00;for(i=1;i<n+1;i+) temp=temp*(xx-xi-1);newton=newton+yii*temp;printf("求得的結(jié)果為：N(%.4f)=%9fn",xx,newton);附錄C:牛頓插值法的Matlab程

32、序function f = Newton(x,y,x0)syms t;if(length(x) = length(y)n = length(x);c(1:n) = 0.0;elsedisp('x 和y的維數(shù)不相等！');return;endf = y(1);y1 = 0;l = 1;for(i=1:n-1)for(j=i+1:n)y1(j) = (y(j)-y(i)/(x(j)-x(i);endc(i) = y1(i+1);l = l*(t-x(i);f = f + c(i)*l;simplify(f);y = yi;if(i=n-1)if(nargin = 3)f = sub

33、s(f,'t',x0);elsef = collect(f);%將插值多項式展開f = vpa(f, 6);endendend附錄D:等距節(jié)點下白前插公式的C語言程序#include<stdio.h> void main() float x11,y1111,xx,temp,newton,t,h;int i,j,n；printf("Newton插值:n請輸入要運算的值:x=");scanf("%f",&xx);printf("請輸入插值的次數(shù)(n<11):n=");scanf("%d",&n);printf("步長:n請輸入要運算的值：h=");scanf("%f",&h);printf(" 請輸入 %dffl 值:n&q