三次样条插值的MATLAB实现

MATLAB 程序设计期中考查

在许多问题中，通常根据实验、观测或经验得到的函数表或离散点上的信息，去研究分析函数的有关特性。其中插值法是一种最基本的方法，以下给出最基本的插值问题——三次样条插值的基本提法：

对插值区间[]b a ,进行划分：b x x x a n ≤<⋯⋯<<≤10，函数()x f y =在节点

i x 上的值()()n i x f y i i ⋯⋯==,2,1,0，并且如果函数()x S 在每个小区间[]1,+i i x x 上

是三次多项式，于[]b a ,上有二阶连续导数，则称()x S 是[]b a ,上的三次样条函数，如果()x S 在节点i x 上还满足条件

()()n i y x S i i ⋯⋯==,1,0

则称()x S 为三次样条插值函数。

三次样条插值问题提法：对[]b a ,上给定的数表如下.

求一个分段三次多项式函数()x S 满足插值条件()()n i y x S i i ⋯⋯==,1,0 式，并在

插值区间[]b a ,上有二阶连续导数。这就需要推导三次样条插值公式：

记()x f '在节点i x 处的值为()i i m x f ='(n i ⋯⋯=,1,0)(这不是给定插值问题数表中的已知值)。在每个小区间[]1,+i i x x 利用三次Hermite 插值公式，得三次插值公式：

()()()()1111+++++++=i i i i i i i i i m m x y x y x x S ββαα，[]1,+∈i i x x x 。为了得到这个公式需要n 4个条件：

(1).非端点处的界点有n 2个；(2).一阶导数连续有1-n 个条件；(3).二阶导数

连续有1-n 个条件，其中边界条件：○1()()n

n m x S m x S ='=' 00

○2()()αα=''=''n

x S x S 00 ○3()()()()1

6500403βααβαα=''+'=''+'n n x S x S x S x S

○4n y y =0 ()()()()000000-''=+''-'=+'n

n x S x S x S x S 其中：()⎩⎨⎧=≠=j i j i x j i ,1,0α ()0='j i x α ()0=j i x β 且(1,0,=j i )。

()⎩⎨

⎧=≠='j i j i x j i

,1,0β ，i m 为对应变量的一阶导数。其推导过程如下： 为了确定i m 的值，把()x S 展开为：

()()()[]

()()[]

13

123

2122+++-+-+

-+-=

i i

i i i i i

i i i y h x x h x x y h x x h x x x S

+

()()

()()

,12

122

21+++--+

--i i

i i i i

i i m h x x x x m h x x x x

这里i i i x x h -=+1，对()x S 连续求两次导，得：

()()

()i i i

i i i i

i i i i

i i y y h x x x m h x x x m h x x x x S --++

--+

--=

''+++++13

112

1

2

1

26246426。于是

考虑()x S ''在节点i x 处的右极限值，得： ()()i i i

i i i i i y y h m h m h x S -+--

=+''++1216

240。 同理，在相邻小区间[]i i x x ,1-上可得()x S ''的表达式为：

()()

()12

1

12

1

112

1

126246426----------++

--+

--=

''i i i i i i i i

i i i i

i y y h x x x m h x x x m h x x x x S

及()x S ''在节点i x 处的左极限值为：

()()12

1

1116

420-------+=

-''i i i i i i i i y y h m h m h x S 。利用()x S 二阶导数于节点i x 处的连续性条件()()00-''=+''i i x S x S ，这里1,2,1-⋯⋯=n i ，有下式成立：

⎪⎪⎭⎫ ⎝⎛-+-=+⎪⎪⎭⎫ ⎝⎛++--++---211211111311121i i i i

i i i i i i i i i h y y h y y m h m h h m h ，用i i h h 111+-除等式两边，并注意

[]11,,

++=-=i i i

i

i i i x x f h y y f y ，上式可简记为： (),1,2,1211-⋯⋯==+++-n i g m m m i i i i i i μλ 且[][]()1111

1,,31+----+=+=-=+=

i i i i i i i i

i i i i i i i i x x f x x f g h h h h h h μλλμλ

最后求得n m m ⋯⋯1的线性方程组为：

⎥⎥⎥⎥⎥⎥⎦

⎤⎢⎢⎢⎢⎢⎢⎣⎡⋯=⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎣⎡⋯⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎣⎡⋯⋯⋯⋯⋯⋯⋯⋯⋯⋯⋯⋯⋯⋯----n n n n n n n n g g g g m m m m 1211211122112000200000020002 λμμλμλλμ (**) 通过以上复杂的求解和迭代，就可以求解出插值函数的近似表达式。得出来的表达式就可以用MATLAB 软件来求解。具体求解过程如下：

已知n 对数据点()()()(),,,,,,,,332211n n y x y x y x y x ⋯⋯，假设函数关系为

()x f y =，但解析式不确定，数据插值就是构造函数关系式()x g y =，使()n i x i ,,3,2,1⋯⋯=，满足关系()()i i x f x g =。

例题：求满足下面函数表所给出的插值条件的三次自然样条函数。

分析：表中所列出的是函数对点，首先要把对应的插值函数求出来，再用

MATLAB 软件来求区间[]5,1上间隔为0.5的各点的值。 求解过程如下：

因自然样条插值函数的边界条件为