新边界条件下的三次样条插值函数

于洋++袁健华++钱江++王美玲

摘要：本文讨论了三次样条插值函数（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），这对大型矩阵计算是有益的。

正如很多实际问题需要从一维推广到二维、三维情况，我们期望进一步地工作可以在高维样条插值函数的构造上展开，将本文提出的新边界条件构造三次样条函数的方法应用于三维情况。综上，本文对三次样条函数的构造提供了新思路，为工业样条的计算提供了理论依据。