三次樣條插值的MATLAB實(shí)現(xiàn)

1、MATLAB程序設(shè)計(jì)期中考查在許多問(wèn)題中，通常根據(jù)實(shí)驗(yàn)、觀測(cè)或經(jīng)驗(yàn)得到的函數(shù)表或離散點(diǎn)上的信息， 去研究分析函數(shù)的有關(guān)特性。其中插值法是一種最基本的方法，以下給出最基本 的插值問(wèn)題一一三次樣條插值的基本提法：對(duì)插值區(qū)間a,b進(jìn)行劃分：a x0 x1xn b，函數(shù)y f x在節(jié)點(diǎn)Xi上的值 y f Xi i 0,1,2, n，并且如果函數(shù) S x在每個(gè)小區(qū)間 Xi, Xi i上 是三次多項(xiàng)式，于a,b上有二階連續(xù)導(dǎo)數(shù)，則稱Sx是a,b上的三次樣條函數(shù)， 如果S x在節(jié)點(diǎn)xi上還滿足條件S Xiyii 0,1, n則稱S x為三次樣條插值函數(shù)。三次樣條插值問(wèn)題提法：對(duì) a,b上給定的數(shù)表如下XX0

2、X1Xnyy。y1yn求一個(gè)分段三次多項(xiàng)式函數(shù)S x滿足插值條件S xi yi i 0,1, n式，并在 插值區(qū)間a,b上有二階連續(xù)導(dǎo)數(shù)。這就需要推導(dǎo)三次樣條插值公式：記f x在節(jié)點(diǎn)xi處的值為f Xj mi ( i 0,1, n )(這不是給定插值問(wèn)題數(shù) 表中的已知值)。在每個(gè)小區(qū)間Xi,Xi 1利用三次Hermite插值公式，得三次插值 公式：S x i x yi i 1 x yi 1i Xi mi 11， x Xi, Xi 1。為了得到這個(gè)公式需要4n個(gè)條件：(1).非端點(diǎn)處的界點(diǎn)有2n個(gè)；(2). 一階導(dǎo)數(shù)連續(xù)有n 1個(gè)條件；(3).二階導(dǎo) 數(shù)連續(xù)有n 1個(gè)條件，其中邊界條件：O1 S

3、 X0m。S Xn S x00S xnYnXoS xnS x00 S xn其中:Xj0,1,i xj 0i xj0 且（i,j 0,1）。i Xj0,1,mi為對(duì)應(yīng)變量的一階導(dǎo)數(shù)。其推導(dǎo)過(guò)程如下:為了確定mi的值,展開(kāi)為:Sx2X Xi 1 h 2 Xhi32xix xi hi 2 xi 1 xYi h?Yi i2X1 X Xihi22X XiXmi-廠hiXi 1mi 1,這里hi連續(xù)求兩次導(dǎo)，得:6x 2xi4Xi 1hi2 mi6x 4xi 2xi 1hi2mi6 xixi 1 2xLYi1yi。于是考慮S在節(jié)點(diǎn)Xi處的右極限值，得:S xi0 mi mihihi6hi2Yi iYi 。

4、同理,在相鄰小區(qū)間Xi 1,Xi上可得Sx的表達(dá)式為:6x 2xi 1 4xihTmi 16x 4xhi121 2xi6 xi 1 xi 2xmiYiYi i及S x在節(jié)點(diǎn)Xi處的左極限值為:2S xi0mi 1hi 1hi62 Yi1Yi。利用S x二階導(dǎo)數(shù)于節(jié)點(diǎn)Xi處的連續(xù)性條件S Xi 0S xi這里有下式成立:1c 11mi 12mihi 1hi 1hi1mi 1 hiYi 1hi2YYiYi 1hi 12'用丄-除等式兩hi 1 hi邊，并注意yi £ , _yy f xi, xi !,上式可簡(jiǎn)記為:hiEi 12mim 1gii 1,2,n 1 ,且ihii1hi

5、i1gi 3i f Xi 1, Xii f Xi , Xihiihihi 1hi最后求得m1mn的線性方程組為：2i00 0 1葉g12 220 0 0m2g2(*)0 00n 1 2n 1mn 1gn 1n 000n 2mngn通過(guò)以上復(fù)雜的求解和迭代，就可以求解出插值函數(shù)的近似表達(dá)式。 得出來(lái)的表 達(dá)式就可以用MATLAB軟件來(lái)求解。具體求解過(guò)程如下：已知n對(duì)數(shù)據(jù)點(diǎn)Xi,比，X2,y2，X3,y3，Xn，yn，假設(shè)函數(shù)關(guān)系為y f x，但解析式不確定，數(shù)據(jù)插值就是構(gòu)造函數(shù)關(guān)系式y(tǒng) g x，使Xi i 1,2,3, ,n，滿足關(guān)系 g Xif Xi。例題：求滿足下面函數(shù)表所給出的插值條件的三

6、次自然樣條函數(shù)。X1 245y f x1342分析：表中所列出的是函數(shù)對(duì)點(diǎn)，首先要把對(duì)應(yīng)的插值函數(shù)求出來(lái)，再用MATLAB軟件來(lái)求區(qū)間1,5上間隔為的各點(diǎn)的值 求解過(guò)程如下：因自然樣條插值函數(shù)的邊界條件為S Xo S Xn 0,這里n 3，故確定mo,mi,m2,m3的方程組形式形如上面的(*)式，其中系數(shù)hi :ho2 1 1,h142h12h21;i :1J2Jhoh13h1h：2 31i1112 1239gi:g131 f X1,X21 f X0,X12；7g232f X2,X32 f X1,X22；h°go3fX0,X12y。6;h2g33fX2,X32y36.將上述參數(shù)帶入

7、(*)式，得到以下方程組：和gi可按如下步驟進(jìn)行:h2235 41;210 0621m°92-033葉212m2702 -33m32001 26解得:5m2-49m38177mo,mi,84由公式2 .xxi 1hi2 x Xjhiyi2 .x xhi2 xi 1hixyi 1可知，1 3227,x 1,4XXx 1,S x8843 345 2103x 4,5XXx 33,884hi2mi2X Xi XX"-2 一 mi 1, hi2xx, 1 xxi由所求出的表達(dá)式可知區(qū)間1,5可分為1,44,5，對(duì)兩個(gè)區(qū)間分別用MATLAB命令即可: 針對(duì)第一個(gè)區(qū)間：y 8x3 8x2 r 1;其圖像如下3.5y2.51.51.52.53.5x0.5 1命令如下：x=1:4;y=(-1/8)*x43+(2/8)*x.A2+(7/4)*x-1;x