1、实验四 三次样条插值的应用一、问题描述The upper portion of this noble beast is to be approximated using clamped cubic spline interpolants. The curve is drawn on a grid from which the table is constructed. Use Algorithm 3.5 to construct the three clamped cubic splines. 二、模型建立三次样条插值给定一个列表显示的函数yi=y(xi),i=0,1,2,.,N-1。特别注意
2、在 xj 和 xj+1 之间的一个特殊的区间。该区间的线性插值公式为：(3.3.1)式和(3.3.2)式是拉格朗日插值公式(3.1.1)的特殊情况。因为它是(分段)线性的，(3.3.1)式在每一区间内的二阶导数为零，在横坐标为 xj 处的二阶导数不定义或无限。三次样条插值的目的就是要得到一个内插公式，不论在区间内亦或其边界上，其一阶导数平滑，二阶导数连续。做一个与事实相反的个假设，除 yi 的列表值之外，我们还有函数二阶导数y“的列表值，即一系列的 yi“值，则在每个区间内，可以在(3.3.1)式的右边加上一个三次多项式，其二阶导数从左边的 yj“值线性变化到右边的 yj+1“值，这么做便得到
3、了所需的连续二阶导数。如果还将三次多项式构造在 xj 和 xj+1 处为零，则不会破坏在终点 xj 和 xj+1 处与列表函数值 yj 和 yj+1 的一致性。进行一些辅助计算便可知，仅有一种办法才能进行这种构造，即用注意,(3.3.3)式和(3.3.4)式对自变量 x 的依赖，是完全通过 A 和 B 对 x 的线性依赖，以及 C 和 D(通过 A 和 B)对 x 的三次依赖而实现。可以很容易地验证，y“事实上是该插值多项式的二阶导数。使用 ABCD 的定义对 x 求(3.3.3) 式的导数，计算 dA/dx dB/dx dC/dx dD/dx,结果为一阶导数因为 x=xj 是 A=1，x=x
4、(i+1)时 A=0，而 B 正相反，则(3.3.6) 式表明 y“恰为列表函数的二阶导数。而且该二阶导数在两个区间(xj-1, xj)和(xj, xj+1)上是连续的。现在唯一的问题是，假设 yj“是已知的。而实际上并不知道。然而，仍不要求从(3.3.5)式算出的一阶导数在两个区间的边界处是连续的。三次样条的关键思想就在于要求这种连续性，并用它求得等式的二阶导数 yi“。设(3.3.5)式在区间(xj-1, xj)上对 x=xj 求得的值，等于同一等式在区间 (xj,xj+1)上对 x=xj 求得的值，便可得到所求方程，是新整理得到(对 j=1,N-2)这意味着有 N-2 个线性方程，但却有
5、 N 个未知数 yi“,i=0,N-1。因此，存在一个具有两个参数的可能解集。为求得唯一解，需要给出两个进一步的条件，一般取 x0 和 xn-1处的边界条件。最常见的做法有：1设 y0“和 yn-1“之一或两个都为零，得到所谓的自然三次样条函数，其一个或两个边界的二次导数为零。2设 yn“和 yn-1“为(3.3.5)式计算得到的值，使得该插值函数的一阶导数一个或两个边界处有特定的值。三次样条插值特别实用的原因之一在于，有两个附加边界条件的方程组(3.3.7)，它不仅是线性的，而且是三对角的。每个 yj“仅与其最邻近的 j+-1 的值有关。因此，方程可以用三对角算法在 O(N)次运算内求解。该
6、算法非常简明，很容易正确地构造出样条计算的程序。但是这使得程序不像(3.3.7)式的实现那样完全透明。三、模型求解借助 Matlab 软件进行编程求解，Matlab 代码如下X = 1, 2, 5, 6, 7, 8, 10, 13, 17, 20, 23, 24, 25, 27, 27.7, 28, 29, 30;Y = 3.0, 3.7, 3.9, 4.2, 5.7, 6.6, 7.1, 6.7, 4.5, 7.0, 6.1, 5.6, 5.8, 5.2, 4.1, 4.3, 4.1, 3.0cs = csapi(X,Y); %三次样条函数fnplt(cs);hold onplot(X,Y,o)legend(cubic spline,data)hold off最终生成的插值函数图像为四、实验感悟三次样条插值优点是计算简单，稳定性好，收敛性有保证且易在计算机上实现。缺点是它只能保证各小段曲线在连接点的连续性，却无法保证整条曲线的光滑性，这就不能满足某些工程技术的要求。但它依旧很强大，我们要好好掌握其精髓。
