插值法-第二次程序题

插值法

题目1：对Runge 函数2

2511)(x x R +=在区间[-1,1]作下列插值逼近，并和R(x)的图像进行比较，并对结果进行分析。

(1)用等距节点，200,1.0,-1≤≤=+=i h ih x i 绘出它的20次Newton

插值多项式的图像。

(2)用节点）20，，2,1，0（,）4212cos(

⋅⋅⋅=+=i i x i π，绘出它的20次Lagrange 插值多项式的图像。

(3)用等距节点，200,1.0,-1≤≤=+=i h ih x i 绘出它的分段线性插

值函数的图像。

(4)用等距节点，200,1.0,-1≤≤=+=i h ih x i 绘出它的三次自然样

条插值函数的图像。

程序及分析：

(1)用等距节点，200,1.0,-1≤≤=+=i h ih x i 绘出它的20次Newton

插值多项式的图像。

Matlab 程序如下：

%计算均差

x=[-1::1];

n=length(x);

syms z

for i=1:n

y(i)=1/(1+25*x(i)*x(i));

end

N=zeros(n,n);

N(:,1)=y';

for j=2:n

for k=j:n

N(k,j)=(N(k,j-1)-N(k-1,j-1))/(x(k)-x(k-j+1));

end

end

for t=1:n

c(t)=N(t,t)

end

%构造插值多项式

f=N(1,1);

for k=2:n

a=1;

for r=1:(k-1)

a=a*(z-x(r));

end

f=f+N(k,k)*a;

end

%作图

a=[-1::1];

n=length(a);

for i=1:n

b(i)=1/(1+25*a(i)*a(i)); end

fx=subs(f,z,a);

subplot(2,1,1);

plot(a,b,'k',a,fx,'r');

c=[::];

n=length(c);

for i=1:n

d(i)=1/(1+25*c(i)*c(i));

end

fx=subs(f,z,c);

subplot(2,1,2);

plot(c,d,'k',c,fx,'r');

结果与分析：

由下图可以看出，在区间[,]上，插值多项式可以很好的逼近被插值函数。而在边界附近，插值多项式与被插值函数的差别很大。即出现了Runge现象。

主要原因是被插值函数的任意阶导数不能达到一致有界。其插值余项)()!1()()（1)1(x n f x R n n n +++=ϖξ不趋近零。插值多项式不能收敛到被插值函数。

(2) 用节点）20，，2,1，0（,）4212cos(

⋅⋅⋅=+=i i x i π，绘出它的20次Lagrange 插值多项式的图像。

Matlab 程序如下：

clear;

%插值点

Runge 函数

插值多项

式

for i=1:21

x(i)=cos((2*(i-1)+1)*pi/42);

end

n=length(x);

for i=1:n

y(i)=1/(1+25*x(i)*x(i));

end

%构造插值基函数

syms z;

temp=1;

for i=1:n

lx=1;

for j=1:n

if i~=j

temp=(z-x(j))/(x(i)-x(j));

lx=lx*temp;

end

end

l(i)=lx;

end

%插值多项式

l=l';

L=y*l;

%作图

a=[-1::1];

n=length(a);

for i=1:n

b(i)=1/(1+25*a(i)*a(i)); end

fx=subs(L,z,a);

subplot(2,1,1);

plot(a,b,'k',a,fx,'x r');

结果与分析：

如下图所示，使用Chebyshev 多项式零点构造的Lagrange 插值多项式比较接近原函数，没有出现Runge 现象。 主要原因是其多项式误差为∞++≤)1()!1(21)（-)(n n n f n x L x f 。

(3)用等距节点，200,1.0,-1≤≤=+=i h ih x i 绘出它的分段线性插

值函数的图像。

Matlab程序如下：

clc;clear;

x=[-1::1];

n=length(x);

syms z

for i=1:n

y(i)=1/(1+25*x(i)*x(i));

end

%构造分段线性插值多项式

for i=1:n-1

l(i)=(z-x(i+1))/(x(i)-x(i+1))*y(i)+(z-x(i))/(x(i+1)-x(i))*y (i+1)

% l(i)=y(i)+(y(i+1)-y(i))/(x(i+1)-x(i))*(z-x(i))

end