三次样条插值的MATLAB实现

1、MATLAB程序设计期中考查在许多问题中，通常根据实验、观测或经验得到的函数表或离散点上的信息， 去研究分析函数的有关特性。其中插值法是一种最基本的方法，以下给出最基本 的插值问题一一三次样条插值的基本提法：对插值区间a,b进行划分：a x0 x1xn b，函数y f x在节点Xi上的值 y f Xi i 0,1,2, n，并且如果函数 S x在每个小区间 Xi, Xi i上 是三次多项式，于a,b上有二阶连续导数，则称Sx是a,b上的三次样条函数， 如果S x在节点xi上还满足条件S Xiyii 0,1, n则称S x为三次样条插值函数。三次样条插值问题提法：对 a,b上给定的数表如下XX0

2、X1Xnyy。y1yn求一个分段三次多项式函数S x满足插值条件S xi yi i 0,1, n式，并在 插值区间a,b上有二阶连续导数。这就需要推导三次样条插值公式：记f x在节点xi处的值为f Xj mi ( i 0,1, n )(这不是给定插值问题数 表中的已知值)。在每个小区间Xi,Xi 1利用三次Hermite插值公式，得三次插值 公式：S x i x yi i 1 x yi 1i Xi mi 11， x Xi, Xi 1。为了得到这个公式需要4n个条件：(1).非端点处的界点有2n个；(2). 一阶导数连续有n 1个条件；(3).二阶导 数连续有n 1个条件，其中边界条件：O1 S

3、 X0m。S Xn S x00S xnYnXoS xnS x00 S xn其中:Xj0,1,i xj 0i xj0 且（i,j 0,1）。i Xj0,1,mi为对应变量的一阶导数。其推导过程如下:为了确定mi的值,展开为:Sx2X Xi 1 h 2 Xhi32xix xi hi 2 xi 1 xYi h?Yi i2X1 X Xihi22X XiXmi-厂hiXi 1mi 1,这里hi连续求两次导，得:6x 2xi4Xi 1hi2 mi6x 4xi 2xi 1hi2mi6 xixi 1 2xLYi1yi。于是考虑S在节点Xi处的右极限值，得:S xi0 mi mihihi6hi2Yi iYi 。

4、同理,在相邻小区间Xi 1,Xi上可得Sx的表达式为:6x 2xi 1 4xihTmi 16x 4xhi121 2xi6 xi 1 xi 2xmiYiYi i及S x在节点Xi处的左极限值为:2S xi0mi 1hi 1hi62 Yi1Yi。利用S x二阶导数于节点Xi处的连续性条件S Xi 0S xi这里有下式成立:1c 11mi 12mihi 1hi 1hi1mi 1 hiYi 1hi2YYiYi 1hi 12'用丄-除等式两hi 1 hi边，并注意yi £ , _yy f xi, xi !,上式可简记为:hiEi 12mim 1gii 1,2,n 1 ,且ihii1hi

5、i1gi 3i f Xi 1, Xii f Xi , Xihiihihi 1hi最后求得m1mn的线性方程组为：2i00 0 1叶g12 220 0 0m2g2(*)0 00n 1 2n 1mn 1gn 1n 000n 2mngn通过以上复杂的求解和迭代，就可以求解出插值函数的近似表达式。 得出来的表 达式就可以用MATLAB软件来求解。具体求解过程如下：已知n对数据点Xi,比，X2,y2，X3,y3，Xn，yn，假设函数关系为y f x，但解析式不确定，数据插值就是构造函数关系式y g x，使Xi i 1,2,3, ,n，满足关系 g Xif Xi。例题：求满足下面函数表所给出的插值条件的三

6、次自然样条函数。X1 245y f x1342分析：表中所列出的是函数对点，首先要把对应的插值函数求出来，再用MATLAB软件来求区间1,5上间隔为的各点的值 求解过程如下：因自然样条插值函数的边界条件为S Xo S Xn 0,这里n 3，故确定mo,mi,m2,m3的方程组形式形如上面的(*)式，其中系数hi :ho2 1 1,h142h12h21;i :1J2Jhoh13h1h：2 31i1112 1239gi:g131 f X1,X21 f X0,X12；7g232f X2,X32 f X1,X22；h°go3fX0,X12y。6;h2g33fX2,X32y36.将上述参数带入

7、(*)式，得到以下方程组：和gi可按如下步骤进行:h2235 41;210 0621m°92-033叶212m2702 -33m32001 26解得:5m2-49m38177mo,mi,84由公式2 .xxi 1hi2 x Xjhiyi2 .x xhi2 xi 1hixyi 1可知，1 3227,x 1,4XXx 1,S x8843 345 2103x 4,5XXx 33,884hi2mi2X Xi XX"-2 一 mi 1, hi2xx, 1 xxi由所求出的表达式可知区间1,5可分为1,44,5，对两个区间分别用MATLAB命令即可: 针对第一个区间：y 8x3 8x2 r 1;其图像如下3.5y2.51.51.52.53.5x0.5 1命令如下：x=1:4;y=(-1/8)*x43+(2/8)*x.A2+(7/4)*x-1;x