slide1
Download
Skip this Video
Download Presentation
Interpolation

Loading in 2 Seconds...

play fullscreen
1 / 28

Interpolation - PowerPoint PPT Presentation


  • 135 Views
  • Uploaded on

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;

loader
I am the owner, or an agent authorized to act on behalf of the owner, of the copyrighted work described.
capcha
Download Presentation

PowerPoint Slideshow about 'Interpolation' - hisa


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
slide1
Interpolation

向 华

武汉大学数学与统计学院

slide2
整体多项式插值

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

slide3
#include

#include

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

}

slide4
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.

slide7
% 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

slide8
分段多项式插值

三阶Hermite插值

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

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

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

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

slide17
B样条函数
  • 为了定义B样条曲线,首先给出n次截幂函数和n阶B样条函数的定义。我们称 为n次截幂函数,即
  • 称Mn(x)为n阶B样条函数,即
slide18
从 Bezier 曲线到B样条曲线

Bezier 曲线在应用中的不足:

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

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

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

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

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

slide20
三阶Hermite插值

三次样条插值

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

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

区间 [xi-1, xi],

slide22
端点条件:
  • 固定斜率
  • 自然端点
  • 非节点

代入d3

另:三弯矩法

slide24
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.

slide25
MATLAB 函数

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

method = nearest

linear

spline

pchip

spline

pchip

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

slide27
yt = pchip(x,y,xt);

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

yt = spline(x,y,xt);

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

slide28
“It is important to prove,

but it is more

important to improve.”

ad