Interpolation

1 / 28

# Interpolation - PowerPoint PPT Presentation

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;

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

## PowerPoint Slideshow about 'Interpolation' - hisa

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

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

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

}

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

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

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

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

Bezier 曲线在应用中的不足：

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

• 固定斜率
• 自然端点
• 非节点

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.

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