Transcript E:WHU教学2010Lecture2008lec8_interpLec8Interpolation

Interpolation
向
华
武汉大学数学与统计学院
整体多项式插值
Lagrangian polynomial interpolation
f ( n1) ( ) n
f ( x)   n f ( x) 
( x  xi )

(n  1)! i 0
x = linspace(-5,5,13);
y = 1./(1+x.^2);
c = polyfit(x,y,12);
t = linspace(x(1),x(end),100);
p = polyval(c,t);
plot(t,p,x,y,'o')
#include <stdio.h>
#include <math.h>
void main()
{
int i, j, kk;
float xa, yans, z;
static n = 3;
/* n+1 is the number of data points. */
static float x[11]={1. , 2. , 3. , 4. };
static float f[11]={.671, .620, .567, .512};
printf( "\nInput x ? " ); scanf( "%f", &xa );
yans = 0;
for( i = 0; i <= n; i++ ) {
z = 1.0;
for( j = 0; j <= n; j++ ) {
if( i != j ) z = z*(xa - x[j])/(x[i] - x[j]);
}
yans = yans + z*f[i];
}
printf( "Answer: g( %g ) = %g \n", xa, yans );
}
Chebyshev interpolation
ab ba
i
xi 

cos , i  0,..., n.
2
2
n
n = 8;
a = -5;
b = 5;
x = (a+b)/2 + (b-a)/2*
(-cos(pi*[0:n]/n));
f = '1./(1+x.^2)';
y = eval(f);
c = polyfit(x,y,n);
t = linspace(x(1),x(end),100);
p = polyval(c,t);
plot(t,p)
hold on;
fplot('1/(1+x^2)',[-5,5]);
Chebyshev interpolation :Runge's phenomenon can be avoided if a suitable
distribution of nodes is used.
牛顿插值
P3 ( x)  c1  c2 ( x  x1 )  c3 ( x  x1 )(x  x2 )  c4 ( x  x1 )(x  x2 )(x  x3 )
c1  yi
c2  f [ x1 , x2 ]
c3  f [ x1 , x2 , x3 ]
c4  f [ x1 , x2 , x3 , x4 ]
P3 ( x)  P2 ( x)  f [ x1, x2 , x3 , x4 ](x  x1 )(x  x2 )(x  x3 )
xi
f[]
f [,]
f [,,]
x1
f[x1]
x2
f[x2]
f[x1,x2]
x3
f[x3]
f[x2,x3]
f[x1,x2,x3]
x4
f[x4]
f[x3,x4]
f[x2,x3,x4]
f [,,,]
f[x1,x2,x3,x4]
P3 ( x)  c1  ( x  x1 )c2  ( x  x2 )c3  c4 ( x  x3 )
% Interpolation with Newton polynomials
% Input: x,y: y=f(x)
%
xt: where interpolant is evaluated.
% Output: yt=f(xt)
function yt = Newton_interpol(x,y,xt)
n = length(y); if length(x)~=n, error('x and y are not compatible'); end
c=y(:);
for j=2:n
for i=n:-1:j
c(i)=( c(i)-c(i-1) ) / ( x(i) - x(i-j+1) ); % note that x(i)-x(i-1) is not right.
end
end
% Nested evaluation of the polynomial
yt = c(n);
for i= n-1 :-1 :1
yt = yt.*(xt - x(i)) + c(i);
end
分段多项式插值
三阶Hermite插值
区间 [xi, xi+1], s=x-xi, h=hi= xi+1-xi (Moler’s book,p.100)
3hs 2  2s 3
h 3  3hs 2  2s 3
Pi ( x) 
f ( xi 1 ) 
f ( xi )
3
3
h
h
s 2 ( s  h) '
s ( s  h) 2 '

f ( xi 1 ) 
f ( xi )
2
2
h
h
区间 [xi, xi+1]
Pi ( x)   i   i ( x  xi )  i ( x  xi ) 2  i ( x  xi )3
 i  f ( xi )
 i  f ' ( xi )
3 f [ xi , xi 1 ]  2 f ' ( xi )  f ' ( xi 1 )
i 
xi 1  xi
f ' ( xi )  2 f [ xi , xi 1 ]  f ' ( xi 1 )
i 
2
xi 1  xi 
function yt = hermite_interpol(x,f,fp,xt)
n = length(x);
if length(f)~=n | length(fp)~=n, error('not compatible'); end
dx = diff(x);
divDiff = diff(f)./dx;
a = f(1:n-1);
b = fp(1:n-1);
c = (3*divDiff - 2*fp(1:n-1) - fp(2:n) ) ./dx;
d = ( fp(1:n-1) - 2*divDiff + fp(2:n) ) ./dx.^2;
% Locate each xt value in the x vector
K = zeros(size(xt));
for m = 1: length(xt)
K(m)=biSearch(x,xt(m)); % Binary search to find index i s.t. x(i)<= xt <=x(i+1)
end
% vectorized evaluation of the piecewise polynomials
xx = xt-x(K);
yt = a(K) + xx.* ( b(K) + xx.*(c(K)+xx.*d(K)) );
f ( x)  xe x
70年代，法国雷诺汽车公司的工程师贝齐尔（Bezier）
创造出一种适用于几何体外形设计的新的曲线表示法。
这种方法的优越性在于：对于在平面上随手勾画出的
一个多边形（称为特征多边形），只要把其顶点坐标
输入计算机，经过不到一秒钟的计算，绘图机就会自
动画出同这个多边形很相像、又十分光滑的一条曲线。
这种方法被人们称为贝齐尔(Bezier)方法(以下统称为
Bezier方法)。
贝齐尔曲线的形状是通过一组多边折线（也称为贝齐尔
控制多边形）的各顶点惟一地定义出来的。在该多边
折线的各顶点中，只有第一点和最后一点在曲线上，
其余的顶点则用来定义曲线的形状。图中列举了一些
Bezier多边折线和相应的Bezier曲线的形状关系。
Bezier 曲线
注意恒等式
n k
nk


1＝（t  1  t )    t (1  t )
k 0  k 
假定f  c0,1，将它写成
n
n
n
f (t )  
k 0
n k
nk
f (t ) t (1  t )
k 
n
n!
其中 是二项式系数
，我们称
k!(n  k )!
k 
n
Bn f  
k 0
 k  n  k
nk
f   t (1  t )
 n  k 
k
  f  Bk ,n (t )
(0  t  1)
k 0  n 
为f (t )的n次Bernstein 多项式, 其中
n
n k
Bk ,n (t )   t (1  t ) n  k
k 
称为n次Bernstein多项式的基函数，
Bezier曲线就是以此为基础构造出来的。
2到8次Bezier
Bezier
曲线图例
B样条函数

为了定义B样条曲线，首先给出n次
截幂函数和n阶B样条函数的定义。我们
n
称 x 为n次截幂函数，即
n

x ,
n
x  
 0,
x0
x0
 称Mn(x)为n阶B样条函数，即
n
1
n
k n
n 1


M n ( x) 
(1)  ( x   k ) 

(n  1)! k 0
2
k 
从 Bezier 曲线到Ｂ样条曲线
Bezier 曲线在应用中的不足：
缺乏灵活性。一旦确定了特征多边形的顶点数(m个)，也就决定了曲线的
阶次(m-1次)，无法更改；
控制性差。当顶点数较多时，曲线的阶次将较高，此时，特征多边形对曲
线形状的控制将明显减弱；
不易修改。由曲线的混合函数可看出，其值在开区间 ( 0 , 1 ) 内均不
为零。因此，所定义之曲线在 ( 0 < t < 1)的区间内的任何一点均要受
到全部顶点的影响，这使得对曲线进行局部修改成为不可能。而在外形设
计中，局部修改是随时要进行的。
为了克服 Bezier 曲线存在的问题，Gordon 等人拓展了 Bezier曲线，就
外形设计的需求出发，希望新的曲线:易于进行局部修改;更逼近特征多边
形;是低阶次曲线。
于是，用 n次Ｂ样条基函数替换了伯恩斯坦基函数，构造了称之为Ｂ样条
曲线的新型曲线。
B样条函数图如图所示。
B样条函数
三阶Hermite插值
Pi ( x)   i   i ( x  xi )  i ( x  xi ) 2  i ( x  xi )3
三次样条插值
Pi ' ( x)   i  2i ( x  xi )  3i ( x  xi ) 2
 i  f ( xi )
i  ?
3 f [ xi , xi 1 ]  2 i   i 1
i 
xi 1  xi
i 
 i  2 f [ xi , xi 1 ]   i 1
xi 1  xi 2
区间 [xi, xi+1], (Moler’s book,p.102)
3hs 2  2s 3
h3  3hs 2  2s 3
s 2 ( s  h)
s ( s  h) 2
Pi ( x) 
yi 1 
yi 
d i 1 
di
3
3
2
2
h
h
h
h

6h  12s  i  6s  2h d i 1  6s  4h d i
P ( x) 
"
i
Pi" ( xi ) 
h2
6 i  2d i 1  4d i
hi
i 
yi 1  yi
xi 1  xi
 6 i  4d i 1  2d i
Pi ( xi 1 ) 
hi
区间 [xi-1, xi],
Pi"1 ( xi ) 
s=x-xi, yi=f(xi),
di=f’(xi)
 6 i 1  4di  2di 1
hi 1
"
Pi" ( xi )  Pi"1 ( xi )
hi di 1  2hi 1  hi di  hi 1di 1  3hii 1  hi 1 i 
di 1  4di  di 1  3i 1  3 i
i  2,, n  1
端点条件:
1. 固定斜率
2d1  d1  31 dn1  2dn  3 n1
2. 自然端点
P1" ( x1 )  0
3. 非节点
P1''' ( x2 )  P2''' ( x2 )
1  2
d1  21  d2  d2  2 2  d3
代入d3
另：三弯矩法
2d1  4d2  51   2
t = linspace(0, pi/2, 4);
x = cos(t);
y = sin(t);
xt = linspace(0,1, 40);
plot(x,y,'s', xt, [pchip(x,y,xt); spline(x,y,xt) ] );
Cubic splines in general don't preserve
monotonicity between neighbouring nodes;
but Hermite piecewise cubic interpolation does.
Spline is more accurate if the data are values
of a smooth function. Pchip has no overshoots
and less oscillation if the data are not smooth.
MATLAB 函数
interp1(x,y,xt,’method’)
method = nearest
linear
spline
pchip
spline
pchip
x = linspace(0, 1.5 ,10);
y = humps(x);
xt = linspace(min(x), max(x));
yt = interp1(x,y,xt,'nearest');
plot(x,y,'o',xt,yt,'-');
yt = interp1(x,y,xt,'linear');
plot(x,y,'o',xt,yt,'-');
yt = pchip(x,y,xt);
plot(x,y,'o',xt,yt,'-');
yt = spline(x,y,xt);
plot(x,y,'o',xt,yt,'-');
“It is important to prove,
but it is more
important to improve.”