三次样条插值函数

4.5三次样条插值函数4.5.1三次样条差值函数的力学背景在工程和数学应用中常有这么一类数据处理问题：在平面上给定了一组有序的离散点列，要求用一条光滑的曲线把这些点按次序连接起来。在过去很长一段时间内，工程技术人员为了得到这条光滑的曲线，常常使用一条富有弹性的均匀细木条（或是有机玻璃条），一次经过这些点，并用亚铁在若干点处压住，然后沿这条细木条画出一条光滑的曲线，并形象地称之为“样条曲线”。在力学上，如果把细木条看成弹性细梁，亚铁看成作用在梁上的集中载荷，“样条曲线”就可模拟为弹性细梁在外加集中载荷作用下的弯曲变形曲线。如果用A表示细梁的刚度系数，M表示弯矩，在建立坐标系后，由于样条是均匀细木条，在两个相邻压铁之间无任何外力，所以M是的线性函数，A为常数，由力学知识可得𝐴𝑘𝑥=𝑀(𝑥)其中为“样条曲线”y=y(x)的曲率，由数学知识对于一条光滑曲线kx=y′′/(1+y’2)3/2。一般来讲，上述样条曲线所适合的微分方程是非线性的，它的解是无法用初等函数表示的，但在通常称为“小挠度”的情况下，即细梁弯曲不大，𝑦′≪1时，可以忽略y′的影响，从而得到近似方程Ay′′=M(𝑥)。由于M为线性函数，故有y4(𝑥)=0，即“样条曲线”为分段多项式，且曲线的函数值、一阶导数、二阶导数都是连续的，而三阶导数是间断的，这就是三次样条差值函数的力学背景。4.5.2三次样条差值函数定义4.5.1设在区间𝑎,𝑏上给定一个分割𝑎=𝑥0𝑥1⋯𝑥𝑛−1𝑥𝑛=𝑏如果定义在𝑎,𝑏上的一个函数S(𝑥)满足下列条件：(1)在每个小区间𝑥𝑛−1,𝑥𝑛𝑖=1,2,⋯,𝑛上，𝑆(𝑥)是三次多项式；(2)在整个区间𝑎,𝑏上，𝑆(𝑥)为二阶连续可导函数，也就是说在每个节点𝑥𝑖𝑖=1,2,⋯,𝑛−1处𝑆(𝑘)𝑥𝑖−0=𝑆(𝑘)𝑥𝑖+0,𝑘=0,1,2(4.5.2)则称𝑆(𝑥)为三次样条插值函数。对定义在区间𝑎,𝑏上的函数𝑓𝑥，如果存在三次样条插值函数𝑆(𝑥)，使得在节点处满足𝑆𝑥𝑖=𝑓𝑥𝑖,𝑖=0,1,2,⋯,𝑛(4.5.3)就称𝑆(𝑥)为插值于𝑓𝑥的三次样条函数，简称三次样条插值函数。问题一：如果并不知道函数𝑓𝑥的解析表达式，而只知道其在节点处的值𝑦𝑖=𝑓𝑥𝑖𝑖=0,1,⋯,𝑛，如何估计𝑓𝑥？一个很自然的方法就是求插值于𝑦𝑖三次样条插值函数𝑆(𝑥)，以𝑆(𝑥)作为对𝑓𝑥的逼近。问题二：那么如何求𝑆(𝑥)呢？可利用𝑦𝑖及其一阶、二阶导数来建立𝑆(𝑥)的表达式及连续方程。M连续方程与𝑆(𝑥)的表达式记𝑀𝑖=𝑆′′𝑥𝑖，𝑓𝑥𝑖=𝑦𝑖，根据三次样条插值函数的定义可知𝑆(𝑥)的二阶导数𝑆′′(𝑥)在每一个子区间𝑥𝑖−1,𝑥𝑖上都是线性函数，于是在𝑥𝑖−1,𝑥𝑖上𝑆(𝑥)的二阶导数可表示为𝑆′′𝑥=𝑀𝑖−1𝑥𝑖−𝑥ℎ𝑖+𝑀𝑖𝑥−𝑥𝑖−1ℎ𝑖,𝑥∈𝑥𝑖−1,𝑥𝑖，(4.5.4)对𝑆′′𝑥连续积分两次得𝑆′𝑥=−𝑀𝑖−1𝑥𝑖−𝑥22ℎ𝑖+𝑀𝑖𝑥−𝑥𝑖−122ℎ𝑖+A𝑖,𝑥∈𝑥𝑖−1,𝑥𝑖，(4.5.5)𝑆𝑥=𝑀𝑖−1𝑥𝑖−𝑥36ℎ𝑖+𝑀𝑖𝑥−𝑥𝑖−136ℎ𝑖+A𝑖𝑥+B𝑖,𝑥∈𝑥𝑖−1,𝑥𝑖，(4.5.6)由插值条件(4.5.5)和(4.5.6)得𝐴𝑖=𝑦𝑖−𝑦𝑖−1ℎ𝑖−ℎ𝑖𝑀𝑖−𝑀𝑖−16,𝐵=𝑦𝑖−1ℎ𝑖−ℎ𝑖𝑀𝑖−16𝑥𝑖+−𝑦𝑖ℎ𝑖+ℎ𝑖𝑀𝑖6𝑥𝑖−1,由式(4.5.5)得𝑆′𝑥𝑖−0=𝑦𝑖−𝑦𝑖−1ℎ𝑖−ℎ𝑖(𝑀𝑖−𝑀𝑖−1)6+𝑀𝑖ℎ𝑖2𝑆′𝑥𝑖+0=𝑦𝑖+1−𝑦𝑖ℎ𝑖+1−ℎ𝑖+1𝑀𝑖+1−𝑀𝑖6−𝑀𝑖ℎ𝑖+12由于一阶导数连续，S′𝑥𝑖−0=𝑆′𝑥𝑖+0(𝑖=1,…,𝑛−1)，所以由上式可得n-1个等式。𝜇𝑖𝑀𝑖−1+2𝑀𝑖+𝜆𝑖𝑀𝑖+1=𝑑𝑖，𝑖=1,…,𝑛−1该方程组称为S(x)的M连续方程。𝜆𝑖=ℎ𝑖+1ℎ𝑖+ℎ𝑖+1𝜇𝑖=ℎ𝑖ℎ𝑖+ℎ𝑖+1𝑑𝑖=6ℎ𝑖+ℎ𝑖+1𝑦𝑖+1−𝑦𝑖ℎ𝑖+1−𝑦𝑖−𝑦𝑖−1ℎ𝑖=6𝑓[𝑥𝑖−1,𝑥𝑖,𝑥𝑖+1]𝜇𝑖𝑀𝑖−1+2𝑀𝑖+𝜆𝑖𝑀𝑖+1=𝑑𝑖，𝑖=1,…,𝑛−1𝜆𝑖=ℎ𝑖+1ℎ𝑖+ℎ𝑖+1𝜇𝑖=ℎ𝑖ℎ𝑖+ℎ𝑖+1𝑑𝑖=6ℎ𝑖+ℎ𝑖+1𝑦𝑖+1−𝑦𝑖ℎ𝑖+1−𝑦𝑖−𝑦𝑖−1ℎ𝑖=6𝑓[𝑥𝑖−1,𝑥𝑖,𝑥𝑖+1]相邻区间的长度比插值数据在处的二阶中心差商的3倍𝑥𝑖故式子说明插值函数的二阶导数在𝑥𝑖−1,𝑥𝑖,𝑥𝑖+1三点处的加权平均（权因子分别为𝜇𝑖/3，2/3，𝜆𝑖/3）为被插数据在𝑥𝑖处的二阶中心差商。这就是力学上的“三弯矩”关系。m连续方程与𝑆(𝑥)的表达式利用Hermite插值公式，在区间𝑥𝑖,𝑥𝑖+1上𝑆𝑥=𝑚𝑖(𝑥𝑖+1−𝑥)2(𝑥−𝑥𝑖)ℎi+12−𝑚𝑖+1𝑥−𝑥𝑖2(𝑥𝑖+1−𝑥)ℎi+12+𝑦𝑖(𝑥𝑖+−𝑥)2[2𝑥−𝑥𝑖+ℎ𝑖+1]ℎ𝑖+13+𝑦𝑖+1(𝑥−𝑥𝑖)2[2𝑥𝑖+1−𝑥+ℎ𝑖+1]ℎ𝑖+13连续求导两次得𝑆′′𝑥=2𝑚𝑖3𝑥−2𝑥𝑖+1−𝑥𝑖ℎ𝑖+12+2𝑚𝑖+13𝑥−𝑥𝑖+1−2𝑥𝑖ℎ𝑖+12+6𝑦𝑖+1−𝑦𝑖ℎ𝑖+13(𝑥𝑖+1+𝑥𝑖−2𝑥)𝑆′′𝑥𝑖−0=2𝑚𝑖−1ℎ𝑖+4𝑚𝑖ℎ𝑖−6𝑦𝑖−𝑦𝑖−1ℎ𝑖2𝑆′′𝑥𝑖+0=−4𝑚𝑖ℎ𝑖+1−2𝑚𝑖+1ℎ𝑖+1+6𝑦𝑖+1−𝑦𝑖ℎ𝑖+12从而得由于二阶导数连续，S′′𝑥𝑖−0=𝑆′′𝑥𝑖+0𝑖=1,…,𝑛−1，式子可化为𝜆𝑖𝑚𝑖−1+2𝑚𝑖+𝜇𝑖𝑚𝑖+1=𝑐𝑖，𝑖=1,…,𝑛−1该式子为S(x)的m连续方程其中𝜆𝑖，𝜇𝑖同上𝑐𝑖=3𝜆𝑖𝑦𝑖−𝑦𝑖−1ℎ𝑖+𝜇𝑖𝑦𝑖+1−𝑦𝑖ℎ𝑖+1，𝑖=1,…,𝑛−1边界条件方程组都为n+1个未知数、n-1个方程的线性方程组无法保证唯一解加两个条件按具体问题的要求在区间端点给出约束条件，称为边界条件边界条件的分类第一类边界条件给定区间[𝑎,𝑏]两端点的斜率𝑚0,𝑚𝑛𝑆′𝑥0=𝑦0′=𝑚0，𝑆′𝑥𝑛=𝑦𝑛′=𝑚𝑛第二类边界条件给定区间[𝑎,𝑏]两端点的二阶导数𝑀0,𝑀𝑛𝑆′′𝑥0=𝑦0′′=𝑀0，𝑆′′𝑥𝑛=𝑦𝑛′′=𝑀𝑛第三类边界条件假设𝑦=𝑓(𝑥)是以𝑏−𝑎为周期的周期函数，则要求三次样条插值函数𝑆(𝑥)也是周期函数，对𝑆(𝑥)加上周期条件𝑆(𝑝)𝑥0+0=𝑆(𝑝)𝑥𝑛−0，𝑝=0,1,2讨论M连续方程的各类边界条件对于第一类边界条件，可利用M连续方程由𝑆′𝑥=−𝑀𝑖−1𝑥𝑖−𝑥22ℎ𝑖+𝑀𝑖𝑥−𝑥𝑖−122ℎ𝑖+A𝑖,𝑥∈𝑥𝑖−1,𝑥𝑖得到边界上的关系式𝑚0=−𝑀0ℎ12+𝐴1𝑚𝑛=𝑀𝑛ℎ𝑛2+𝐴𝑛带入A的表达式可得2𝑀0+𝑀1=6ℎ1(𝑦1−𝑦0ℎ1−𝑚0)𝑀𝑛−1+2𝑀𝑛=6ℎ𝑛(𝑚𝑛−𝑦𝑛−𝑦𝑛−1ℎ𝑛)对于第二类边界条件，由条件式可得𝑀0=𝑦0′′，𝑀𝑛=𝑦𝑛′′结合2𝑀0+𝑀1=6ℎ1(𝑦1−𝑦0ℎ1−𝑚0)𝑀𝑛−1+2𝑀𝑛=6ℎ𝑛(𝑚𝑛−𝑦𝑛−𝑦𝑛−1ℎ𝑛)可得2𝑀0+𝜆0𝑀1=𝑑1μ𝑛𝑀𝑛−1+2𝑀𝑛=𝑑𝑛其中𝑑0=6𝜆0ℎ1𝑓𝑥0,𝑥1−𝑚0+2(1−𝜆0)𝑀0𝑑𝑛=6𝜇𝑛ℎ𝑛𝑚0−𝑓𝑥𝑛−1,𝑥𝑛+2(1−𝜇𝑛)𝑀𝑛当𝜆0=𝜇𝑛=0时，式子2𝑀0+𝜆0𝑀1=𝑑1μ𝑛𝑀𝑛−1+2𝑀𝑛=𝑑𝑛即为𝑀0=𝑦0′′，𝑀𝑛=𝑦𝑛′′当𝜆0=𝜇𝑛=1时，式子2𝑀0+𝜆0𝑀1=𝑑1μ𝑛𝑀𝑛−1+2𝑀𝑛=𝑑𝑛即为2𝑀0+𝑀1=6ℎ1(𝑦1−𝑦0ℎ1−𝑚0)𝑀𝑛−1+2𝑀𝑛=6ℎ𝑛(𝑚𝑛−𝑦𝑛−𝑦𝑛−1ℎ𝑛)则可写成一个矩阵形式2𝜇1⋱𝜆02⋱𝜆1⋱𝜇𝑛−12𝜇𝑛𝜆𝑛−12𝑀0𝑀1⋮𝑀𝑛−1𝑀𝑛=𝑑0𝑑1⋮𝑑𝑛−1𝑑𝑛由于该方程组的系数是强对角占优，因此有唯一解对于第三类边界条件，由𝑦0=𝑦𝑛，𝑀0=𝑀𝑛，𝑆′𝑥0+0=𝑆′𝑥0−0可推出2𝑀0+𝜆0𝑀1+𝜇0𝑀n−1=𝑑0其中𝜆0=ℎ1ℎ1+ℎ𝑛，𝜇0=ℎ𝑛ℎ1+ℎ𝑛，𝑑0=6ℎ1+ℎ𝑛(𝑓𝑥0,𝑥1−𝑓𝑥𝑛−1,𝑥𝑛)2𝜆0𝜇0𝜇12𝜆1⋱⋱⋱𝜇𝑛−22𝜆𝑛−2𝜆𝑛−1𝜇𝑛−12𝑀0𝑀1⋮𝑀𝑛−2𝑀𝑛−1=𝑑0𝑑1⋮𝑑𝑛−2𝑑𝑛−1联立得方程组，且该方程组仅有唯一解。三次样条插值函数的性质性质一极小模性质设𝑓(𝑥)∈𝐶2[𝑎,𝑏]是任一被插值函数，𝑆(𝑥)是自然三次样条插值函数，则有[𝑆′′(𝑥)]2𝑏𝑎𝑑𝑥≤[𝑓′′(𝑥)]2𝑏𝑎𝑑𝑥式中等号当且仅当𝑓(𝑥)≡𝑆(𝑥)时成立例4.5.1设𝑓(𝑥)为定义在区间[0,3]上的函数，分割基点为𝑥1𝑖=0,1,2,3，并给出𝑓𝑥0=0，𝑓𝑥1=0.5，𝑓𝑥2=2.0，𝑓𝑥3=1.5和𝑓′𝑥0=0.2，𝑓′𝑥3=−1，试求三次样条插值函数𝑆(𝑥)，并使其满足边界条件性质二最佳逼近性质设𝑓(𝑥)∈𝐶2[𝑎,𝑏]是任一被插值函数，𝑆𝑓(𝑥)是带有斜率边界条件的三次样条插值函数，𝑆(𝑥)是与𝑆𝑓(𝑥)有相同分割的任一三次样条插值函数，则有[𝑓′′𝑥−𝑆𝑓′′(𝑥)]2𝑑𝑥≤[𝑓′′𝑥−𝑆′′(𝑥)]2𝑑𝑥𝑏𝑎𝑏𝑎性质三误差估计设函数𝑓(𝑥)∈𝐶4[𝑎,𝑏]，∆是区间[𝑎,𝑏]的一个分割，𝑆(𝑥)是关于𝑓𝑥的带有斜率边界条件或二阶导数边界条件的插值函数，则有插值估计max𝑎≤𝑥≤𝑏[𝑓𝑥−𝑆(𝑥)](𝑟)≤𝐶𝑟𝑀4ℎ4−𝑟其中𝐶0=5384，𝐶1=124，𝐶2=38，𝐶3=𝛽+𝛽−12，ℎ=max𝑖ℎ𝑖，𝛽=max𝑖ℎ𝑖min𝑖ℎ𝑖是分割比，并且系数𝐶0与𝐶1是最优估计