2015年4月7日 星期二

三次樣條插值 Cubic spline interpolation

沒有留言:
 

「內插」就是找一個函數,完全符合手邊的一堆數據。此函數稱作「內插函數」。換句話說,找到一個函數,穿過所有給定的函數值。外觀就像是在相鄰的函數值之間,插滿函數值,因而得名「內插」。


ㄧ、樣條插值定義
樣條插值(spline interpolation)使用分段的多項式進行插值,樣條插值可以使用低階多項式樣條實現較小的插值誤差,避免使用高階多項式所出現的龍格現象(誤差)。

假設有n+1個不同的節點xi

x_0 < x_1 < ... < x_{n-1} < x_n, \,\!

以及n+1個節點值yi,我們要找到一個n階樣條函數


S(x) := \left\{\begin{matrix}
    S_0(x) & x \in [x_0, x_1] \\
    S_1(x) & x \in [x_1, x_2] \\
    \vdots & \vdots \\
S_{n-1}(x) & x \in [x_{n-1}, x_n]
\end{matrix}\right.

其中每個Si(x)都是一個k次的多項式。



二、線性插值(Linear interpolation)
線性樣條插值是最簡單的樣條插值。數據點使用直線進行連接,如果使用線性插值,則插值函數的一階導數在表列點上會沒有定義


整個公式會像Lagrange polynomial interpolation的形式。




三、三次樣條插值(Cubic spline interpolation)
形式上是線性插值的改善,而最高次數為三次而不會過度轉彎發生抖動(龍格現象),處處圓滑跨過表列點時也一樣。

對於n+1個給定點的數據集{xi},我們可以用n段三次多項式在數據點之間構建一個三次樣條。表示對函數f進行插值的樣條函數,需要符合
  • 插值特性(曲線必通過座標點)
    • S(xi) = f(xi)
  • 樣條相互連接且兩次連續可導(曲線的銜接點必須平滑)
    • Si-1(xi)   = Si(xi)  , i=1,...,n-1
    • S'i-1(xi)  = S'i(xi) , i=1,...,n-1
    • S''i-1(xi) = S''i(xi), i=1,...,n-1

也就是說,預期插值函數的二階導數函數 f''(x) 乃由各表列值的二階導數值所線性內插構成。


因此,從推廣線性插值出發,除了表列的fi值, 我們也需要表列fi值的二次微分f''的值。希望加上一些高階項,讓函數的線條變得更圓滑,但這樣做的同時,又希望保留原本己經在線性插值式中滿足之表列點代入會得表列值的特性,我們因此要求後面加入的東西至少在表列點處其值要為零。
(用 A3 - A = A(A+1)(A-1) 這樣的形式的項,的確可以保證 C、D 在 x = xk+1 及 x = xk 都會等於零)


微分兩次,果真得到其二階導數函數恰呈線性插值的形式(倒果為因的證明,比較漂亮)。


透過上列 A、B、C、D 的定義,有了可以用的插值公式,但是這個公式還不能用,理由是公式中的 f''k 我們並不知道。現在,要用以下的策略把它們定出來:

我們看cubic spline函數的一階導數,它的使用範圍是點 xk 與點 xk+1 之間。xk+1 是這個區間的右端點,自然有在其描述的範圍之內,然而 xk+1 同時也是 xk+1 到 xk+2 的另一條曲線段的左端點,有它自已不一樣的斜率函數公式,為了確保一階導數的圓滑度,我們可以要求在 xk+1 這個點上兩邊的斜率函數值要一樣。 這個條件用下去所寫出來的一條新的等式,把所有未知的 f''k 通通移到等號左邊,則整理得以下的形式。注意所未知數形成一個線性代數方程式:


我們因此可以透過數值方法求解。一旦解出各 y''i ,則 cubic spline 插值公式的建立就確定並完成了。



[未知數個數分析]

Si(X)=aiX^3+biX^2+ciX+di,i=1~N-1。

每兩個相鄰點就有一個多項式,每個多項式有a、b、c、d4個未知數,共有N-1道多項式,所以共4N-4個未知數。我們只要把這4N-4個未知數求出來即大功告成。

1.  曲線必通過座標點:Si(Xi)=Yi,Si(Xi+1)=Yi+1,i=1~N-1。
共2N-2個式子。

2.  曲線的銜接點必須平滑,及一階導函數和二階導函數值相等。
S'i-1(Xi)=S'i(Xi),i=2~N-1。
S''i-1(Xi)=S''i(Xi),i=2~N-1。
共2N-4個式子。

3.  兩端點的二階導函數給定初值
給定初值的方法很多,我們採用natural spline的方式
令S''1(X1)=0,S''N-1(XN)=0。
共2個式子。

由以上1、2、3我們可得到4N-4個式子,即可求出全部的4N-4個變數。



[三對角矩陣 tri-diagonal matrix
上面的線性代數方程式事實上是一個簡化的特殊形式,寫成矩陣式會看到只有對角線與其近鄰的上下兩條元素不是零,其他都是。這種矩陣叫三對角矩陣 (tri-diagonal matrix),有專用的快速演算法可以解出。



[端點條件]
由i的取值範圍可知,共有n-1個公式, 但卻有n+1個未知量m 。要想求解該方程組,還需另外兩個式子。所以需要對兩端點x0和xn的微分加些限制。 選擇不是唯一的,3種比較常用的限制如下。

a. 自由邊界(Natural)
首尾兩端沒有受到任何讓它們彎曲的力。

b. 固定邊界(Clamped)
首尾兩端點的微分值是被指定的,這裡分別定為A和B。

c. 非節點邊界(Not-A-Knot)
指定樣條曲線的三次微分匹配

下圖可以看出不同的端點邊界對樣條曲線的影響:


Reference

wiki - 樣條插值

插值法(內插與外插)Interpolation and Extrapolation

Spline Interpolation 的數學概念

Chapter 10. 數值分析IV─差值法及曲線近似(2)

演算法筆記 - Interpolation

三次樣條插值(Cubic Spline Interpolation)及代碼實現(C語言)



沒有留言:

張貼留言

技術提供:Blogger.