插值法-第二次程序题

插值法
题目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（,）42
1
2cos(
⋅⋅⋅=+=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:0.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:0.001: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=[-0.6:0.001:0.6]; 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');
结果与分析：
由下图可以看出，在区间[-0.6,0.6]上，插值多项式可以很好的逼近被插值函数。

而在边界附近，插值多项式与被插值函数的差别很大。

即出现了Runge 现象。

主要原因是被插值函数的任意阶导数不能达到一致有界。

其插值余项
)()!
1()()（1)1(x n f x R n n n +++=
ϖξ不趋近零。

插值多项式不能收敛到被插值函数。

(2) 用节点）20，，2,1，0（,）42
1
2cos(
⋅⋅⋅=+=i i x i π，绘出它的20次Lagrange 插值多项式的图像。

Matlab 程序如下：
clear; %插值点
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));
Runge 函数 插值多项式
lx=lx*temp; end end l(i)=lx; end %插值多项式 l=l'; L=y*l;
%作图
a=[-1:0.01: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 。

Runge 函数 X
L 插值多项式
Runge 函数 X
L 插值多项式 Newton 插值多项式
(3)用等距节点，200,1.0,-1≤≤=+=i h ih x i 绘出它的分段线性插值函数
的图像。

Matlab 程序如下：
clc;clear; x=[-1:0.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
%作图
for i=1:n-1
a=[x(i):0.01:x(i+1)]; f=subs(l(i),z,a) plot(a,f,'k') hold on end
结果与分析：
如下图所示，分段线性插值多项式比较接近原函数，没有出现Runge 现象。

利用线性插值多项式的误差估计：
)2(2
max 8
h )（f x R n ≤
(4)用等距节点，200,1.0,-1≤≤=+=i h ih x i 绘出它的三次自然样条插值
函数的图像。

Matlab 程序如下：
clc;clear; x=[-1:0.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 h(i)=x(i+1)-x(i); end
for i=1:n-2
u(i)=h(i)/(h(i+1)+h(i)); r(i)=1-u(i); end
G=zeros(n,n); for i=1:n
……..
分段线性
Runge 函数
G(i,i)=2;
end
for i=2:n-1
G(i,i-1)=u(i-1);
G(i,i+1)=r(i-1);
end
G(n,n-1)=1;
G(1,2)=1;
d=zeros(1,n);
for i=2:n-1
d(i)=6*((y(i+1)-y(i))/h(i)-(y(i)-y(i-1))/h(i-1))/(h(i)+h(i-1)); end
syms u v;
u=diff(1/(1+25*v*v),v);
a=subs(u,v,x(1));
b=subs(u,v,x(n));
d(1)=((y(2)-y(1))/h(1)-a)/h(1)*6;
d(n)=(b-(y(n)-y(n-1))/h(n-1))/h(n-1)*6;
d=d';
M=inv(G)*d;
for i=1:n-1
s(i)=M(i)*(x(i+1)-z)^3/0.6+M(i+1)*(z-x(i))^3/0.6+(y(i)-M(i)*0.01/6)*( x(i+1)-z)/0.1+(y(i+1)-M(i+1)*0.01/6)*(z-x(i))/0.1;
end
for i=1:n-1
a=[x(i):0.01:x(i+1)];
f=subs(s(i),z,a);
plot(a,f,'x r')
hold on
end
结果与分析：
三次样条插值函数得到的图像如下：
可以看出，三次样条插值函数的曲线及其光滑。

得到的函数十分接近被插值函数。

XXXX Runge 函数
三次样条
题目2：对函数：
在区间[-1,1]作下列插值逼近，并和被插值函数的图像进行比较，并对结果进行分析。

(1)用等距节点，200,1.0,-1≤≤=+=i h ih x i 绘出它的20次Newton 插值
多项式的图像。

(2)用节点）20，，2,1，0（,）42
1
2cos(
⋅⋅⋅=+=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 程序如下：
clc;clear; %计算均差
x=[-1:0.1:1]; n=length(x); syms z ; y=zeros(1,n) for i=1:10
y(i)=sin(pi*x(i)); end
for i=11:15
y(i)=cos(pi*x(i)); end for i=15:n y(i)=0; 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
%作图
v=linspace(-1,0,50); u=sin(pi*v);
plot(v,u,'k')
hold on
v=linspace(0,0.5,25);
u=cos(pi*v);
plot(v,u,'k')
hold on
v=linspace(0.5,1,10000);
u=0;
plot(v,u,'k')
hold on
a=[-1:0.001:1];
fx=subs(f,z,a);
plot(a,fx,'r');
结果与分析：
等距节点20次Newton插值得到的函数图像如下：可以看出，在整个区间上，插值多项式精度都不是很高。

出现了Runge现象。

见下图
被插值函数
插值多项式
(2)用节点）20，，2,1，0（,）42
12cos(⋅⋅⋅=+=i i x i π，绘出它的20次Lagrange 插值多项式的图像。

Matlab 程序如下：
clc;clear;
%求插值节点
for i=1:21
x(i)=cos((2*(i-1)+1)*pi/42);
end
n=length(x);
y=zeros(1,n);
for i=1:n
if x(i)<0
y(i)=sin(pi*x(i));
elseif x(i)>0.5
y(i)=0;
else
y(i)=cos(pi*x(i));
end
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:0.01:1];
fx=subs(L,z,a);
plot(a,fx,'x r');
结果与分析：
如下图所示，使用Chebyshev多项式零点构造的Lagrange插值多项式比Newton插值多项式接近原函数，没有出现Runge现象。

被插值函数
插值多项式
(3)用等距节点，200,1.0,-1≤≤=+=i h ih x i 绘出它的分段线性插值函数
的图像。

Matlab 程序如下：
clc;clear;
x=[-1:0.1:1];
n=length(x);
syms z;
for i=1:10
y(i)=sin(pi*x(i));
end
for i=11:15
y(i)=cos(pi*x(i));
end
for i=15:n
y(i)=0;
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
%作图
for i=1:n-1
a=[x(i):0.01:x(i+1)];
f=subs(l(i),z,a);
plot(a,f,'x r')
hold on
end
结果与分析：
如下图所示，分段线性插值多项式比较接近原函数，没有出现Runge 现象。

但是在间断点处及导数不存在的点误差较大。

主要是因为这些地方构造的线性函数斜率较大，不能较好的趋近原函数。

(4)用等距节点，200,1.0,-1≤≤=+=i h ih x i 绘出它的三次自然样条插值
函数的图像。

Matlab 程序如下：
clc;clear;
x=[-1:0.1:1];
n=length(x);
syms z
for i=1:10
y(i)=sin(pi*x(i));
end
for i=11:15
y(i)=cos(pi*x(i));
end
for i=15:n
y(i)=0;
end
for i=1:n-1
h(i)=x(i+1)-x(i);
end
被插值函数
插值多项式
for i=1:n-2
u(i)=h(i)/(h(i+1)+h(i));
r(i)=1-u(i);
end
G=zeros(n,n);
for i=1:n
G(i,i)=2;
end
for i=2:n-1
G(i,i-1)=u(i-1);
G(i,i+1)=r(i-1);
end
G(n,n-1)=1;
G(1,2)=1;
d=zeros(1,n);
for i=2:n-1
d(i)=6*((y(i+1)-y(i))/h(i)-(y(i)-y(i-1))/h(i-1))/(h(i)+h(i-1)); end
syms u v;
u=diff(sin(pi*v),v);
a=subs(u,v,x(1));
b=0;
d(1)=((y(2)-y(1))/h(1)-a)/h(1)*6;
d(n)=(b-(y(n)-y(n-1))/h(n-1))/h(n-1)*6;
d=d';
M=inv(G)*d;
for i=1:n-1
s(i)=M(i)*(x(i+1)-z)^3/0.6+M(i+1)*(z-x(i))^3/0.6+(y(i)-M(i)*0.01/6)*( x(i+1)-z)/0.1+(y(i+1)-M(i+1)*0.01/6)*(z-x(i))/0.1;
end
for i=1:n-1
a=[x(i):0.01:x(i+1)];
f=subs(s(i),z,a);
plot(a,f,'x r')
hold on
end
结果与分析：
三次样条插值函数得到的图像如下：
可以看出，三次样条插值函数在间断点处也有较大误差。

被插值函数
插值多项式。