三次自然样条插值的统一表示及计算方法.docx

三次自然样条插值的统一表示及计算方法程义军 1 , 2 ,孙海燕 1( 1. 武汉大学 测绘学院 ,湖北 武汉 430079;2. 75719 部队 ,湖北 武汉 430074 )摘要 :三次自然插值样条常用分段表示法 ,通过三弯距方程进行求解 . 叙述了自然三次自然样条的统一表示方n 1法 , 即 g ( x) = a1 + a2 x + 12 6 i | x - xi | , 探讨了两种表示之间的联系与性质 , 论述了统一表示的计算方法 .3i = 0关键词 :插值 ;三次样条 ;自然样条中图分类号 : o241. 3文献标识码 : a文章编号 : 1008 - 8423 ( 2007) 01 - 0012 - 04样条函数具有光滑滑连续的特点 ,因而在工程应用及数据插值拟合中应用非常广泛 . 但通常情况下 ,三次样条都是采用分段表示的方法 ,即利用它具有二阶连续导数的特点 ,通过建立三弯距方程求得各节点上的 二阶导数 1 , 2 ,进而得到样条多项式在每一段上的系数. 这种方法具有很好的计算稳定性 ,但由于是分段表示 ,直观性不强 ,而且要求各节点按顺序排列 ,如果在中间增加一个节点 ,整个样条就得重新建立 .事实上 , 由于三次样条在整个区间上为三次多项式 , 从而有一种更为通用的表示方式 :n 1g ( x ) = a1 + a2 x + 12 6 i | x - xi | ,3i = 0这种表示适合于整个区间 , 而且它的一个最大优点就是节点不必按顺序排列 , 从而也更加灵活 , 也容易向多维扩展 35 . 而根据插值样条的唯一性理论 , 这种表示方式与分段表示法是等价的 , 因而两者之间也必存在 一定的联系 . 本文对两种表示方法间的联系与有关性质进行了分析 , 并探讨了 g ( x ) 的计算方法 .1 三次自然样条的分段表示定义 1 1, 2 设 a = x x x = b, 若函数满足 :0 1n( i) s ( x )在每个 xi - 1 , xi ( i = 1, 2, n )上是次数至多为三次的多项式 ;( ii) s ( x )在区间 a, b 上二阶导数连续 ;则称 s ( x )为 a, b 上关于给定分划的三次样条函数 . 特别在区间两端点处有 s( x0 ) = s( xn )为 a, b 上的三次自然样条函数 .= 0, 则称 s ( x )对于给定的函数值 yi = f ( xi ) ( i = 0, 1,三次样条插值函数., n ) , 若 s ( x )满足 s ( xi ) = yi ( i = 0, 1, n ) , 则称s ( x )为 f ( x )的由定义 1, s ( x )在每个小区间上为三次多项式 , 则 s( x )在每个小区间上为线性函数 . 设 s( xi ) =m i ( i =0, 1, n ) , 对于自然样条 m 0 =m n = 0, 利用插值条件 , 可得 s ( x )在每个小区间上 xi - 1 , xi ( i = 1, 2, n ) 的表达式为 1 :m i - 1 m i yi - 1 m i - 1 yi m i33s ( x ) =( xi- x ) +( x -xi - 1 )( xi - x ) +( x -xi - 1 )( 1 )+-6 hi6 hihi 6hi 6式中 hi = xi - xi - 1 6 由定义 1中的第二个条件 , s( x )在 xi 连续 , 则:收稿日期 : 2006 - 10 - 13.基金项目 :国家自然科学基金资助项目 ( 40274005 ) .作者简介 :程义军 ( 1971 - ) ,男 ,博士 ,工程师 ,主要从事测量数据处理理论与方法的研究.13第 1期程义军等 :三次自然样条插值的统一表示及计算方法 hi hi + hi +1 hi +1 1 1 11m i - 1+m i+m i +1=yi - 1+yi +1( 2 )-636hihi +1hihi +1( i = 1, 2, n - 1 )因为 m 0 =m n = 0, 上述方程组可表示成矩阵形式 :t( 3 )rm = q ytt其中 r 为 ( n - 1 )的方阵 , m = m 1 , m 2 , m n - 1 , y = y0 , y1 , yn , q 为 ( n + 1 ) ( n - 1 ) 的矩阵 , 且 q ij- 1- 1- 1- 1( i = 1, 2, n + 1; j = 1, 2, n - 1 )定义为 q ij = hj, q j + 1, j =- hj- hj + 1 , q j + 2, j = hj + 1且当 | i - j | 2 时 q ij = 1 10. 而 r ij ( i, j = 1, 2, n - 1 )的定义为 r ij = 3 ( hi + hi + 1 ) ( i = 1, 2, n - 1 ) , r i, i + 1 = r i + 1, i = 6 hi + 1 ( i = 1,n - 2 ) , 且当 | i - j| 2时 r ij = 0.显然 r 为严格对角占优的正定阵 , 由式 ( 3 )可解得 m 的惟一值 :- 1 t( 4 )m = r q y将 m 代入式 ( 1 ) , 即得样条函数 s ( x )的分段表达式.2 三次自然样条的统一表示由第 1部分可知 , 因为 r 是非奇异阵 , m 有惟一解 , 从而插值样条函数 s ( x ) 是惟一的 . 但 s ( x ) 的表示形式却可以是多样的. 与g ( x )表示.s ( x )的分段表示相对应的 , 有 s ( x )在整个区间上的统一的表达式 . 为示区别 , 下面用定理 1 3 x , a , a ,= , , t 为实数 , 若 g ( x )满足 :设 x x 0 1n1 20 1nn 1g ( x ) = a1 + a2 x + 12 6 i | x - xi | ,3i = 0nn且 6 i6 i xi= 0, 则 g ( x ) 为定义区间 x0 , xn 上的三次自然样条.i = 0证明i = 0显然 g ( x ) 在定义区间 ( x0 , xn ) 上为三次多项式 , 且有连续的二阶导数 , 从而 g ( x ) 为 ( x0 , xn ) 上n 1三次样条 . 当 x x0 时 , g ( x ) = a1 + a2 x + 12 6 i |3x -xi | , 从而 :i = 0nnng ( x ) = - 1 6 ( x -= - 1 ( x 6 + 6 x ) .x )iiii i22i = 0i = 0i = 0同理 , 当 x xn 时 , 所以 g ( x ) 为区间 x0 , xn 上的三次自然样条 . 13定义 e 为 ( n + 1 ) 维的方阵 , 且 eij = 12 | xi - 1 -, ( i, j = 1, 2, n + 1 ) , 显然有 eii = 0 ( i = 1, 2,xj- 1 |, n + 1 ) ; t 为 2 ( n + 1 ) 的矩阵 , 即 :1x01x1, n ) , 记1t =xnt设 g ( x )在各节点上的值为 gi = g ( xi ) ( i = 0, 1,g = g0 , g1 , gn , 则可表示为矩阵形式 :g = e + tt a( 5 )tt其中 = 0 ,1 ,n,n , a = a1 , a2 .n6 i6 i xi= 0也可表示为 t= 0, 因此有 :而=i = 0i = 0tetg0( 6 )=t0a因为各节点 xi 互不相同 , 由 e 的定义知 e 是满秩的 , 而阵 , 有惟一解 :a = ( t e tt ) - 1 t e- 1 gt 也是行满秩的 , 从而式 ( 6 ) 左端系数阵为满秩方( 7 ) = e- 1 ( g - tt a ) = e- 1 i - tt ( t e tt ) - 1 t e- 1 g14湖北民族学院学报 (自然科学版 )第 25卷求得 , a 后代入 g ( x )表达式 , 即为所求的三次自然插值样条 .根据插值多项式的惟一性 , 按上述方法所得的 g ( x )与 s ( x )是等价的 . 事实上有如下性质 .+-+-性质 1 0 = g ?( x0 ) ,n = - g ?( xn ) ,j = g ?( xj- g ?( xj ) ( j = 1, 2,根据 g ( x ) 的定义 , 有 :), n - 1 ) .nn+1+= a + a x + 6 ( x1 6 ( x+g ( x+ )x ) 3- x ) 3 ,-j1 2 ji jii ji1212i = 0i = j+1j- 1n-1-= a + a x + 6 ( x1 6 ( x+g ( x- )x ) 3- x ) 3 ,-j1 2 ji jii ji1212j = 0i = j分别求三次导数 , 得 :j- 1n1 6 + 11 6 ,g ?( x+ ) =-jiji222i = 0i = j+1j- 1ng ?( x- ) = 1 6 -11 6 ,-jiji222i = 0i = j+1+-两式相减得 j = g ?( xj )- g ?( xj ) . 由于在区间外g ( x ) = g ?( x ) = 0, 故有 :+-0 = g ?( x0 ) ,n= - g ?( xn ) .性质 2 = qm , q、m 定义同前 .因为在每个小区间上 g ?( x )为常数 , 从而有 :g ( xi + 1 ) - g ( xi ) g ( xi ) - g ( xi - 1 )=, g ?( x - ) =,g ?( x + )iihhi + 1i由性质 1即得 : g ( xi + 1 ) - g ( xi ) g ( xi ) - g ( xi - 1 )1=-, ( i = 1, 2, n - 1 ) .hi + 1 g ( xn - 1 )hi g ( x1 )同时有 0 =,n =. 将 i 的表达式进行适当的变换 :h1hn 11 1 1 i =g ( xi - 1 ) +-g ( xi ) +g ( x + 1 ) , ( i = 1, 2, n - 1 ) ,hihi + 1hihi + 1顾及 0 、n 及 q 的定义 , 则有 = qm.性质 3r = q t eq.由式 ( 3 )知 : rm = q t g = q t ( e+ tt a ) , 根据 q 与 t 的定义知 q t tt = 0, 顾及性质 2, 从而 : rm = q t e= q t eqm ,由于 m 的惟一性 , 所以必有 r = q t eq.上述 3个性质反映了三次自然插值样条两种表达方式之间的关系 . 性质 1表明了 i 为节点 xi 两边 g ( x )三次导数的增量 ; 性质 2表明若已知 m 则可计算 ( a 的计算见式 ( 9 ) ) , 反之若已知 , 则可利用 q 的特点递 推计算出 m. 由于 m 、r、e 均只与节点有关 , 性质 3则建立了这三个矩阵之间的关系 .3 , a的计算虽然与 s ( x )相比 , g ( x )的表达式更简洁 、直观 , 但 e 的性质却没有 r 好 . r 是一个严格对角占优的带状 矩阵 , 因此方程 ( 3 )很容易解算 (三对角阵的追赶法 1 ) . 虽然 e 是满秩阵 , 但其主对角线元素全为 0, 且当| i - j|越大时 , eij也越大 , 容易导致 e 的条件数变坏 , 而很坏的条件数会引起矩阵计算的不稳定 .上述性质虽然提供了一种计算 , a 的方法 , 但需要先求出 q、r; 虽是一种可行的办法 , 但不符合建立统 一表达式的本意. 因为在 s ( x )的定义与求解中要求各节点按大小顺序排列 , 而对 g ( x ) 却不必要 , 各节点可以是任意顺序 (但不能有相同的节点 ) , 从而灵活性更好 . 为解方程 ( 6 ) , 需要采用矩阵分解 .先对矩阵进行 q r 分解 6, 7 , 即 tt = fg, 其中 f 为 ( n + 1 ) ( n + 1 )的正交阵 , g 为 ( n + 1 ) 2 的上三角 阵 . 取 f1 为 f 的前 2列 , f2 为余下的 n - 1列 ; 并取 g1 为 g的前 2 行 , 因为 g为上三角阵 , 则剩余的均为 0,因此有 :15第 1期程义军等 :三次自然样条插值的统一表示及计算方法g1tt= f1 g1( 8 )f1 , f20因此由 t= 0, 有 gt ft= 0, 因为 g 为满秩的 , 从而得 ft= 0. 由于f 为正交阵 , 即 ft f = i、ft f = 0,1 1111 11 2令 = ft0, 因为 f ft 为单位阵 , 则有 = f . 设 w = ft g, w = ft g, 式 ( 6 ) 分别左乘ft 与ft , 并注意到22 221 1 2 212tf1= 0、= f2及 f 为正交阵 , 可得 :tg1 a = w 1 - f1 e( 9 )( 10 )tf2 e f2 = w 2将 ft ef 分解成 ud2 u t ( svd 分解 8 ) , 其中 u 为正交阵 , d 为对角阵且对角线上各元素满足d d 2 21 2dn - 3 0. 从而可得 : = u (d2 ) - 1 u tw( 11 )2进而得 : = f2因为 g1 的 2 2 上三角阵 , 由式 ( 9 )很容易求得 a.参考文献 : 1 2 3 4 5 6 朱方生 ,李大美 ,李素贞. 计算方法 m . 武汉 :武汉大学出版社 , 2003.程正兴 ,李水根. 数值逼近与常微分方程数值解 m . 西安 :西安交通常大学出版社 , 2000.green p j , silve rm an b w. nonp a ram e tric r egre ssion and gene ra lized l inea r mode lsm . london: cha pman & h a ll, 1994.彭振中. 半参数估计薄板样条法 d . 武汉 :武汉大学 , 2005.孙海燕 ,丁咚. 薄板样条函数及复杂曲面的数学表示 j . 测绘工程 , 2006 , 15 ( 2 ) : 4 - 14.p re ss w h , f lanne ry b p, teuko lsky s a , e t a l. n um e rica l r ec ip e s in for tran: the a rt of scien tific comp u ting ( 2 nd ed ) m . cam b ridge: cam b ridge u n iversity p ress, 1992: 107 - 110.刘丁酉. 矩阵分析 m . 武汉 :武汉大学出版社 , 2003.a nde rson e, b a i z, b ischof c, e t a l. la pack u se rs gu ide ( h ttp: / /www. netlib. o r