实验四 三次样条插值.doc

实验四 三次样条插值的应用一、问题描述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。特别注意在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值，这么做便得到了所需的连续二阶导数。如果还将三次多项式构造在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(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个线性方程，但却有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)次运算内求解。该算法非常简明，很容易正确地构造出样条计算的程序。但是这使得程序不像(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最终生成的插值函数图像为四、实验感悟三次样条插值优点是计算简单