新邊界條件下的三次樣條插值函數

于洋++袁健華++錢江++王美玲

摘要：本文討論了三次樣條插值函數（cubic spline）邊界條件的更一般情形。將邊界條件的“端點”導數條件換成“中間插值節點”的導數，從而將求樣條函數的“三對角矩陣”進行了降階并向“上（下）三角矩陣”的轉化。在理論上證明了新邊界條件下樣條函數構造的唯一性，并通過數值實驗驗證了新邊界條件下的樣條函數與元函數有較好地擬合度。

關鍵詞：計算數學；三次樣條插值；邊界條件；三對角矩陣；降階

中圖分類號：0241.3

文獻標識碼：A

DOI：10.3969/j.issn.1003-6970.2016.02.006

引言

在科學研究與工程設計的諸多問題中，經常會要求計算函數值，或者給出表征某一現象的函數解析式。當遇到復雜問題時，常常找不到簡單函數顯式地表出物理量間的關系，這時根據觀測的若干離散數據點構造簡單的函數代替未知函數或復雜函數是高效的解決辦法，即插值法。插值法是函數逼近的一個重要方法，是計算數學中常用的技巧，在實際應用問題中也有著廣泛的應用。本文討論的核心內容是在處理實際問題中常用的樣條插值函數。樣條函數是計算數學以及計算機輔助幾何設計的重要工具。近年來，隨著計算機技術的不斷進步，樣條函數的應用研究得到了迅速的發展和廣泛的應用。

樣條（spline）是繪圖員用來描繪光滑曲線的一種簡單工具。為了把一些指定點（樣點）連接成一條光滑曲線，往往用一條富有彈性的細長材料，如木條（樣條）固定在樣點上，其他地方讓它自由彎曲。這樣連接樣點的曲線稱為樣條曲線。樣條插值函數就是對樣條曲線進行數學模擬得到的。下文我們主要針對三次樣條插值函數進行分析和論證。

設f（x）是區間[a，b]上的一個二次連續可微函數。在區間[a，b]上給定一組基點：

是二次連續可微的，Si（x），（i=1，2，…，n）均是不高于三次的多項式或零多項式，且滿足條件

則稱S（x）為函數f（x）的三次樣條插值函數，簡稱三次樣條。

記mi=s"（xi），f（xi）=fi，根據三次樣條的定義知，S（x）的二階導數S"（x）在每個子區間[xi，xi+1]上都是線性函數。所以在[xi，Xi+l]上S（x）=Si（x）的二階導數

對（4）式連續積分兩次，并根據條件（2）和（3）可得

hi-1mi-1+2（hi-1+hi）mi+mi+1=6{f[xi，Xi+1]-f[xi-1，xi]}（5）

其中i=2，…，n，這是一個含有n+l個未知量m1，m2，…，mn+1，而只有n-l個方程的線性方程組。通常要補加兩個約束條件，稱為邊界條件或端點條件，使（5）式成為由n+l個方程構成的方程組。

一般，為使上述方程組有唯一解，在區間[a，b]的端點a，b（即x1，xn+l）處有以下三種邊界條件：

（a.）給定f（x）在兩端點d，6的二階導數；

（b.）給定f（x）在兩端點a，b的一階導數；

（c.）未給任何導數時，求S'（x）在兩端點的近似值。以X1，x2，x3，X4為基點做一個三次Newton插值多項式Na（X）；再以點xn+1，xn，xn-1，xn-2可得Nb（x），令S'（a）=Na（a），s（b）=Nb（b），剩下的工作同（b.），這種情況也稱作Lagrange三次樣條插值函數。

1 新邊值條件及理論證明

根據上述論述知：（5）式是一個含有n+l個未知量m1，m2，…，mn+1，而只有n-l個方程的線性方程組。由矩陣知識易得，該方程組的解包含兩個自由未知量。不同于（a.）（b.）（c.）三種“端點處的邊界條件”，我們給出以下四組“關于中間插值節點的新邊界條件”：

對比新舊邊界條件易得：（A）（B.）分別是對（a.）（b.）的一般化推廣，（D.）是對（C.）的一般化推廣。對于（A.）（B.）（C.）（D.）中后k，j的關系，有如下12中情形：對（A.）邊界情況有三類：2

為了避免過多贅述，下面僅給出以上12種情況中比較有代表性的情形的理論證明：即（D.）類中3

由（4）式Si（x）求-次積分并化簡可得：

這里，且方程的系數矩陣是主對角線占優的線性方程組，即存在唯一解。則由（II）式可解mj，mj+1…，mk-2，mk-1。

接下來，問題轉化成：一方面，已知S"（xj=mj，S"（xj+1）=mj+1求關于未知量m1，m2，…，mj-1的三線型上三角方程組；另一方面，已知s"（Xk-1）=mk-1，s"（Xk）=mk求關于未知量mk+1，mk+2…，mn+1的三線型下三角方程組。

其中（12）（13）式子中Ψi，δi記法同（II）式。

綜上，即得到了所有mi（i=1，2，…，n+l）的值，從而得到唯一確定的樣條函數。

2 數值算例

下面給出新邊界條件下樣條函數的構造算法和計算實例。應用matlab編寫函數程序，得到樣條函數對原函數的擬合對比曲線。

2.1 算法

輸入：進行插值的全部節點的橫坐標xi（i=1，…，n+l），對應節點的縱坐標fi（i=l，…，n+1）和邊界條件。

輸出：三次樣條插值函數si（x）（i=1，…，n）的系數。

Step l：計算步長hi=Xi+1-xi，其中i=l，…，n；

Step 2：根據Step l的步長計算差商并形成三對角方程組（或三線型上三角方程組，三線型下三角方程組）；

Step 3：解Step 2中的方程組，找出所有mi（i=1，…，n+l）；

Step 4：將mi（i=l，…，n+l）帶回到樣條函數Si（x）（i=1，…，n）中；

Step 5：將預測點x帶入f（x）的樣條函數S（x）得近似值，matlab作圖得插值曲線。

2.2 算例

計算f（x）=X?在以下集合中的點{-5，-4，-3，-2，-1，0，1，2，3，4，5}處的三次樣條插值函數，并且滿足邊界條件f（-2）=-4，f”（2）=2。

分析 ：同上我們有記法δi=hi+hi+1，Ψi= f[x，+1，xi+2]-f[Xi，xi+1]。題目已知M4=-4，m8=2，步長hi=h=l，首先需要解如下三對角方程組：

利用追趕法懈（14）式的三對角方程組，可得m4，m5，m6，m77的值。

下面分別解m1，m2，m3對應的三線型上三角方程組：

以及解m8，m9，m10，m11對應的三線型下三角方程組：解（15）和（16）式即得剩下的所有mi的值。最后，將Si（x）整理成關于（x-xi）的三次多項式：

由（17）式中（x-Xi）?，（x-Xi）?，（X-xi）的系數及常數項關于mi，hi，fi的表達式，可計算出S（x）的系數矩陣，帶回（17）式即可得插值函數的表達式。

最后，通過matlab程序得圖1：新邊界條件下樣條函數與原函數的擬合對比圖，從圖1可見新邊界條件下的插值虛線與原函數圖形的走向一致，曲線在坐標系下的位置也相近。

與圖1相對應地，通過matlab軟件得表1：新邊界條件下樣條函數與原函數的擬合對比數據，由表1中數據可知：當步長變小時，相對誤差變小。由此可見新邊界條件下構造樣條函數的可行性和相對可靠性。

3 結論

本文給出了三次樣條插值函數更一般化的構造條件，并且由上述新邊界條件的論證過程及數值實驗可知，任給樣條插值節點的一階或者二階導數的值中的兩個，便可以將整個樣條函數確定。特別地，由（B.）情況可知，可以用“中間節點的導數值”取代“邊界節點的單側導數值”。進一步看，（c.）情況中選擇Newton插值節點時可將“端點”換成“中間節點”，即插值區間的中部，從而降低了插值誤差，提高了精度。精度的提高，會直接影響算法對生產實際問題處理的有效性。工程中模擬機翼下輪廓曲線等的插值問題可能會要求中間某點的曲率是固定值，或二階導數為固定值的情況，則本文提出的新邊界條件能有效處理這一情形。此外，邊界條件取中間節點相對于取端點的情形，對線性方程組的系數矩陣“三對角矩陣”進行了維數上的降階和形式上向“三線型上三角”（或“三線型下三角”）的轉化。同階的“三對角矩陣”的計算量是O（n?），而“三線型上三角”（或“三線型下三角”）的計算量是O（n），這對大型矩陣計算是有益的。

正如很多實際問題需要從一維推廣到二維、三維情況，我們期望進一步地工作可以在高維樣條插值函數的構造上展開，將本文提出的新邊界條件構造三次樣條函數的方法應用于三維情況。綜上，本文對三次樣條函數的構造提供了新思路，為工業樣條的計算提供了理論依據。