實驗名稱插值法

1、探索實驗5 插值法一、 實驗?zāi)康?了解插值問題及其適用的場合，理解并掌握常用的插值算法的構(gòu)造和計算，了解差商概念、Runge現(xiàn)象及樣條插值方法，學(xué)習(xí)用計算機(jī)求近似函數(shù)的一些科學(xué)計算方法和簡單的編程技術(shù)。 二、概念與結(jié)論1. 插值問題與插值函數(shù)： 由實驗或測量的方法得到所求函數(shù) y=f(x) 在互異點x0 , x1, . , xn 處的值 y0 , y1 , , yn ,構(gòu)造一個簡單函數(shù) j(x) 作為函數(shù) y=f(x) 的近似表達(dá)式 y= f(x) » j(x)使 j(x0)=y0 , j(x1)=y1 , ¼, j(xn)=yn ,(1)這類問題稱為插值問題。 f(x)

2、稱為被插值函數(shù)，j(x) 稱為插值函數(shù)， x0 , x1, . , xn 稱為插值節(jié)點。(1)式稱為插值條件。常用的插值函數(shù)是多項式函數(shù)。且當(dāng)n=1時是稱為線性插值，n=2時稱為Simpson插值或拋物線插值。2插值定理： 假設(shè)x0 ,x1,xn 是n+1個互異節(jié)點,函數(shù)f(x)在這組節(jié)點的值f(xk)(k=0,1,n)是給定的，那么存在唯一的n 次次多項式pn (x)滿足pn (xk)=f(xk), k=0,1,n3.插值的截斷誤差 設(shè)jn(x)是過點x0 ，x1 ，x2 ，xn的 n 次插值多項式， f(n+1)(x)在（a，b）上存在，其中a，b是包含點x0 ，x1 ，x2 ，,xn的任

3、一區(qū)間，則對任意給定的xÎa，b，總存在一點xÎ（a，b）（依賴于x）使 其中wn+1(x)=(x x0) (x - x1)(x-xn ) ，f(n+1)(x) 是f(x)的n+1階微商在 x 的值。 4. 差商：給定一個函數(shù)表 x | x0 x1 . xn - - y | y0 ,y1 . yn其中當(dāng)i ¹j時 ,xi¹xj 記 fxi=f(xi) ,定義f(x)關(guān)于xi,xj的一 階差商 一般的, f(x)關(guān)于xi,xi+1,xi+k的k 階差商定義為:5. 分段線性插值： 在區(qū)間a,b上給定一組節(jié)點: a=x0<x1<x2<<

4、;xn=b且在節(jié)點上的函數(shù)值f(xi )=yi, i=0,1,2,n, 如果函數(shù)j(x) 滿足：1) j (x)在a , b上連續(xù)2) j (xi )=yi, i=0,1,2,n3) j (x)在每個子區(qū)間xi , xi+1(i=0,1,2,n-1)上是次數(shù)為1的多項式則稱j (x)是定義在a ,b上的分段線性插值函數(shù)。6. m次樣條插值函數(shù)： 設(shè)f(x)是區(qū)間a,b上的一個連續(xù)可微函數(shù),在區(qū)間a,b上給定一組節(jié)點: a=x0<x1<x2<<xn=b且y = f(x)在點 x0，x1，x2, xn的值為y0，y1，y2, yn，如果函數(shù)s(x)滿足條件 (1) s(x)

5、在每個子區(qū)間xi , xi+1(i=0,1,2,n-1)上是次數(shù)不超過m的多項式;(2) s(x)在區(qū)間a , b上有m-1階連續(xù)導(dǎo)數(shù);(3) s(xi)=f(xi) =yi , i=0,1,2,n則稱s(x) 為函數(shù)f(x)的在a ,b上的m次樣條插值函數(shù)。x0，x1，x2, 稱為樣條結(jié)點,其中x1,xn-1稱為內(nèi)結(jié)點, x0 , xn 稱為邊界結(jié)點。當(dāng)m=3時,則稱s(x)為函數(shù)f(x)的三次樣條插值函數(shù)。三、程序中Mathematica語句解釋：1. Product f(i) , i ，imin，imax，h 計算乘積 f(imin)´f(imin +h) ´f(im

6、in +2h) ´´f(imin +nh) imax h £ imin + nh £ imax , h>02Simplifyexpr 對expr進(jìn)行化簡3. Appendlist, elem 在表list的后面加一個元素 elem 4. Prependlist, elem 在表list的前面加一個元素 elem四、方法與程序 插值法是函數(shù)逼近的重要方法之一, 它是求近似函數(shù)的一種方法，有著廣泛的應(yīng)用。插值法有很多種,其中以拉格朗日(Lagrange)插值和牛頓(Newton)插值為代表的多項式插值最有特點,常用的插值還有Hermit插值,分段插值和

7、樣條插值等。這里只給出Lagrange插值、Newton插值 、分段線性插值和樣條插值的構(gòu)造過程及程序。1Lagrange插值 Lagrange插值是將待求的n次多項式插值函數(shù)Pn(x)改寫成另一種表示方式，再利用插值條件確定其中的待定插值基函數(shù)，從而求出插值多項式。Lagrange插值是多項式插值，它成功地用構(gòu)造插值基函數(shù)的方法解決了求多項式插值函數(shù)出現(xiàn)的病態(tài)問題。1.1 Lagrange插值的構(gòu)造過程： 令插值多項式為 Ln(x)=l0n(x)y0 + l1n(x)y1 + l2n(x)y2 + lnn(x)yn 式中 l0n(x), l1n(x), l2n(x), ,lnn(x) 都是次

8、數(shù)£ n 的待定多項式，利用插值條件 Pn(x0)=y0 , Pn(x1)=y1 , . , Pn(xn)=yn可以得到 l0n(x0) = 1, l1n(x0) = 0, l2n(x0) = 0, , lnn(x0) = 0 l0n(x1) = 0, l1n(x1) = 1, l2n(x1) = 0 , , lnn(x1) = 0 . l0n(xn) = 0, l1n(xn) = 0, l2n(xn) = 0, , lnn(xn) =1考慮函數(shù) l0n(x),由其在 x1 , x2 . xn 的零值，可以知道它含有 x-x1 , x-x2 ,. ,x-xn n個因子，于是可令 l0

9、n(x) = a(x-x1)(x-x2). (x-xn)利用 l0n(x0)=1 確定其中的系數(shù)a，得 (x-x1)(x-x2)(x-xn )l0n(x) = - (x0-x1)(x0-x2 )(x0-xn ) 類似的可以得出 l1n(x) , l2n(x), lnn(x) (x-x0)(x-x1) (x-xi-1 ) (x-xi+1) (x-xn )lin(x)= - (xi-x0)(xi-x1 )(xi-xni-1 ) (xi-xi+1) (xi-xn ) l0n(x), l1n(x) , lnn(x)稱為以x0 , x1 ,. ,xn為節(jié)點的插值基函數(shù)。于是，將這些插值基函數(shù)代入前面的多

10、項式Ln(x)即得到如下所求的n次Lagrange插值多項式： 1.2 求n次Lagrange插值多項式算法： 1 輸入n+1個插值點： (xi, yi),i=0,1,n2 計算插值基函數(shù)l0n(x), l1n(x) , lnn(x)3 給出n次Lagrange插值多項式：Ln(x)= y0 l0n(x)+ y1 l1n(x) +yn lnn(x)1.3 求Lagrange插值多項式程序：Clearlag,xi,x,yi;xi=Input"xi="yi=Input"yi="n=Lengthxi-1;p=Sumyii*(Product(x-xij)/(xi

11、i-xij),j,1,i-1*Product(x-xij)/(xii-xij),j,i+1,n+1),i,1,n+1;lagx_=Simplifyp說明：本程序用于求n次Lagrange插值多項式。程序執(zhí)行后，按要求通過鍵盤輸入插值基點xi:x0 , x1, . , xn 和對應(yīng)函數(shù)值yi: y0 , y1 , , yn 后，程序即可給出對應(yīng)的n次Lagrange插值多項式lagx。程序中變量說明：xi：存放插值基點x0 , x1, . , xn yi: 存放對應(yīng)函數(shù)值y0 , y1 , , ynlagx: 存放求出的n次Lagrange插值多項式Ln(x)注：語句lagx_=Simplify

12、p用簡化形式給出對應(yīng)的n次Lagrange插值多項式。 1.4 例題與實驗 例1.給定數(shù)據(jù)表x 0 1 2 3y=f(x) 1 3 5 12用Lagrange插值法求三次插值多項式,并給出函數(shù)f(x)在x =1.4的近似值。解： 執(zhí)行Lagrange插值程序后，在輸入的兩個窗口中按提示分別輸入0, 1, 2, 3、1, 3, 5, 12，每次輸入后用鼠標(biāo)點擊窗口的“OK”按扭，得如下插值函數(shù)： 6 + 22 x - 15 x2 + 5 x3- 6所以得到三次插值多項式L3(x)=1+11 x/3-5 x2/2+5 x3/6接著鍵入： lag1.4得輸出:3.52，因此f(x)在x =1.4的近

13、似值為3.52，即f(1.4)»3.52.2. Newton插值 Newton插值是多項式插值的另一種表示形式，它在增加插值節(jié)點時具有靈活性，只要在原有的插值多項式中增加一項就得到新的插值函數(shù)。2.1 Newton插值構(gòu)造過程： 為了使插值Newton多項式具有承襲性，令插值函數(shù)具有下列形式：式中稱為Newton插值基函數(shù)。為求出Nn(x),利用插值條件和差商，得出如下n次Newton插值多項式：計算 Newton插值多項式Nn(x)常用如下差商表：2.2 求Newton插值多項式算法： 1. 輸入n+1個插值點： (xi, yi),i=0,1,n2. 計算差商表3.給出n次Newt

14、on插值多項式。2.3 求Newton插值多項式程序：Clearnewt,s,x;xi=Input"xi="yi=Input"yi="n=Lengthxi;(*計算差商表*)f=Table0,n,n;Dofi,1=yii,i,1,nDofi,j+1=(fi,j-fi+1,j)/(xii-xii+j), j,1,n-1,i,1,n-jPrint"差商表"DoPrintxii," ",fi,i,1,n(*求Newton插值多項式*)fa=1;s=f1,1;Dofa=(x-xik)*fa;s=s+fa*f1,k+1,k,

15、1,n-1newtx_=sSimplify%說明：本程序用于求n次Newton插值多項式。程序執(zhí)行后，按要求通過鍵盤輸入插值基點xi:x0 , x1, . , xn 和對應(yīng)函數(shù)值yi: y0 , y1 , , yn 后，程序依次給出輸入的數(shù)據(jù)表、計算出的差商表、Newton插值多項式、Newton插值多項式的簡化形式。程序中變量說明：xi：存放插值基點x0 , x1, . , xn yi: 存放對應(yīng)函數(shù)值y0 , y1 , , ynf:存放函數(shù)值y0 , y1 , , yn及所有差商newtx: 存放求出的n次newton插值多項式Nn(x)注：1)語句f=Table0,n,n用于產(chǎn)生一個n&

16、#180;n的矩陣變元用于存放函數(shù)值y0 , y1 , , yn及所有差商。2)在Mathematica中有一個求n次插值多項式的命令，命令形式為： InterpolatingPolynomialx01,y0，x1,y1,x2,y2,xn,yn, x 它可以求過n+1個插值點x01,y0，x1,y1,x2,y2,xn,yn的n次插值多項式Pn(x)。2.4例題與實驗例2.給定數(shù)據(jù)表 x 4.0002 4.0104 4.0233 4.0294 y 0.6020817 0.6031877 0.6045824 0.60524041) 計算差商表2) 用Newton插值法求三次插值多項式Nn(x)3)

17、 求f(4.011)的近似值解： 執(zhí)行Newton插值程序后，在輸入的兩個窗口中按提示分別輸入4.0002, 4.0104, 4.0233, 4.0294、0.6020817, 0.6031877, 0.6045824, 0.6052404,每次輸入后用鼠標(biāo)點擊窗口的“OK”按扭，得如下插值函數(shù)：4.0002, 4.0104, 4.0233, 4.02940.6020817, 0.6031877, 0.6045824, 0.6052404 差商表4.0002 0.6020817, 0.108431, -0.0136404, 0.02116294.0104 0.6031877, 0.108116

18、, -0.0130225, 04.0233 0.6045824, 0.107869, 0, 04.0294 0.6052404, 0, 0, 00.6020817 + 0.108431 (-4.0002 + x) - 0.0136404 (-4.0104 + x) (-4.0002 + x) + 0.0211629 (-4.0233 + x) (-4.0104 + x) (-4.0002 + x)-1.41642 + 1.23926 x - 0.268313 x2 + 0.0211629 x3于是我們得到本題的差商表：xi yi f, f , , f , , ,4.0002 0.6020817

19、, 0.108431, -0.0136404, 0.02116294.0104 0.6031877, 0.108116, -0.01302254.0233 0.6045824, 0.107869 4.0294 0.6052404 和三次插值多項式N3(x)= -1.41642 + 1.23926 x - 0.268313 x2 + 0.0211629 x3接著鍵入： newt4.011得輸出0.603253，因此f(4.011)» 0.603253例3. 多項式插值的誤差估計式中可以看到，當(dāng)插值節(jié)點越多時誤差會越小，這個結(jié)論正確嗎？通過實驗說明該結(jié)論的正確性。 解: 考慮函數(shù)f(x)

20、 = (1+x2)-1 在區(qū)間-4,4內(nèi)選取不同個數(shù)的等距插值節(jié)點做觀察，這里分別選-4,4 內(nèi)的9個和11個的等距節(jié)點來做實驗，將對應(yīng)的插值函數(shù)圖與被插函數(shù)f(x) = (1+x2)-1畫在一起做觀察，為簡單起見，這里用Mathematica 命令做實驗，對應(yīng)命令為： In1:= u=Tablex,(1+x2)-1,x,-4,4 ; （*采取f(x) 在-4,4 內(nèi)的9個插值點 In2:= g=ListPlotu, PlotStyle->PointSize0.04 （*將散點圖圖形文件存放在變量g中In3:= s=InterpolatingPolynomialu , x ; （*將插值

21、函數(shù)存放在變量s中In4:= t= Plots, (1+x2)-1, x,-4,4, PlotStyle->Thickness0.005,Thickness0.006 （*將插值函數(shù)s與f(x)畫在一起的圖形文件存放在變量t中In5:= Showt, g （*將散點圖, 插值函數(shù)s, f(x)畫在一個坐標(biāo)系中圖 1 在-4,4中選9個等距節(jié)點的插值函數(shù)與被插函數(shù)圖，粗線為被插函數(shù)圖In6:= u=Tablex,(1+x2)-1,x,-4,4，0.8 ; （*采取f(x) 在-4,4 內(nèi)的11插值點 In7:= g=ListPlotu, PlotStyle->PointSize0.0

22、4 In8:= s=InterpolatingPolynomialu , x ; In9:= t= Plots, (1+x2)-1, x,-4,4, PlotStyle->Thickness0.005,Thickness0.006 In10:= Showt, g 圖 2 在-4,4中選11個等距節(jié)點的插值函數(shù)與被插函數(shù)圖 圖中細(xì)的曲線是插值函數(shù), 粗的曲線是原來函數(shù)(被插函數(shù))。從圖形上觀察, 說明增加插值節(jié)點不能保證使插值誤差變小，本題圖象可以看到，節(jié)點的增多得到的插值函數(shù)產(chǎn)生的誤差在靠近端點的函數(shù)值不但沒有減少反而增大了。因此當(dāng)插值節(jié)點越多時誤差會越小，這個結(jié)論不正確。這個實驗也說明

23、高次插值函數(shù)效果不好。插值點個數(shù)增加，出現(xiàn)誤差增大的現(xiàn)象稱為Runge現(xiàn)象。Runge現(xiàn)象的發(fā)現(xiàn)使人們知道了高次插值函數(shù)誤差較大, 這也是為什么實用中一般不使用高次多項式插值的原因。3. 分段線性插值 分段線性插值是最簡單的分段插值，它在每個子區(qū)間xi , xi+1(i=0,1,2,n-1)上用直線段連接數(shù)據(jù)點構(gòu)成的折線函數(shù)作為插值函數(shù)的。分段線性插值函數(shù)具有計算簡單，且有每個插值子區(qū)間趨于零時，分段線性插值函數(shù)可以克服Runge現(xiàn)象，且具有隨著插值節(jié)點的增多，它可以是插值誤差不斷減小、一致逼近被插函數(shù)的特點。 3.1 分段線性插值函數(shù)的構(gòu)造過程： 在每個子區(qū)間xi , xi+1(i=0,1,

24、2,n-1)上，用對應(yīng)的兩點：x xi xi+1y yi yi+1做過此兩點的直線： L(x)=k(x -xi)+ yi這里k=(yi+1 -yi)/(xi+1 -xi)，然后將所有直線連起來即可。顯然這樣得到的折線滿足插值條件。 分段線性插值函數(shù)是一個分段函數(shù)。3.2 用分段線性插值函數(shù)做近似計算的算法：1. 輸入n個插值點： (xi, yi),i=1,n2. 輸入要近似計算的自變量點xa3.尋找包含xa的小區(qū)間xi , xi+14.用 xi , xi+1 上的線性插值函數(shù)在xa處的值作為f(xa)的函數(shù)值3.3 用分段線性插值函數(shù)做近似計算的程序：Clearx,a,b;lia_,b_,x_

25、:=(b2-a2)/(b1-a1)*(x-a1)+a2xi=Input"xi="yi=Input"yi="xa=Input"xa="n=Lengthxi;Ifxa<xi1|xa>xin,Print"超限"Break, DoIfxa>=xik,m=k,k,1,n-1; Print"m=",m,"x=",xim; lixim,yim,xim+1,yim+1,xa說明：本程序用分段線性插值做近似計算。程序執(zhí)行后，按要求通過鍵盤輸入插值基點xi: x1, . , x

26、n 和對應(yīng)函數(shù)值yi: y1 , , yn 和要做近似計算的點xa后，程序?qū)⒔o出包含xa小區(qū)間xi , xi+1的下標(biāo)i和xi ,和函數(shù)f(xa)的近似值。如果xa不在x1, xn 內(nèi)程序給出超限提示。程序中變量說明：xi：存放插值基點 x1, . , xn yi: 存放對應(yīng)函數(shù)值 y1 , , ynxa: 存放要做近似計算的自變量值m:存放包含xa小區(qū)間xi , xi+1的下標(biāo)i注：在Mathematica中有一個求分段線性插值函數(shù)的命令，命令形式為： Interpolation x1,y1,x2,y2,.,xn,yn，InterpolationOrder -> 1用如上Mathema

27、tica命令求出的分段線性插值函數(shù)沒有給出具體的分段函數(shù)表達(dá)式, 而是用“InterpolatingFunctionx1, xn, <> ”作為所求的分段插值函數(shù), 通?？梢杂?變量=Interpolation 數(shù)據(jù)點集合把所得的分段插值函數(shù)存放在變量中, 如果要計算插值函數(shù)在某一點的如p點的值,只要輸入 “變量p” 即可，如果要表示這個插值函數(shù)應(yīng)該用 “變量x” 。此外，命令 Interpolation x1,y1,x2,y2,.,xn,yn，InterpolationOrder -> 3可以給出更好的分段插值函數(shù)。 3.4例題與實驗例4.給定數(shù)據(jù)表xi 1 2 3 4 y

28、i 0 -5 -6 31） 用分段線性插值函數(shù)求函數(shù)f(x)在x =1.4的近似值2） 用Mathematica命令畫出分段線性插值圖形解：執(zhí)行求分段線性插值函數(shù)程序后，在輸入的兩個窗口中按提示分別輸入1，2，3，4、0，-5，-6，3、1.4每次輸入后用鼠標(biāo)點擊窗口的“OK”按扭，得如下插值函數(shù)：m=1 x=1-2此結(jié)果說明f(x)在x =1.4的近似值為2，它對應(yīng)的小區(qū)間為x1 , x2。為完成第二個問題，鍵入命令：d=Interpolation 1,0,2,-5,3,-6,4,3,InterpolationOrder -> 1Plotdx,x,1,4得輸出如下：Interpolat

29、ingFunction1, 4, <>圖3例4的分段線性插值函數(shù)圖形。例5.在例3中，用分段插值函數(shù)代替n次多項式插值函數(shù)，觀察插值節(jié)點增多時誤差的變化情況。 解: 對函數(shù)f(x) = (1+x2)-1 在區(qū)間-4,4內(nèi)選取不同個數(shù)的等距插值節(jié)點做觀察，這里分別選-4,4 內(nèi)的9個和21個的等距節(jié)點來做實驗，將對應(yīng)的插值函數(shù)圖與被插函數(shù)f(x) = (1+x2)-1畫在一起做觀察，為簡單起見，這里用Mathematica 命令做實驗，對應(yīng)命令為： In1:= u=Tablex,(1+x2)-1,x,-4,4 ; （*采取f(x) 在-4,4 內(nèi)的9個插值點 In2:= g=List

30、Plotu, PlotStyle->PointSize0.04 （*將散點圖圖形文件存放在變量g中In3:= s= Interpolationu,InterpolationOrder -> 1 ; （*將線性插值函數(shù)存放在變量s中In4:=t=Plotsx, (1+x2)-1, x,-4,4,PlotStyle->Thickness0.005,Thickness0.006 （*將線性插值函數(shù)s(x)與f(x)畫在一起的圖形文件存放在變量t中In5:= Showt, g （*將散點圖, 線性插值函數(shù)s(x), f(x)畫在一個坐標(biāo)系中圖 4 在-4,4中選9個等距節(jié)點的線性插值

31、函數(shù)與被插函數(shù)圖，粗線為被插函數(shù)圖In6:= u=Tablex,(1+x2)-1,x,-4,4，0.4 ; （*采取f(x) 在-4,4 內(nèi)的21插值點 In7:= g=ListPlotu, PlotStyle->PointSize0.04 In8:= s= Interpolationu,InterpolationOrder -> 1 In9:=t=Plotsx, (1+x2)-1, x,-4,4,PlotStyle->Thickness0.005,Thickness0.006 In10:= Showt, g （*將散點圖, 線性插值函數(shù)s(x), f(x)畫在一個坐標(biāo)系中圖

32、 5 在-4,4中選21個等距節(jié)點的插值函數(shù)與被插函數(shù)圖，粗線為被插函數(shù)圖 從圖形上觀察, 說明增加插值節(jié)點能保證使插值誤差變小。因此當(dāng)插值節(jié)點越多時分段線性插值誤差會越小。4樣條插值 分段插值會導(dǎo)致插值函數(shù)在子區(qū)間的端點(銜接處)不光滑,即導(dǎo)數(shù)不連續(xù),對于一些實際問題,不但要求一階導(dǎo)數(shù)連續(xù),而且要求二階導(dǎo)數(shù)連續(xù)。三次樣條插值可以滿足這些要求的概念。 所謂“樣條”(SPLINE)是工程繪圖中的一種工具,它是有彈性的細(xì)長木條,繪圖時,用細(xì)木條連接相近的幾個結(jié)點,然后再進(jìn)行拼接,連接全部結(jié)點,使之成為一條光滑曲線,且在結(jié)點處具有連續(xù)的曲率。樣條函數(shù)是對這樣的曲線進(jìn)行數(shù)學(xué)模擬得到的。它除了要求給出各

33、個結(jié)點處的函數(shù)值外,只需提供兩個邊界點處導(dǎo)數(shù)信息,便可滿足對光滑性的不同要求項式。這里給出最常用的三次樣條插值的構(gòu)造與算法。4.1 三次樣條插值函數(shù)的構(gòu)造過程： 設(shè)三次樣條插值函數(shù)為S(x)。記Mi = S(xi), f(xi)= fi= yi ,考慮它在任一區(qū)間xi,xi+1上的形式.根據(jù)三次樣條的定義可知 ,S(x)的二階導(dǎo)數(shù)S(x)在每一個子區(qū)間xi,xi+1 ( i=0,1,2,n-1)上都是線性函數(shù).于是在xi,xi+1 上S(x)=Si(x)的二階導(dǎo)數(shù)表示成其中 hi= xi+1xi . 對S(x)連續(xù)積分兩次,并利用插值條件S(xi)= yi ,得到因此，只要能求出所有的M i，

34、就能求出三次樣條插值函數(shù)S(x)。 為求出Mi，利用三次樣條插值函數(shù)S(x)在節(jié)點處的連續(xù)性S'(xi-)= S'(xi+) ,(i=1,2,n-1)可得iMi-1+2Mi+iMi+1= dI (i=1,2,n-1) 其中 上面的方程組有n-1個方程，但有n+1個變量Mi，故需兩個方程才能求唯一解，為此引入下列常用的邊界條件。第一型邊界條件： 已知f(x)在兩端點的導(dǎo)數(shù)f(a)和f(b) ，要求 S(a) = f(a) , S(b) = f(b)第二型邊界條件： 已知f(x)在兩端點的二階導(dǎo)數(shù)f"(a)和f"(b) ，要求 S(a)=M0 = f(a) ,

35、S(b)=Mn= f(b)特別當(dāng) S(a)= S(b) =0時，S(x)稱為自然三次樣條第三型邊界條件： 已知f(x)是以b -a為周期的周期函數(shù) ，要求S(x)滿足周期條件S (a) = S(b) , S(a+)= S(b-) , S(a+)= S(b-)三次樣條插值問題加上第i型邊界條件稱為第i型插值問題（i，）?？梢宰C明第i型插值問題的解是存在且唯一的。對第一型邊界條件和第二型邊界條件，它們的參數(shù) Mk 對應(yīng)如下的三對角方程組：其中，在第一型邊界條件下：l0=1, d0=6(fx0,x1-f¢(x0)/h0mn=1, dn=6(f¢(xn)-fxn-1,xn) /hn

36、-1在第二型邊界條件下：l0=0, d0=2f¢¢(x0)mn=0, dn=2f¢¢(xn)以上各組條件與方程組（*）聯(lián)立，可以解出未知參數(shù)M0，M1 ,Mn，然后代入S(x) 表達(dá)式，即可求得樣條函數(shù) 。上面構(gòu)造方法中Mi相應(yīng)于力學(xué)中細(xì)梁在xi處截面的彎矩，每一個方程中又至多出現(xiàn)相鄰的三個Mi，通常稱為三彎矩法。 4.2 求三次樣條插值函數(shù)算法：1.確定邊界條件，判定是第幾型插值問題；2.根據(jù)所確定的條件計算各值，形成方程組(*)；3.解三對角方程組(*)，求得M0， M1 ， M2, Mn ；4.將求得的Mi值代回S(x)的表達(dá)式中4.3 求三次樣條

37、插值函數(shù)程序：mar=Input"邊界條件類型=1或2(1:第一型,2：第二型)"f1=Input"輸入左邊界條件函數(shù)值="f2=Input"輸入右邊界條件函數(shù)值="xi=Input"輸入插值節(jié)點xi="yi=Input"輸入對應(yīng)函數(shù)值yi="n=Lengthxi-1;h=Tablexii+1-xii,i,1,n;mu=Tablehi/(hi+hi+1),i,1,n-1;lam=1-mu;d=Table 6(yii+1-yii)/hi-(yii-yii-1)/hi-1)/(hi-1+hi), i

38、,2,n;Ifmar=2, d=Prependd,2f1; d=Appendd,2f2; mu=AppendPrependmu,0,0; lam=AppendPrependlam,0,0, d=Prependd,6*(yi2-yi1)/h1-f1)/h1; d=Appendd,6*(f2-(yin+1-yin)/hn); mu=AppendPrependmu,0,1; lam=AppendPrependlam,1,0c=lam;b=Table2,n+1;a=mu;n=n+1;eps1=0.000000001;Do ak=ak/bk-1;bk=bk-ak*ck-1; dk=dk-ak*dk-1,

39、 k,2,n;x=Table0,n;xn=dn/bn;Doxk=(dk-ck*xk+1)/bk, k,n-1,1,-1;Print"Mk的值為 " , x/N;Clearx1,x2,y1,y2,m1,m2,h1,t,skskx1_,x2_,y1_,y2_,m1_,m2_,h1_,t_:=m1*(x2-t)3/(6*h1)+m2*(t-x1)3/(6*h1)+(y1- m1*h12/6)*(x2-t)/h1+(y2-m2*h12/6)*(t-x1)/h1;sp1=Tableskxii,xii+1,yii,yii+1,xi,xi+1,hi,t, i,1,n-1;Print&qu

40、ot;所求樣條函數(shù)S(t)為"DoPrint"S(t)=",sp1i; Print" ",xii,"<=t<",xii+1, i,1,n-1說明：本程序用于求具有第一型邊界條件和第二型邊界條件的三次樣條插值函數(shù)。程序執(zhí)行后，按要求通過鍵盤輸入：邊界條件類型=1或2(1:第一型,2：第二型)、左邊界條件函數(shù)值f¢(x0)或f¢¢(x0)、右邊界條件函數(shù)值f¢(xn)或f¢¢(xn)、插值基點xi:x0 , x1, . , xn 和對應(yīng)函數(shù)值yi: y0 ,

41、 y1 , , yn 后，程序即可給出M0， M1 ， M2, Mn和三次樣條插值函數(shù)。程序中變量說明：xi：存放插值基點 x0, x1, . , xn yi: 存放對應(yīng)函數(shù)值 y0 , y1, , ynh: 存放小區(qū)間長度hk mu: 存放數(shù)值mk lam 存放數(shù)值lk d: 存放數(shù)值dkx：存放求出的Mk值M0 , M1, . ,Mn sp1: 存放求出的三次樣條插值函數(shù)s(x)skx1,x2,y1,y2,m1,m2,h1,t:在區(qū)間x1,x2上的三次樣條插值函數(shù)s(t)注：1）語句sp1=Tableskxii,xii+1,yii,yii+1,xi,xi+1,hi,t, i,1,n-1是以

42、表的形式給出三次樣條插值函數(shù)s(t)的，其中在sp1中，第i個分量表示三次樣條插值函數(shù)s(t)的第i段.4.4例題與實驗例6. 已知函數(shù)f(x）的取值如下表： x 27.7 28 29 30f(x) 4.1 4.3 4.1 3.0 f(x) 3.0 -4.0試求f(x) 在27.7,30上的三次樣條插值函數(shù)并計算f(28.4)的近似值。解：這是第一類邊界條件的問題，執(zhí)行三次樣條插值程序后，在輸入的五個窗口中按提示分別輸入：1、3.0 、4.0、27.7, 28, 29, 30、4.1, 4.3, 4.1, 3，每次輸入后用鼠標(biāo)點擊窗口的“OK”按扭，得如下輸出結(jié)果： Mk的值為 -23.531

43、4, 0.39604, 0.829703, -9.11485所求樣條函數(shù)S(t)為S(t)=14.8432 (28 - t) - 13.073 (28 - t)3 + 14.3135 (-27.7 + t) + 0.220022 (-27.7 + t)3 27.7£t<28S(t)=4.23399 (29 - t) + 0.0660066 (29 - t)3 + 3.96172 (-28 + t) + 0.138284 (-28 + t)3 28£t<29S(t)=3.96172 (30 - t) + 0.138284 (30 - t)3 + 4.51914 (

44、-29 + t) -1.51914 (-29 + t)3 29£t<30因為28.4Î28,29接著鍵入： sp12/.t->28.4得輸出: 4.14819，因此f(28.4)» 4.14819例7. 已知函數(shù)f(x）的取值如下表： x 0 1 2 3f(x) 0 0.5 2 1.5 試求1) 給定邊界條件f¢(0)=0.2，f¢(3)=-1，求相應(yīng)的三次樣條插值函數(shù)s(x)并畫出圖形2) 給定邊界條件f²(0)=-0.3，f²(3)=3.3，求相應(yīng)的三次樣條插值函數(shù)s(x)并畫出圖形 解：問題1是第一類邊界條

45、件的問題，執(zhí)行三次樣條插值程序后，在輸入的五個窗口中按提示分別輸入：1、0.2 、1、0, 1, 2, 3、0, 0.5, 2, 1.5，每次輸入后用鼠標(biāo)點擊窗口的“OK”按扭，得如下輸出結(jié)果：Mk的值為 -0.36, 2.52, -3.72, 0.36所求樣條函數(shù)S(x)為S(t)=0.06 (1 - t) - 0.06 (1 - t)3 + 0.08 t + 0.42 t3 0£t<1S(t)=0.08 (2 - t) + 0.42 (2 - t)3 + 2.62 (-1 + t) - 0.62 (-1 + t)3 1£t<2S(t)=2.62 (3 - t) - 0.62 (3 - t)3 +