数值分析实验3计算142许莎莎-2014014517

数值分析实验三专业班级：信计142姓名：许莎莎学号：20140145171）实验目的掌握三种插值方法：牛顿多项式插值，三次样条插值，拉格朗日插值，并比较三种插值方法的优劣。编写牛顿多项式插值，三次样条插值，拉格朗日插值的程序编码，并在MATLAB软件中去实现2）实验题目、内容、结果1.已知函数在下列各点的值为xi0.20.40.60.81.0f(xi)0.980.920.810.640.38试用4次牛顿插值多项式P4（x）及三次样条函数S（x）（自然边界条件）对数据进行插值，用图给出{(xi,yi),xi=0.2+0.08i,i=0,1,11,10},P4(x),及S(x)。牛顿插值多项式functionvarargout=newtonliu(varargin)clear,clcx=[0.20.40.60.81.0];fx=[0.980.920.810.640.38];newtonchzh(x,fx);functionnewtonchzh(x,fx)n=length(x);fprintf('*****************²î·Ö±í*****************************\n');FF=ones(n,n);FF(:,1)=fx';fori=2:nforj=i:nFF(j,i)=(FF(j,i-1)-FF(j-1,i-1))/(x(j)-x(j-i+1));endendfori=1:nfprintf('%4.2f',x(i));forj=1:ifprintf('%10.5f',FF(i,j));endfprintf('\n');end计算结果为：xif(xi)一阶差商二阶差商三阶差商四阶差商0.200.9800.400.920-0.300000.600.810-0.55000-0.625000.800.640-0.85000-0.75000-0.208331.000.380-1.30000-1.12500-0.62500-0.52083所以有四次插值牛顿多项式为：P4（x）=0.98-0.3(x-0.2)-0.62500(x-0.2)(x-0.4)-0.20833(x-0.2)(x-0.4)(x-0.6)-0.52083(x-0.2)(x-0.4)(x-0.6)(x-0.8）三次样条插值用三次样条插值函数由上题分析知,要求各点的M值：06.7500-4.5000-3.7500-0MMMMM2.5000000020.50000000.50002.500000.50002.5000000243210functiontgsanci(n,s,t)x=[0.20.40.60.81.0];y=[0.980.920.810.640.38];n=5forj=1:1:n-1h(j)=x(j+1)-x(j);endforj=2:1:n-1r(j)=h(j)/(h(j)+h(j-1));endforj=1:1:n-1u(j)=1-r(j);endforj=1:1:n-1f(j)=(y(j+1)-y(j))/h(j);endforj=2:1:n-1d(j)=6*(f(j)-f(j-1))/(h(j-1)+h(j));endd(1)=0d(n)=0a=zeros(n,n);forj=1:1:na(j,j)=2;endr(1)=0;u(n)=0;forj=1:1:n-1a(j+1,j)=u(j+1);a(j,j+1)=r(j);endb=inv(a)m=b*d'p=zeros(n-1,4);forj=1:1:n-1p(j,1)=m(j)/(6*h(j));p(j,2)=m(j+1)/(6*h(j));p(j,3)=(y(j)-m(j)*(h(j)^2/6))/h(j);p(j,4)=(y(j+1)-m(j+1)*(h(j)^2/6))/h(j);endp得到m=(0-1.6071-1.0714-3.10710)T即M0=0;M1=-1.6071;M2=-1.0714;M3=-3.1071;M4=0则根据三次样条函数定义，可得：S(x)=]0.1,8.0[x)8.0(9.10.13036.38.00-x0.12.5893-]8.0,6.0[x6.0x3036.3x8.00857.4.60-x5893.2-x.800.8929-]6.0,4.0[x4.0x7508.4x6.04.6536.40-x8929.0x.601.3393-]4.0,2.0[x)2.0(6536.44.0900.42.01.33934.003333333，）（）（）（），（）（）（）（），（）（）（）（，）（）（）（xxxxxxx画牛顿插值以及三次样条插值图形x=[0.20.40.60.81.0];y=[0.980.920.810.640.38];plot(x,y)holdonfori=1:1:5y(i)=0.98-0.3*(x(i)-0.2)-0.62500*(x(i)-0.2)*(x(i)-0.4)-0.20833*(x(i)-0.2)*(x(i)-0.4)*(x(i)-0.6)-0.52083*(x(i)-0.2)*(x(i)-0.4)*(x(i)-0.6)*(x(i)-0.8)endk=[011011]x0=0.2+0.08*kfori=1:1:4y0(i)=0.98-0.3*(x(i)-0.2)-0.62500*(x(i)-0.2)*(x(i)-0.4)-0.20833*(x(i)-0.2)*(x(i)-0.4)*(x(i)-0.6)-0.52083*(x(i)-0.2)*(x(i)-0.4)*(x(i)-0.6)*(x(i)-0.8)endplot(x0,y0,'o',x0,y0)holdony1=spline(x,y,x0)plot(x0,y1,'o')holdons=csape(x,y,'variational')fnplt(s,'r')holdongtext('三次样条自然边界')gtext('原图像')gtext('4次牛顿插值')2.在区间[-1,1]上分别取难0、20用两组等距节点对龙格函数f(x)=211+25x作多项式插值及三样条插值，对每个n值，分别画出插值函数及f(x)的图形。建立M文件多项式插值function[C,D]=newpoly(X,Y)n=length(X);D=zeros(n,n)D(:,1)=Y'forj=2:nfork=j:nD(k,j)=(D(k,j-1)-D(k-1,j-1))/(X(k)-X(k-j+1));endendC=D(n,n);fork=(n-1):-1:1C=conv(C,poly(X(k)))m=length(C);C(m)=C(m)+D(k,k);end当n=20时，在CommandWindow中输入以下的命令：clear,clf,holdon;X=-1:0.2:1;Y=1./(1+25*X.^2);[C,D]=newpoly(X,Y);x=-1:0.01:1;y=polyval(C,x);plot(x,y,X,Y,'.');gridon;xp=-1:0.2:1;z=1./(1+25*xp.^2);plot(xp,z,'r')当n=20时，在CommandWindow中输入以下的命令：clear,clf,holdon;X=-1:0.1:1;Y=1./(1+25*X.^2);[C,D]=newpoly(X,Y);x=-1:0.01:1;y=polyval(C,x);plot(x,y,X,Y,'.');gridon;xp=-1:0.1:1;z=1./(1+25*xp.^2);plot(xp,z,'r')三样条插值当n=10时，在CommandWindow中输入以下的命令：clear,clcX=-1:0.2:1;Y=1./(25*X.^2+1);dx0=0.0739644970414201;dxn=-0.0739644970414201;S=csfit(X,Y,dx0,dxn)x1=-1:0.01:-0.5;y1=polyval(S(1,:),x1-X(1));x2=-0.5:0.01:0;y2=polyval(S(2,:),x2-X(2));x3=0:0.01:0.5;y3=polyval(S(3,:),x3-X(3));x4=0.5:0.01:1;y4=polyval(S(4,:),x4-X(4));plot(x1,y1,x2,y2,x3,y3,x4,y4,X,Y,'.')当n=20时，在CommandWindow中输入以下的命令：clear,clcX=-1:0.1:1;Y=1./(25*X.^2+1);dx0=0.0739644970414201;dxn=-0.0739644970414201;S=csfit(X,Y,dx0,dxn)x1=-1:0.01:-0.5;y1=polyval(S(1,:),x1-X(1));x2=-0.5:0.01:0;y2=polyval(S(2,:),x2-X(2));x3=0:0.01:0.5;y3=polyval(S(3,:),x3-X(3));x4=0.5:0.01:1;y4=polyval(S(4,:),x4-X(4));plot(x1,y1,x2,y2,x3,y3,x4,y4,X,Y,'.')3.下列数据点的插值X01491625364964Y012345678可以得到平方根数的近似，在区间[0,64]上作图。(1)用这9个点作8次多项式插值L8(x).(2)用三次样条(第一边界条件)程序求S(x)。从得到结果看在[0,64]上，哪个插值更精确；在区间[0,1]上，两种插值哪个更精确？（1）8次多项式插值x1=[01491625364964];y1=[012345678];n=length(y1);a=ones(n,2);a(:,2)=-x1';c=1;fori=1:nc=conv(c,a(i,:));endq=zeros(n,n);r=zeros(n,n+1);fori=1:n[q(i,:),r(i,:)]=deconv(c,a(i,:))endDw=zeros(1,n);fori=1:nDw(i)=y1(i)/polyval(q(i,:),x1(i));endp=Dw*q;symsxL8;fori=1:nL8(i)=p(n-i+1)*x^(i-1);enddisp('');L8=vpa(collect((sum(L8))),5)L8=-0.00000000032806*x^8+0.000000067127*x^7-0.0000054292*x^6+0.00022297*x^5-0.0049807*x^4+0.060429*x^3-0.38141*x^2+1.3257*x（2）三次样条插值functiontgsanci(n,s,t)y=[012345678];x=[01491625364964];n=9forj=1:1:n-1h(j)=x(j+1)-x(j);endforj=2:1:n-1r(j)=h(j)/(h(j)+h(j-1));endforj=1:1:n-1u(j)=1-r(j);endforj=1:1:n-1f(j)=(y(j+1)-y(j))/h(j);endforj=2:1:n-1d(j)=6*(f(j)-f(j-1))/(h(j-1)+h(j));endd(1)=0d(n)=0a=zeros(n,n);forj=1:1:na(j,j)=2;endr(1)=0;u(n)=0;forj=1:1:n-1a(j+1,j)=u(j+1);a(j,j+1)=r(j);endb=inv(a)m=b*d't=ap=zeros(n-1,4);forj=1:1:n-1p(j,1)=m(j)/(6*h(j));p(j,2)=m(j+1)/(6*h(j));p(j,3)=(y(j)-m