Interpolation
This presentation is the property of its rightful owner.
Sponsored Links
1 / 28

Interpolation PowerPoint PPT Presentation


  • 91 Views
  • Uploaded on
  • Presentation posted in: General

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;

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.While downloading, if for some reason you are not able to download a presentation, the publisher may have deleted the file from their server.


- - - - - - - - - - - - - - - - - - - - - - - - - - E N D - - - - - - - - - - - - - - - - - - - - - - - - - -

Presentation Transcript


Interpolation

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')


Interpolation

#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 );

}


Interpolation

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.


Interpolation

牛顿插值


Interpolation

% 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


Interpolation

分段多项式插值

三阶Hermite插值

区间 [xi, xi+1], s=x-xi, h=hi= xi+1-xi(Moler’s book,p.100)


Interpolation

区间 [xi, xi+1]


Interpolation

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


Interpolation

70年代,法国雷诺汽车公司的工程师贝齐尔(Bezier)创造出一种适用于几何体外形设计的新的曲线表示法。这种方法的优越性在于:对于在平面上随手勾画出的一个多边形(称为特征多边形),只要把其顶点坐标输入计算机,经过不到一秒钟的计算,绘图机就会自动画出同这个多边形很相像、又十分光滑的一条曲线。这种方法被人们称为贝齐尔(Bezier)方法(以下统称为Bezier方法)。

贝齐尔曲线的形状是通过一组多边折线(也称为贝齐尔控制多边形)的各顶点惟一地定义出来的。在该多边折线的各顶点中,只有第一点和最后一点在曲线上,其余的顶点则用来定义曲线的形状。图中列举了一些Bezier多边折线和相应的Bezier曲线的形状关系。


Interpolation

Bezier 曲线


Interpolation

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


Interpolation

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

Bezier曲线图例


Interpolation

B样条函数

  • 为了定义B样条曲线,首先给出n次截幂函数和n阶B样条函数的定义。我们称 为n次截幂函数,即

  • 称Mn(x)为n阶B样条函数,即


Interpolation

从 Bezier 曲线到B样条曲线

Bezier 曲线在应用中的不足:

缺乏灵活性。一旦确定了特征多边形的顶点数(m个),也就决定了曲线的阶次(m-1次),无法更改;

控制性差。当顶点数较多时,曲线的阶次将较高,此时,特征多边形对曲线形状的控制将明显减弱;

不易修改。由曲线的混合函数可看出,其值在开区间 ( 0 , 1 ) 内均不为零。因此,所定义之曲线在 ( 0 < t < 1)的区间内的任何一点均要受到全部顶点的影响,这使得对曲线进行局部修改成为不可能。而在外形设计中,局部修改是随时要进行的。

为了克服 Bezier 曲线存在的问题,Gordon 等人拓展了Bezier曲线,就外形设计的需求出发,希望新的曲线:易于进行局部修改;更逼近特征多边形;是低阶次曲线。

于是,用 n次B样条基函数替换了伯恩斯坦基函数,构造了称之为B样条曲线的新型曲线。


Interpolation

B样条函数图如图所示。

B样条函数


Interpolation

三阶Hermite插值

三次样条插值


Interpolation

区间 [xi, xi+1], (Moler’s book,p.102)

s=x-xi, yi=f(xi), di=f’(xi)

区间 [xi-1, xi],


Interpolation

  • 端点条件:

  • 固定斜率

  • 自然端点

  • 非节点

代入d3

另:三弯矩法


Interpolation

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.


Interpolation

MATLAB 函数

interp1(x,y,xt,’method’)

method = nearest

linear

spline

pchip

spline

pchip


Interpolation

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


Interpolation

yt = pchip(x,y,xt);

plot(x,y,'o',xt,yt,'-');

yt = spline(x,y,xt);

plot(x,y,'o',xt,yt,'-');


Interpolation

“It is important to prove,

but it is more

important to improve.”


  • Login