1 / 28

Interpolation

Interpolation. 向 华 武汉大学数学与统计学院. 整体多项式插值. Lagrangian polynomial interpolation. 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;

hisa
Download Presentation

Interpolation

An Image/Link below is provided (as is) to download presentation Download Policy: Content on the Website is provided to you AS IS for your information and personal use and may not be sold / licensed / shared on other websites without getting consent from its author. Content is provided to you AS IS for your information and personal use only. Download presentation by click this link. While downloading, if for some reason you are not able to download a presentation, the publisher may have deleted the file from their server. During download, if you can't get a presentation, the file might be deleted by the publisher.

E N D

Presentation Transcript

  1. Interpolation 向 华 武汉大学数学与统计学院

  2. 整体多项式插值 Lagrangian polynomial interpolation 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')

  3. #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 ); }

  4. Chebyshev interpolation 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.

  5. 牛顿插值

  6. % 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

  7. 分段多项式插值 三阶Hermite插值 区间 [xi, xi+1], s=x-xi, h=hi= xi+1-xi(Moler’s book,p.100)

  8. 区间 [xi, xi+1]

  9. 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)) );

  10. 70年代，法国雷诺汽车公司的工程师贝齐尔（Bezier）创造出一种适用于几何体外形设计的新的曲线表示法。这种方法的优越性在于：对于在平面上随手勾画出的一个多边形（称为特征多边形），只要把其顶点坐标输入计算机，经过不到一秒钟的计算，绘图机就会自动画出同这个多边形很相像、又十分光滑的一条曲线。这种方法被人们称为贝齐尔(Bezier)方法(以下统称为Bezier方法)。70年代，法国雷诺汽车公司的工程师贝齐尔（Bezier）创造出一种适用于几何体外形设计的新的曲线表示法。这种方法的优越性在于：对于在平面上随手勾画出的一个多边形（称为特征多边形），只要把其顶点坐标输入计算机，经过不到一秒钟的计算，绘图机就会自动画出同这个多边形很相像、又十分光滑的一条曲线。这种方法被人们称为贝齐尔(Bezier)方法(以下统称为Bezier方法)。 贝齐尔曲线的形状是通过一组多边折线（也称为贝齐尔控制多边形）的各顶点惟一地定义出来的。在该多边折线的各顶点中，只有第一点和最后一点在曲线上，其余的顶点则用来定义曲线的形状。图中列举了一些Bezier多边折线和相应的Bezier曲线的形状关系。

  11. Bezier 曲线

  12. 称为n次Bernstein多项式的基函数，Bezier曲线就是以此为基础构造出来的。称为n次Bernstein多项式的基函数，Bezier曲线就是以此为基础构造出来的。

  13. 2到8次Bezier曲线的图例。 Bezier曲线图例

  14. B样条函数 • 为了定义B样条曲线，首先给出n次截幂函数和n阶B样条函数的定义。我们称 为n次截幂函数，即 • 称Mn(x)为n阶B样条函数，即

  15. 从 Bezier 曲线到Ｂ样条曲线 Bezier 曲线在应用中的不足： 缺乏灵活性。一旦确定了特征多边形的顶点数(m个)，也就决定了曲线的阶次(m-1次)，无法更改； 控制性差。当顶点数较多时，曲线的阶次将较高，此时，特征多边形对曲线形状的控制将明显减弱； 不易修改。由曲线的混合函数可看出，其值在开区间 ( 0 , 1 ) 内均不为零。因此，所定义之曲线在 ( 0 < t < 1)的区间内的任何一点均要受到全部顶点的影响，这使得对曲线进行局部修改成为不可能。而在外形设计中，局部修改是随时要进行的。 为了克服 Bezier 曲线存在的问题，Gordon 等人拓展了Bezier曲线，就外形设计的需求出发，希望新的曲线:易于进行局部修改;更逼近特征多边形;是低阶次曲线。 于是，用 n次Ｂ样条基函数替换了伯恩斯坦基函数，构造了称之为Ｂ样条曲线的新型曲线。

  16. B样条函数图如图所示。 B样条函数

  17. 三阶Hermite插值 三次样条插值

  18. 区间 [xi, xi+1], (Moler’s book,p.102) s=x-xi, yi=f(xi), di=f’(xi) 区间 [xi-1, xi],

  19. 端点条件: • 固定斜率 • 自然端点 • 非节点 代入d3 另：三弯矩法

  20. Cubic splines in general don't preserve monotonicity between neighbouring nodes; but Hermite piecewise cubic interpolation does. 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) ] ); 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.

  21. MATLAB 函数 interp1(x,y,xt,’method’) method = nearest linear spline pchip spline pchip

  22. 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,'-');

  23. yt = pchip(x,y,xt); plot(x,y,'o',xt,yt,'-'); yt = spline(x,y,xt); plot(x,y,'o',xt,yt,'-');

  24. “It is important to prove, but it is more important to improve.”

More Related