插值法部分习题.doc

20.给定数据表如下：

0.25

0.30

0.39

0.45

0.53

0.5000
0.5477
0.6425
0.6708
0.7280
试求三次样条插值S（x），并满足条件：
解：（1）在编辑窗口输入：
>> x=[0.25,0.30,0.39,0.45,0.53];
>> y=[0.5000,0.5477,0.6245,0.6708,0.7280];
>> dx0=1.0000;dxn=0.6868;
>> s=csfit(x,y,dx0,dxn)
s =
1.8863 -1.0143 1.0000 0.5000
0.7952 -0.7314 0.9127 0.5477
0.6320 -0.5167 0.8004 0.6245
0.3151 -0.4029 0.7452 0.6708

1、已知函数在下列各点的值

0.2 0.4 0.6 0.8 1.0

0.98 0.92 0.81 0.64 0.38
试用4次牛顿插值多项式及三次样条函数（自然边界条件）对数据进行插值。用图给出{（，），}，及.
解：
（1）在编辑窗口输入：
>> x=[0.2,0.4,0.6,0.8,1.0];
>> y=[0.98,0.92,0.81,0.64,0.38];
>> N=newpoly(x,y)
N =
-0.5208 0.8333 -1.1042 0.1917 0.9800
由此可以得出牛顿插值多项式为：
=
>> S=ThrSample2(x,y,0,0,0)
S =
125*(46/25*t-69/125)*(t-3/5)^2+125*(567/500-81/50*t)*(t-2/5)^2+25*(-57/140*t+57/350)*(3/5-t)^2-25*(-81/200+27/40*t)*(t-2/5)^2/5)^2
（2）

i
0 1 11 10

0.2 0.28 1.08 1.0

0.98 0.9596 0.2403 0.38

在编辑窗口输入：
>> x=[0.2,0.28,0.4,0.6,0.8,1.0,1.08];
>> y=[0.98 0.9596 0.92 0.81 0.64 0.38 0.2403];
>> cs=spline(x,[0 y 0]);xx=linspace(0.2,1.08,101);
>> plot(x,y,'o',xx,ppval(cs,xx),'-');
>> x1=[0.2,0.28,0.4,0.6,0.8,1.0,1.08];
>> y1=[0.98 0.9596 0.92 0.81 0.64 0.38 0.2403];
>> polyfit(x1,y1,4)
ans =
-0.4126 0.5487 -0.8475 0.1017 0.9893
>> plot(x,y,'o',xx,ppval(cs,xx),'-',x1y1,'-.r'),grid

3．下列数据点的插值
x
0
1
4
9
16
25
36
49
64
y
0
1
2
3
4
5
6
7
8
可以得到平方根函数的近似，在区间[0,64]上作图．
（1）用这9个点作8次多项式插值．
（2）用三次样条（自然边界条件）程序求．
解：
在编辑窗口输入：
>> x=[0 1 4 9 16 25 36 49 64];
>> y=0:1:8;y=sqrt(x);
>> x1=[0 1 4 9 16 25 36 49 64];y1=0:1:8;
>> m=polyfit(x,y,8)
Warning: Polynomial is badly conditioned. Remove repeated data points
or try centering and scaling as described in HELP POLYFIT.
> In polyfit at 81
m =
Columns 1 through 4
-0.00000000032806 0.00000006712680 -0.00000542920942 0.00022297151886
Columns 5 through 8
-0.00498070758014 0.06042943489173 -0.38141003716579 1.32574370075005
Column 9
-0.00000000000311
>> x1=0:1:64;
>> y1=polyval(m,x1);
>> x2=[0 1 4 9 16 25 36 49 64];
>> y2=0:1:8;
>> cs=spline(x,[0 y 0]);
>> xx=linspace(0,64,101);
>> plot(x,y,'o',xx,ppval(cs,xx),'-k',x1,y1,'--b',x,y,':r')

从图中结果可以看出：在区间[0,64]上，三次样条插值更准确些；在区间[0,1]上，拉格朗日插值更准确些.
附：各插值函数的Matlab实现：
(1) csfit函数运行程序
function S=csfit(x,y,dx0,dxn)
n=length(x)-1;
h=diff(x);
d=diff(y)./h;
a=h(2:n-1);
b=2*(h(1:n-1)+h(2:n));
c=h(2:n);
u=6*diff(d);
b(1)=b(1)-h(1)/2;
u(1)=u(1)-3*(d(1)-dx0);
b(n-1)=b(n-1)-h(n)/2;
u(n-1)=u(n-1)-3*(dxn-d(n));
for k=2:n-1
temp=a(k-1)/b(k-1);
b(k)=b(k)-temp*c(k-1);
u(k)=u(k)-temp*u(k-1);
end
m(n)=u(n-1)/b(n-1);
for k=n-2:-1:1
m(k+1)=(u(k)-c(k)*m(k+2))/b(k);
end
m(1)=3*(d(1)-dx0)/h(1)-m(2)/2;
m(n+1)=3*(dxn-d(n))/h(n)-m(n)/2;
for k=0:n-1
S(k+1,1)=(m(k+2)-m(k+1))/(6*h(k+1));
S(k+1,2)=m(k+1)/2;
S(k+1,3)=d(k+1)-h(k+1)*(2*m(k+1)+m(k+2))/6;
S(k+1,4)=y(k+1);
End
(2)ThrSample2函数运行程序
function [f,f0]=ThrSample2(x,y,y2_1,y2_N,x0)
syms t;
f=0.0;
f0=0.0;
if(length(x)==length(y))
n=length(x);
else
disp();
return;
endor i=1:n

for i=1:n
if(x(i)<=x0)&&(x(i+1)>=x0)
index=i;
break;
end
end
f
if(x(i)<=x0)&&(x(i+1)>=x0)
index=i;
return;
end
end
A=diag(2*ones(1,n));
A(1,2)=1;
A(n,n-1)=1;
u=zeros(n-2,1);
lamda=zeros(n-1,1);
c=zeros(n,1);
for i=2:n-1
u(i-1)=(x(i)-x(i-1))/(x(i+1)-x(i-1));
lamda(i)=(x(i+1)-x(i))/(x(i+1)-x(i-1));
c(i)=3*lamda(i)*(y(i)-y(i-1))/(x(i)-x(i-1))+3*u(i-1)*(y(i+1)-y(i))/(x(i+1)-x(i));
A(i,i+1)=u(i-1);
A(i,i-1)=lamda(i);
end
c(1)=3*(y(2)-y(1))/(x(2)-x(1))-(x(2)-x(1))*y2_1/2;
c(n)=3*(y(n)-y(n-1))/(x(n)-x(n-1))x(n)-x(n-1))*y2_N/2;
m=followup(A,c);
h=x(i+1)-x(i);
f=y(i)*(2*(t-x(i))+h)*(t-x(i+1))^2/h/h/h+y(i+1)*(2*(x(i+1)-t)+h)*(t-x(i))^2/h/h/h+m(i)*(t-x(i))*(x(i+1)-t)^2/h/h-m(i+1)*(x(i+1)-t)*(t-x(i))^2/h/h;
f0=subs(f,'t',x0);
%ThrSample2函数运行调用的函数：
function x=followup(A,b)
n=rank(A);
for(i=1:n)
if(A(i,i)==0)
disp();
return;
end
end;
d=ones(n,1);
a=ones(n-1,1);
c=ones(n-1);
for(i=1:n-1)
a(i,1)=A(i+1,i);
c(i,1)=A(i,i+1);
d(i,1)=A(i,i);
end
d(n,1)=A(n,n);
for(i=2:n)
d(i,1)=d(i,1)-(a(i-1,1)/d(i-1,1))*c(i-1,1);
b(i,1)=b(i,1)-(a(i-1,1)/di-1,1))*b(i-1,1);
end
x(n,1)=b(n,1)/d(n,1);
for(i=(n-1):-1:1)
x(i,1)=(b(i,1)-c(i,1)*x(i+1,1))/d(i,1);
end
(3) newpoly函数运行程序
function[c,d]=newpoly(x,y)
n=length(x);
d=zeros(n,n);
d(:,1)=y';
for j=2:n
for k=j:n
d(k,j)=(d(k,j-1)-d(k-1,j-1))/(x(k)-x(k-j+1));
end
end
c=d(n,n);
for k=(n-1):-1:1
c=conv(c,poly(x(k)));
m=length(c);
c(m)=c(m)+d(k,k);
end
(4) lagran函数运行程序
function[c,l]=lagran(x,y)
w=length(x);
n=w-1;
l=zeros(w,w);
for k=1:n+1
v=1;
for j=1:n+1
if k~=j
v=conv(v,poly(x(j)))/(x(k)-x(j));
end
end
l(k,:)=v;
end
c=y*l;
6