数值分析(样条插值)总结

《数值分析》 14
样条插值
多元插值
18:23
1/35
多项式插值是一个极端, 它可以进行无限次的微分, 但它通常不能保持 给定数据所描述的形状 , 特别是在端点附近。分段线性插值是另一个极端 , 它几乎没有任何光滑性。它连续但一阶导数存在跳变。另一方面它保持了 给定数据的局部单调性。
18:23
4 2 4 j ( x j ) [ 2 ( x x j 1 ) ( x x j ) 2 ] x x j h h h 2 ( x ) [ 4 ( x x ) ( x x ) 2 ] j 1 j j j 1 2 2 x x j h h h
S(x) 为在 [xj,xj+1]的三次多项式满足:
(1) Sj (xj-1) = yj-1 , Sj (xj) = yj , ( j = 1,· · · ,n) (2) S'j (xj) = S'j+1 (xj) ( j = 1,· · · ,n-1)
(3) · · ,n-1) 18:23 S'' (x ) = S'' j j j+1 (xj) ( j = 1,· 则称 S(x)为三次样条插值函数。
2
S j ( x ) (1 2
x j 1 x
)(
x xj
)2 y j 1
18:23
8/35
品味
当x∈[xj , xj+ 1] ( j= 0,1,…n-1 )时 Sj(x)= aj + bj x + cj x2 + dj x3
取 xj+1 – xj = h，( j = 0,1,2,· · · ,n-1)。 当 x∈[xj, xj+ 1]时, 分段Hermite插值
18:23
方程数少于未知数个数 ？？
5/35
(1)自然边界条件: S'' (x0)=0, S'' (xn)=0
(2)周期边界条件: S'(x0)=S'(xn), S''(x0)=S'' (xn) (3)反射边界条件: S''(x0)=S''(x1), S''(xn-1)=S'' (xn) (4) 固定边界条件: S'(x0)=f '(x0), S'(xn)=f'(xn)
是否可以在光滑性和局部单调性之间折衷呢 ? 2/35
定义 5.4 给定区间[a , b]上的一个分划:
a = x0 < x 1 < … < x n = b · · ,n), 如果 已知 f(xj) = yj (j = 0,1,·
S1 ( x ), x [ x0 , x1 ] S ( x ), x [ x , x ] 2 1 2 S( x) S n ( x ), x [ x n1 , xn ]
6 6 4 2 S ( x0 0) 2 y0 2 y1 m0 m1 0 h h h h
6 6 2 4 S ( xn 0) 2 yn1 2 yn mn1 mn 0 h h h h
自然样条的导数值满足:
3 2m0 m1 [ y1 y0 ] h
11/35
6 6 4 2 S ( x j 0) 2 y j 2 y j 1 m j m j 1 h h h h
同理有 联立得
6 6 2 4 S ( x j 0) 2 y j 1 2 y j m j 1 m j h h h h
6 6 4 2 2 y j 2 y j 1 m j m j 1 h h h S( x j 0) S( x j 0) h
18:23
14/35
回顾: 严格主对角占优矩阵一定是非奇异的。
严格对角占优矩阵定义 :| aii | | aij |
i j
反证法, 如果A奇异, 则Ax =0存在非零解x。 注意v x / x 仍然为Ax =0的非零解且 v 则存在m ,|vm | 1且|v j| 1。 注意对于任意的i满足aii vi aij v j 0。
18:23
6/35
例 5.7 已知f(–1) = 1, f(0) = 0, f(1) = 1。构造分 段三次多项式是满足自然边界的样条函数。
1 3 3 2 x x , x [1, 0] 2 2 S ( x) 1 x 3 3 x 2 , x [0, 1] 2 2

18:23
6 6 2 4 2 y j 1 2 y j m j 1 m j h h h h 3 m j 1 4m j m j 1 ( y j 1 y j 1 ) h
( j=1, 2, · · · · · · , n-1 )
12/35
设自然边界条件成立即
m j 1 4m j m j 1
3 ( y j 1 y j 1 ) h
3 m n 1 2 m n [ yn yn 1 ] h
18:23
( j=1, 2, · · · · · · , n-1 )
13/35
自然样条的导数值满足:
3 2m0 m1 [ y1 y0 ] h
m j 1 4m j m j 1
3 m n 1 2 m n [ yn yn 1 ] h
3 ( y j 1 y j 1 ) h
( j=1, 2, · · · · · · , n-1 )
2 1 m0 3[ y1 y0 ] / h 1 4 1 m 3( y y ) / h 3 1 1 严格对角占优矩阵 1 4 1 mn 1 3( yn yn 2 ) / h 1 2 mn 3( yn yn 1 ) / h
证:显然 求导数得
S ( 1) 1
S (1) 1 S (0 ) S (0 ) 0
显然
3 2 x 3 x , x [1, 0] 3 x 3, x [1, 0] 2 S ( x ) S ( x ) 3 x 3, x [0, 1] 3 x 2 3 x , x [0, 1] 2 3 3
回顾1: 若 lim f ( x ) f ( x0 ), 则称函数在x0左连续;
18:23
x x0
若 lim f ( x ) f ( x0 ), 则称函数在x0 右连续。
x x0 +
4/35
n个三次多项式, 待定系数共4n个！
当x∈[xj , xj+ 1] ( j= 0,1,…n-1 )时 Sj(x)= aj + bj x + cj x2 + dj x3 插值条件: S(xj) = yj ( j = 0,1,· · · ,n ) 连续性条件: S(xj+0) =S(xj–0) ( j = 1,· · · ,n-1) S'(xj+0) =S' (xj–0) ( j = 1,· · · ,n-1) S'' (xj+0) =S' ' (xj–0) ( j = 1,· · · ,n-1) 由样条定义,可建立方程(4n-2)个！
S ( x j 0) j ( x j ) y j j 1 ( x j ) y j 1 j ( x j )m j j 1 ( x j )m j 1
18:23
6 6 4 2 S ( x j 0) 2 y j 2 y j 1 m j m j 1 h h h h
h x xj 2 j 1 ( x ) ( x x j 1 )( ) h
18:23
j ( x ) ( x x j )(
x j 1 x
)2
10/35
x xj 2 8 j( x j ) [ 3 ( x j 1 x ) (1 2 ) 2 ]x x j 2 h h h h 6 ( x ) [ 8 ( x x ) (1 2 x j 1 x ) 2 ] 2 j 1 j j 3 2 x x j h h h h
3/35
三次样条插值函数满足的连续条件:
(1) S(xj–)= S(xj+) ( j = 1,· · · ,n-1) 连续
(2) S' (xj–)= S'(xj+) ( j = 1,· · · ,n-1)导数连续
(3) S'' (xj–)= S'' (xj+) ( j = 1,· · · ,n-1)二阶导数连续
ji
1。
故amm vm amj v j
jm
18:23 |a 进而 mm |=|amm v m | | amj v j | | amj |矛盾。
j来自百度文库m
jm
15/35
程序片段1:
Matlab Code : 三次样条插值 function v = splinetx(x,y,u) % First derivatives h = diff(x); delta = diff(y)./h; d = splineslopes(h,delta); % Piecewise polynomial coefficients n = length(x); c = (3*delta - 2*d(1:n-1) - d(2:n))./h; b = (d(1:n-1) - 2*delta + d(2:n))./h.^2; % Find subinterval indices k so that x(k) <= u < x(k+1) k = ones(size(u)); for j = 2:n-1 k(x(j) <= u) = j; end
S j ( x ) (1 2 x xj h h h h x j 1 x 2 x xj 2 ( x x j )( ) m j ( x x j 1 )( ) m j 1 h h )( x j 1 x ) y j (1 2
2
x j 1 x
)(
x xj
)2 y j 1
18:23
未知数个数 (n+1)！！
9/35
由S ''(x)连续,有等式 S''(xj + 0)=S''(xj – 0)
当 x∈[xj , xj+1]时, S(x) 由基函数组合而成
j ( x ) (1 2
x xj h h x j 1 x x x j 2 j 1 ( x ) (1 2 )( ) h h )( x j 1 x )2
S ( 1)
18:23
S (1) 0
2 2 S (1) 0 S (0 ) S (0 ) 3
S (1)
S (0 ) S (0 ) 0
7/35
分段Hermite插值公式导出的样条方法
已知函数表
x f(x)
x0 y0
x1 y1
· · · · · ·
xn yn
设 f(x) 在各插值节点 xj 处的一阶导数为 mj
取 xj+1 – xj = h，( j = 0,1,2,· · · ,n-1)。当 x∈[xj, xj+ 1]时,
分段三次Hermite插值
x xj h h h h x j 1 x 2 x xj 2 ( x x j )( ) m j ( x x j 1 )( ) m j 1 h h )( x j 1 x ) y j (1 2
18:23
16/35
程序片段2:
Matlab Code : 三次样条插值 % Evaluate spline s = u - x(k); v = y(k) + s.*(d(k) + s.*(c(k) + s.*b(k))); % ------------------------------------------------------function d = splineslopes(h,delta) % SPLINESLOPES Slopes for cubic spline interpolation. % splineslopes(h,delta) computes d(k) = S'(x(k)). % Uses not-a-knot end conditions. % Diagonals of tridiagonal system n = length(h)+1; a = zeros(size(h)); b = a; c = a; r = a; a(1:n-2) = h(2:n-1); a(n-1) = h(n-2)+h(n-1); b(1) = h(2); b(2:n-1) = 2*(h(2:n-1)+h(1:n-2));b(n) = h(n-2); c(1) = h(1)+h(2);c(2:n-1) = h(1:n-2);