三次樣條插值算法詳解課件

1樣條函數(shù)的定義定義4.1

設(shè)區(qū)間[a,b]上給定一個(gè)節(jié)點(diǎn)劃分a=x0<x1<……<xn-1<xn=b 如果存在正整數(shù)k使得[a,b]上的分段函數(shù)s(x)滿(mǎn)足如下兩條：(1)在[a,b]上有直到k-1階連續(xù)導(dǎo)數(shù)。(2)在每個(gè)小區(qū)間[xi,xi+1]上是次數(shù)不大于k的多項(xiàng)式。則稱(chēng)分段函數(shù)s(x)是以(2.6)為節(jié)點(diǎn)集的k次樣條函數(shù)。1樣條函數(shù)的定義定義4.1設(shè)區(qū)間[a,b]上給定一個(gè)節(jié)點(diǎn)劃2三次樣條插值函數(shù)的定義并且關(guān)于這個(gè)節(jié)點(diǎn)集的三次樣條函數(shù)s(x)滿(mǎn)足插值條件：則稱(chēng)這個(gè)三次樣條函數(shù)s(x)為三次樣條插值函數(shù)。2三次樣條插值函數(shù)的定義并且關(guān)于這個(gè)節(jié)點(diǎn)集的三次樣條函數(shù)s(3三次樣條插值函數(shù)的邊界條件插值條件：連續(xù)性條件：一階導(dǎo)數(shù)連續(xù)條件：二階導(dǎo)數(shù)連續(xù)條件：3三次樣條插值函數(shù)的邊界條件插值條件：連續(xù)性條件：一階導(dǎo)數(shù)連4（1）因?yàn)閟(x)在每個(gè)小區(qū)間上是一個(gè)次小于三次的多項(xiàng)式，故有四個(gè)未知系數(shù)；（2）因?yàn)閟(x)有n分段，從而共有4n個(gè)未知系數(shù)?。?）但插值條件與樣條條件僅給出4n-2個(gè)條件，無(wú)法定出4n個(gè)未知系數(shù)，還差2個(gè)條件！這2個(gè)條件我們用邊界條件給出！

4（1）因?yàn)閟(x)在每個(gè)小區(qū)間上是一個(gè)次小于三次的多項(xiàng)式，5通常我們對(duì)插值多項(xiàng)式在兩端點(diǎn)的狀態(tài)加以要求也就是所謂的邊界條件:第一邊界條件：由區(qū)間端點(diǎn)處的一階導(dǎo)數(shù)給出即 5通常我們對(duì)插值多項(xiàng)式在兩端點(diǎn)的狀態(tài)加以要求也就是第一邊界條6第二邊界條件：由區(qū)間端點(diǎn)處的二階導(dǎo)數(shù)給出即

特殊情況為自然邊界條件：由區(qū)間端點(diǎn)處的二階導(dǎo)數(shù)恒為0給出即

6第二邊界條件：由區(qū)間端點(diǎn)處的二階導(dǎo)數(shù)給出即 特殊情況為自7這樣三次樣條插值問(wèn)題就分成三類(lèi)！其實(shí)不止這三類(lèi)!第三類(lèi)又稱(chēng)周期邊界條件：由區(qū)間端點(diǎn)處的函數(shù)值或?qū)?shù)值滿(mǎn)足周期條件給出

7這樣三次樣條插值問(wèn)題就分成三類(lèi)！其實(shí)不止這三類(lèi)!第三類(lèi)又稱(chēng)8樣條函數(shù)的例子容易驗(yàn)證：是滿(mǎn)足如下數(shù)據(jù)的第一類(lèi)邊界樣條插值問(wèn)題解：x0123y0000y’108樣條函數(shù)的例子容易驗(yàn)證：是滿(mǎn)足如下數(shù)據(jù)的第一類(lèi)邊界樣條插值9樣條函數(shù)的例子9樣條函數(shù)的例子10通常有三轉(zhuǎn)角法、三彎矩法、B樣條基函數(shù)法。三次樣條插值函數(shù)的求法這三種方法的基本思想是類(lèi)似的，都是通過(guò)待定某些參數(shù)來(lái)確定插值函數(shù)，但肯定不是待定4n個(gè)參數(shù)。而是利用已知條件將待定參數(shù)減小到最少。比如：待定一階導(dǎo)數(shù)、待定二階導(dǎo)數(shù)、采用基函數(shù)方法來(lái)確定插值函數(shù)。10通常有三轉(zhuǎn)角法、三彎矩法、B樣條基函數(shù)法。三次樣條插值函11三轉(zhuǎn)角法:待定一階數(shù)為了確定三次樣條插值函數(shù)的表達(dá)式S(x)，我們采用待定系數(shù)法來(lái)求解，我們待定什么系數(shù)呢？考慮到帶一階導(dǎo)數(shù)的分段三次Hermite插值多項(xiàng)式11三轉(zhuǎn)角法:待定一階數(shù)為了確定三次樣條插值函數(shù)的表達(dá)式S12我們采用待定一階導(dǎo)數(shù)的方法即設(shè)因?yàn)榉侄稳蜨ermite插值多項(xiàng)式已經(jīng)至少是一階連續(xù)可導(dǎo)了，為了讓它成為三次樣條函數(shù)只需確定節(jié)點(diǎn)處的一階導(dǎo)數(shù)使這些節(jié)點(diǎn)處的二階導(dǎo)數(shù)連續(xù)即可！12我們采用待定一階導(dǎo)數(shù)的方法即設(shè)因?yàn)榉侄稳蜨ermite131314由于在內(nèi)部節(jié)點(diǎn)處二階導(dǎo)數(shù)連續(xù)條件：整理化簡(jiǎn)后得：14由于在內(nèi)部節(jié)點(diǎn)處二階導(dǎo)數(shù)連續(xù)條件：整理化簡(jiǎn)后得：15稱(chēng)為三轉(zhuǎn)角法基本方程組以上推導(dǎo)還沒(méi)有考慮邊界條件！針對(duì)不同類(lèi)型的三次樣條問(wèn)題，就可以導(dǎo)出不同的方程組！15稱(chēng)為三轉(zhuǎn)角法基本方程組以上推導(dǎo)還沒(méi)有考慮邊界條件！針對(duì)不16第一類(lèi)三次樣條插值問(wèn)題方程組基本方程組化為n-1階方程組由于已知：化為矩陣形式16第一類(lèi)三次樣條插值問(wèn)題方程組基本方程組化為n-1階方程組17這是一個(gè)嚴(yán)格對(duì)角占優(yōu)的三對(duì)角方程組，用追趕法可以求解！17這是一個(gè)嚴(yán)格對(duì)角占優(yōu)的三對(duì)角方程組，18第二類(lèi)三次樣條插值問(wèn)題的方程組由于已知：故得：18第二類(lèi)三次樣條插值問(wèn)題的方程組由于已知：故得：19稍加整理得聯(lián)合基本方程組得一個(gè)n+1階三對(duì)角方程組，化成矩陣形式為：仍然是嚴(yán)格對(duì)角占優(yōu)19稍加整理得聯(lián)合基本方程組得一個(gè)n+1階三對(duì)角方程組，20第三類(lèi)樣條插值問(wèn)題的方程組立即可得下式：由于：20第三類(lèi)樣條插值問(wèn)題的方程組立即可得下式：由于：21其中：聯(lián)合基本方程得一個(gè)廣義三對(duì)角或周期三對(duì)角方程組：這個(gè)方程組的系數(shù)矩陣仍然是嚴(yán)格對(duì)角占優(yōu)陣！21其中：聯(lián)合基本方程得一個(gè)廣義三對(duì)角或周期三對(duì)角方程組：這22求解這些不同類(lèi)型的樣條插值問(wèn)題的方程組，我們可得所要待定的一階導(dǎo)數(shù)：稱(chēng)為三次樣條插值問(wèn)題三轉(zhuǎn)角公式！再代入S(x)的每一段表達(dá)式，就求得三次樣條函數(shù)的表達(dá)式!22求解這些不同類(lèi)型的樣條插值問(wèn)題的方程組，我們可得稱(chēng)為三次23例1.對(duì)于給定的節(jié)點(diǎn)及函數(shù)值解:這是自然邊界條件下的樣條問(wèn)題。23例1.對(duì)于給定的節(jié)點(diǎn)及函數(shù)值解:這是自然邊界條件下的樣24我們可以將上述計(jì)算列于表中：k0123xk1245yk1342mk????Mk0??0hk121*λk*2/31/3*μk*1/32/3*gk69/2-7/2-624我們可以將上述計(jì)算列于表中：k0123xk1245yk125由些得如下方程組：利用三轉(zhuǎn)角公式：25由些得如下方程組：利用三轉(zhuǎn)角公式：26同樣可以求得第三段表達(dá)式！26同樣可以求得第三段表達(dá)式！2727282829三彎矩法:待定二階導(dǎo)數(shù)選擇二階導(dǎo)數(shù)作為待定參數(shù):由于三次樣條S(x)是三次多項(xiàng)式,故它的二階導(dǎo)數(shù)是一次多項(xiàng)式,從而思考:(1)的原因?29三彎矩法:待定二階導(dǎo)數(shù)選擇二階導(dǎo)數(shù)作為待定參數(shù):由于三次303031從而推導(dǎo)出了三次樣條S(x)在第k個(gè)小區(qū)間[xk,xk+1]上的表達(dá)式為:它的系數(shù)都是用二階導(dǎo)數(shù)與函數(shù)值表示!31從而推導(dǎo)出了三次樣條S(x)在第k個(gè)小區(qū)間[xk,xk+32對(duì)所有中間節(jié)點(diǎn)xk,k=1,2,…,n-1,左邊小區(qū)間與右邊小區(qū)間上的三次多項(xiàng)式的一階導(dǎo)數(shù)應(yīng)當(dāng)連續(xù)!確定二階導(dǎo)數(shù)32對(duì)所有中間節(jié)點(diǎn)xk,k=1,2,…,n-1,左邊小區(qū)間與33三彎矩法基本方程注意到這個(gè)基本方程只包括了n-1個(gè)方程!但卻有n個(gè)二階導(dǎo)數(shù)需要待定,這是一個(gè)欠定方程組,還需要根據(jù)邊界條件再確定兩個(gè)方程!33三彎矩法基本方程注意到這個(gè)基本方程只包括了n-1個(gè)方程!34曲率調(diào)整樣條這種樣條的邊界條件是已知兩端點(diǎn)的二階導(dǎo)數(shù)值!這樣從三彎矩基本方程可以導(dǎo)數(shù)確定其它n-2個(gè)待定參數(shù)的方程組:34曲率調(diào)整樣條這種樣條的邊界條件是已知兩端點(diǎn)的二階導(dǎo)數(shù)值!35自然樣條這種樣條的邊界條件是:已知兩端點(diǎn)的二階導(dǎo)數(shù)值為0!這樣從三彎矩基本方程可以導(dǎo)數(shù)確定其它n-2個(gè)待定參數(shù)的方程組:35自然樣條這種樣條的邊界條件是:已知兩端點(diǎn)的二階導(dǎo)數(shù)值為036固支樣條這種樣條的邊界條件是:已知兩端點(diǎn)的一階導(dǎo)數(shù)值!根據(jù)前面推導(dǎo)過(guò)程中得到的樣條函數(shù)S(x)的一階導(dǎo)數(shù)的表達(dá)式(2.11),得方程36固支樣條這種樣條的邊界條件是:已知兩端點(diǎn)的一階導(dǎo)數(shù)值!根37固支樣條這樣從三彎矩基本方程可以導(dǎo)數(shù)確定n個(gè)待定參數(shù)的方程組:37固支樣條這樣從三彎矩基本方程可以導(dǎo)數(shù)確定n個(gè)待定參數(shù)的方38非扭結(jié)樣條這種樣條的邊界條件是:要求樣條S(x)在開(kāi)始的兩個(gè)小區(qū)間[x0,x1],[x1,x2]上的三階導(dǎo)數(shù)相同,在最后兩個(gè)小區(qū)間[xn-2,xn-1],[xn-1,xn]上的三階導(dǎo)數(shù)相同.對(duì)表達(dá)式(2.9)再求一次導(dǎo)數(shù)得方程38非扭結(jié)樣條這種樣條的邊界條件是:要求樣條S(x)在開(kāi)始的39非扭結(jié)樣條再由三彎矩基本方程,可得39非扭結(jié)樣條再由三彎矩基本方程,可得40周期樣條這種樣條的邊界條件是:要求樣條S(x)及其導(dǎo)數(shù)是以區(qū)間長(zhǎng)度xn-x0為周期的函數(shù)即這些條件可以確定如下兩個(gè)方程:40周期樣條這種樣條的邊界條件是:要求樣條S(x)及其導(dǎo)數(shù)是41再由三彎矩基本方程,可得周期樣條41再由三彎矩基本方程,可得周期樣條42在Matlab中數(shù)據(jù)點(diǎn)稱(chēng)之為斷點(diǎn)。如果三次樣條插值沒(méi)有邊界條件，最常用的方法，就是采用非扭結(jié)（not-a-knot）條件。這個(gè)條件強(qiáng)迫第1個(gè)和第2個(gè)三次多項(xiàng)式的三階導(dǎo)數(shù)相等，對(duì)最后一個(gè)和倒數(shù)第2個(gè)三次多項(xiàng)式也做同樣地處理。Matlab中三次樣條插值也有現(xiàn)成的函數(shù)：

y=interp1(x0,y0,x,'spline')；

y=spline(x0,y0,x)；

pp=csape(x0,y0,conds), pp=csape(x0,y0,conds,valconds)，

y=ppval(pp,x)。 其中x0,y0是已知數(shù)據(jù)點(diǎn)，x是插值點(diǎn)，y是插值點(diǎn)的函數(shù)值。對(duì)于三次樣條插值，我們提倡使用函數(shù)csape，csape的返回值是pp形式，要求插值點(diǎn)的近似函數(shù)值，必須調(diào)用函數(shù)ppval。MATLAB中三次樣條函數(shù)法實(shí)現(xiàn)42在Matlab中數(shù)據(jù)點(diǎn)稱(chēng)之為斷點(diǎn)。如果三次樣條插值沒(méi)有邊43pp=csape(x0,y0,conds,valconds) conds指定插值的邊界條件，其值可為：

'complete' 邊界為一階導(dǎo)數(shù)，一階導(dǎo)數(shù)的值在valconds參數(shù)中給出。

'not-a-knot' 非扭結(jié)條件

'periodic' 周期條件

'second' 邊界為二階導(dǎo)數(shù)，二階導(dǎo)數(shù)的值在valconds參數(shù)中給出，若忽略valconds參數(shù)，二階導(dǎo)數(shù)的缺省值為[0,0]。MATLAB中三次樣條函數(shù)法實(shí)現(xiàn)43pp=csape(x0,y0,conds,valcond44例2：第一邊界條件的例題x1245y1342y’17/8-19/8clear;x=[1,2,4,5];y=[1,3,4,2];pp=csape(x,y,'complete',[17/8,-19/8]);pp.coefsMATLAB代碼44例2：第一邊界條件的例題x1245y1342y’17/845第一邊界條件的例題pp=form:'pp'breaks:[1245]coefs:[3x4double]pieces:3order:4dim:1pp.coefs-0.12500 2.12501.0000-0.1250-0.37501.75003.00000.3750-1.1250-1.25004.0000分段多項(xiàng)式結(jié)構(gòu)與系數(shù)矩陣45第一邊界條件的例題pp=分段多項(xiàng)式結(jié)構(gòu)與系數(shù)矩陣46如下代碼求解上述樣條問(wèn)題：x-4-3-2-101234y00.151.122.362.361.460.490.060y’00例3：自然邊界條件例題46如下代碼求解上述樣條問(wèn)題：x-4-3-2-101234y47x=[-4 -3 -2 -1 0 1 2 3 4];y=[0 0.15 1.12 2.36 2.36 1.46 0.49 0.06 0];pp=csape(x,y,'second');xx=-4:0.01:4;yy=ppval(pp,xx);holdon;plot(x,y,'ok');plot(xx,yy,'k-');holdoff;MATLAB程序47x=[-4 -3 -2 -1 0 1 2 3 4];MA48pp=form:'pp'breaks:[-4-3-2-101234]coefs:[8x4double]pieces:8order:4dim:1分段多項(xiàng)式結(jié)構(gòu)：48pp=分段多項(xiàng)式結(jié)構(gòu)：49

pp.coefs=0.180856038291610.00000000000000-0.030856038291610-0.084280191458030.542568114874820.511712076583210.15000000000000-0.393735272459500.289727540500741.344007731958761.120000000000000.14922128129602-0.891478276877760.742256995581742.360000000000000.13685014727540-0.44381443298969-0