第五章任玉杰作业

第五章任玉杰作业
6.2 第六题
解：先保存名为lagranzi.m的M文件如下：
function [y,R]=lagranzi(X,Y,x,M)
n=length(X); m=length(x);
for i=1:m
z=x(i);s=0.0;
for k=1:n
p=1.0; q1=1.0; c1=1.0;
for j=1:n
if j~=k
p=p*(z-X(j))/(X(k)-X(j));
end
q1=abs(q1*(z-X(j)));c1=c1*j;
end
s=p*Y(k)+s;
end
y(i)=s;
end
R=M*q1/c1;
在matlab里输入主程序：
x=2*pi/9; M=1; X=[pi/6 ,pi/4, pi/3];Y=[0.5,0.7071,0.8660]; [y,R]=lagranzi(X,Y,x,M)
运行得结果：
y =

0.6434

R =
8.8610e-004
6.3 第四题
解：先保存牛顿插值调用函数：
function [y,R,A,C,L]=newdscg(X,Y,x,M)
n=length(X); m=length(x);
for t=1:m
z=x(t); A=zeros(n,n);A(:,1)=Y';
s=0.0; p=1.0; q1=1.0; c1=1.0;
for j=2:n
for i=j:n
A(i,j)=(A(i,j-1)- A(i-1,j-1))/(X(i)-X(i-j+1));
end
q1=abs(q1*(z-X(j-1)));c1=c1*j;
end
C=A(n,n);q1=abs(q1*(z-X(n)));
for k=(n-1):-1:1
C=conv(C,poly(X(k)));
d=length(C);C(d)=C(d)+A(k,k);
end
y(k)= polyval(C, z);
end
R=M*q1/c1;L(k,:)=poly2sym(C);
再输入主程序：
syms M,X=[-3,0,1,2]; Y =[27,1,2,17]; x=1.4;
[y,R,P]=newdscg(X,Y,x,M)
运行得结果：
y =

6.0120

R =
77/1250*M

P =
27.0000 0 0 0
1.0000 -8.6667 0 0
2.0000 1.0000 2.4167 0
17.0000 15.0000 7.0000 0.9167

6.4 第五题
给出节点数据(3.15)37.03f，(1.00)7.24f，(0.01)1.05f，03.215.3'f，
06.1700.1'f,
05.2301.0'f
作埃尔米特插值多项式，并写出估计其误差的公式。

解：
保存名为hermite.m的M文件：
function [Hc, Hk,wcgs,Cw]= hermite (X,Y,Y1,M1)
m=length(X); n=m-1;s=0; H=0;q=1;c1=1; L=ones(m,m); G=ones(1,m);
for k=1:n+1
V=1;
for i=1:n+1
if k~=i
s=s+(1/(X(k)-X(i)));
V=conv(V,poly(X(i)))/(X(k)-X(i));
end
h=poly(X(k));g=(1-2*h*s);G=g*Y(k)+h*Y1(k);
end
H=H+conv(G,conv(V,V));b=poly(X(k));b2=conv(b,b);
q=conv(q,b2);t=2*n+2;
Hc=H;Hk=poly2sym(H); Q=poly2sym(q);
end
for i=1:t
c1=c1*i;
end
syms M,wcgs=M*Q/c1; Cw=q/c1;

在MATLAB工作窗口输入程序
>>X=[-3.15,-1.00,0.01];Y=[37.03,7.24,1.05];Y1=[2.03,17.06,23.05];
[Hc,Hk,wcgs,Cw]=hermite(X,Y,Y1)
运行后输出埃尔米特插值多项式Hk及其系数向量Hc，误差公式wcgs及其系数向量Cw如下
Hc =
1.0e+02 *
Columns 1 through 5
0.136316058211219 0.954914816902903 2.153699318501380 1.563807894332014
0.236785931159195
Column 6
0.008078596636877
Hk =
(1918477965514565*x^5)/140737488355328 + (6719615646210121*x^4)/70368744177664 +
(1894413954740411*x^3)/8796093022208 + (5502159882963045*x^2)/35184372088832 +
(6664931445844529*x)/281474976710656 + 7276552960702781/9007199254740992
wcgs =
(M*(x^6 + (207*x^5)/25 + (1643574610260027*x^4)/70368744177664 +
(7226976987537243*x^3)/281474976710656 + (5292828581206899*x^2)/562949953421312 -
(3527858739303659*x)/18014398509481984 + 35749573842067/36028797018963968))/720
Cw =
Columns 1 through 5
0.001388888888889 0.011500000000000 0.032439722222222 0.035660250000000
0.013058267013889
Columns 6 through 7
-0.000271993750000 0.000001378125000

其上机如下图：
6.5 第二题
解：先保存名为xxczhjt1的M文件如下：
function s= xxczhjt1(x0,y0,xi,x,y)
s= interp1(x0,y0,xi);
Sn= interp1(x0,y0,x0);
plot(x0,y0,'o',x0,Sn,'-',xi,s,'*',x,y,'-.')
legend('节点(xi,yi)','分段线性插值函数Sn(x)','插值点(x,s)','被插值函数y')
再输入主程序如下：
>> h=2*pi/13; x0=-pi:h:pi;
y0=sin(3-cos(4*x0)./(1+25*x0.^2));
xi=-pi+h/2:h:pi-h/2;
fi= sin(3-cos(4*xi)./(1+25*xi.^2));
b=max(x0); a=min(x0);
x=a:0.001:b; y= sin(3-cos(4*x)./(1+25*x.^2));
si=xxczhjt1(x0,y0,xi,x,y);
Ri=abs((fi-si)./fi);
i=[xi',fi',si',Ri']
title('y= sin(3-cos(4 x)./(1+25 x^2))的分段线性插值的有关图形')
运行后得到结果：
i =

-2.8999 0.1438 0.1421 0.0115
-2.4166 0.1346 0.1370 0.0182
-1.9333 0.1424 0.1440 0.0117
-1.4500 0.1575 0.1487 0.0554
-0.9666 0.1106 0.1086 0.0186
-0.4833 0.0896 0.2184 1.4368
0 0.9093 0.3640 0.5997
0.4833 0.0896 0.2184 1.4368
0.9666 0.1106 0.1086 0.0186
1.4500 0.1575 0.1487 0.0554
1.9333 0.1424 0.1440 0.0117
2.4166 0.1346 0.1370 0.0182
2.8999 0.1438 0.1421 0.0115