《数值分析(样条插值)总结》由会员分享,可在线阅读,更多相关《数值分析(样条插值)总结(35页珍藏版)》请在金锄头文库上搜索。
1、样条插值 多元插值,数值分析 14, ,是否可以在光滑性和局部单调性之间折衷呢?,多项式插值是一个极端, 它可以进行无限次的微分, 但它通常不能保持给定数据所描述的形状, 特别是在端点附近。分段线性插值是另一个极端, 它几乎没有任何光滑性。它连续但一阶导数存在跳变。另一方面它保持了给定数据的局部单调性。,定义 5.4 给定区间a , b上的一个分划: a = x0 x1 xn = b 已知 f(xj) = yj (j = 0,1,n), 如果,S(x) 为在 xj,xj+1的三次多项式满足: (1) Sj (xj-1) = yj-1 , Sj (xj) = yj , ( j = 1,n) (2
2、) Sj (xj) = Sj+1 (xj) ( j = 1,n-1) (3) Sj (xj) = Sj+1 (xj) ( j = 1,n-1) 则称 S(x)为三次样条插值函数。,三次样条插值函数满足的连续条件: (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)二阶导数连续,回顾1:,当xxj , xj+ 1 ( j= 0,1,n-1 )时 Sj(x)= aj + bj x + cj x2 + dj x3,插值条件: S(xj) = yj
3、( j = 0,1,n) 连续性条件: S(xj+0) =S(xj0) ( j = 1,n-1) S(xj+0) =S (xj0) ( j = 1,n-1) S (xj+0) =S (xj0) ( j = 1,n-1),由样条定义,可建立方程(4n-2)个!,n个三次多项式, 待定系数共4n个!,方程数少于未知数个数 ?,(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)=
4、f(xn),例 5.7 已知f(1) = 1, f(0) = 0, f(1) = 1。构造分段三次多项式是满足自然边界的样条函数。,证:显然,求导数得,显然,分段Hermite插值公式导出的样条方法,已知函数表,设 f(x) 在各插值节点 xj 处的一阶导数为 mj,取 xj+1 xj = h,( j = 0,1,2,n-1)。当 xxj, xj+ 1时, 分段三次Hermite插值,品味,取 xj+1 xj = h,( j = 0,1,2,n-1)。 当 xxj, xj+ 1时, 分段Hermite插值,当xxj , xj+ 1 ( j= 0,1,n-1 )时 Sj(x)= aj + bj
5、x + cj x2 + dj x3,未知数个数 (n+1)!,由S (x)连续,有等式 S(xj + 0)=S(xj 0),当 xxj , xj+1时, S(x) 由基函数组合而成,同理有,联立得,( j=1, 2, , n-1 ),自然样条的导数值满足:,设自然边界条件成立即,( j=1, 2, , n-1 ),自然样条的导数值满足:,( j=1, 2, , n-1 ),严格对角占优矩阵,回顾: 严格主对角占优矩阵一定是非奇异的。,程序片段1:,Matlab Code : 三次样条插值 function v = splinetx(x,y,u) % First derivatives h =
6、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,程序片段2:,Ma
7、tlab 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 syste
8、m 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);,程序片段3:,Matlab Code : 三次样条插值 % Right-hand side r(1) = (h(1)+2*c(1)*h(2)*delta(1)+ . h(1)2*delta(2)/
9、c(1); r(2:n-1) = 3*(h(2:n-1).*delta(1:n-2)+ . h(1:n-2).*delta(2:n-1); r(n) = (h(n-1)2*delta(n-2)+ . (2*a(n-1)+h(n-1)*h(n-2)*delta(n-1)/a(n-1); % Solve tridiagonal linear system d = tridisolve(a,b,c,r);,误差估计,多项式插值是一个极端, 它可以进行无限次的微分, 但它通常不能保持给定数据所描述的形状, 特别是在端点附近。分段线性插值是另一个极端, 它几乎没有任何光滑性。它连续但一阶导数存在跳变。另
10、一方面它保持了给定数据的局部单调性。,三次样条插值函数满足的连续条件: (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)二阶导数连续,S(xj + 0)=S(xj 0),Runge反例 f(x)=1/(x2+1), (-5=x=5),x=-5:5; y=1./(x.2+1); x1=-5:.01:5; y1=1./(x1.2+1); u=-5:.01:5; v1=polyinterp(x,y,u); plot(x1,y1,g-,u,v1,b
11、-) v2=splinetx(x,y,u); plot(x1,y1,g-,u,v2,r-),一般认为Ln(x)的次数n越高逼近f(x)的精度越好?,interpgui,样条插值是插值函数的光滑性与局部单调性之间的折衷方案。,函数的范数,样条插值函数的极性,设f(x)C2a, b, 对于a = x0 x1 xn = b,有 f(xj)=yj(j=0,1,n)。S(x)是满足S(xj)=yj(j=0,1,n) 的三次自然样条。则有 |S (x)|2|f (x)|2,样条函数S(x)在a, b上总曲率最小。理论上确保了样条插值能产生尽可能平稳的函数, 而不会像高阶多项式插值那样出现巨大的震荡。,即,
12、证明:,矩形区域上函数f(x, y)的双线性插值。,插值条件: P(x1, y1) = z1, P(x2, y1) = z2, P(x2, y2) = z3, P(x1, y2) = z4,P(x, y) = ax + by + cxy + d,l1(u,v)= (1 u)(1 v) l2(u,v)= u(1 v) l3(u,v)= u v l4(u,v)= (1 u) v,其中,P(x, y) = z1(1 u)(1 v)+ z2 u(1 v) + z3 u v + z4 (1 u)v,二元函数插值,一维数据的插值 interp1(linear, cubic, spline),数据与插值函数
13、,x=0:2:24; y=22,21,19,18,20,24,27,32,31,28,26,23,22; xi=0:0.1:24; yi=interp1(x,y,xi,spline); plot(xi,yi,-,x,y,o);,二维网格数据的插值 interp2 (linear, cubic, spline),x,y=meshgrid(1:5,1:3); z=82 81 80 82 84;79 63 61 65 81;84 84 82 85 86; mesh(x,y,z) xi,yi=meshgrid(1:0.2:5,1:0.2:3); zi=interp2(x,y,z,xi,yi,splin
14、e); mesh(xi,yi,zi);,数据与插值函数,散乱数据的插值 griddata, ndgrid (scattered data ),rng(0,twister) x = rand(100,1)*4-2; y = rand(100,1)*4-2; z = x.*exp(-x.2-y.2); ti = -2:.25:2; xq,yq = meshgrid(ti,ti); zq = griddata(x,y,z,xq,yq); mesh(xq,yq,zq), hold on, plot3(x,y,z,o),图像放大,I = imread(circuit.tif); J = imresize
15、(I,1.25);% nearest,bilinear ,bicubic imshow(I) figure, imshow(J),总结:,完整多项式插值(拉格朗日形式) Hermite插值(函数值和导数值插值) 分段线性插值 分段三次样条插值 多元函数插值,插值函数的构造化为基函数(building blocks)的构造,所有的现象背后都有一个规律(函数)吗? 插值方法可以用来预测(外推)吗?,美国Wired杂志 Mathematician Predicts Who Will Live and Die in Game of Thrones When Extrapolation Fails Us: Incorrect Mathematical Conjectures,