20.（1）0000.1)25.0(S ='，6868.0)53.0(S =' （2）0)53.0(S )25.0(S =''=''解：（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.67085.0x x 0143.1x 8863.1S 231++-= ]3.0,25.0[x ∈ 5477.0x 9127.0x 7314.0x 7952.0S 232++-= ]39.0,3.0[x ∈6245.0x 8004.0x 5167.0x 632.0S 233++-= ]45.0,39.0[x ∈ 6708.0x 7452.0x 4029.0x 3151.0S 234++-= ]53.0,45.0[x ∈ （2）在编辑窗口输入：>> x=[0.25,0.30,0.39,0.45,0.53];>> y=[0.5000,0.5477,0.6245,0.6708,0.7280]; >> f=ThrSample2(x,y,0,0,0) f =1000000/729*(5477/5000*t-279327/1000000)*(t-39/100)^2+1000000/729*(108663/200000-1249/1000*t)*(t-3/10)^2+10000/81*(8310710142007177/9007199254740992*t-24932130426021531/90071992547409920)*(39/100-t)^2-10000/81*(140377307168768943/450359962737049600-3599418132532537/4503599627370496*t)*(t-3/10)^21试用4次牛顿插值多项式)x (P 4及三次样条函数)x (S （自然边界条件）对数据进行插值。

用图给出{（i x ，i y ），10,11,1,0i ,i 08.02.0x i =+=}，)x (P 4及)x (S .解：（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 由此可以得出牛顿插值多项式为：)x (P 4=98.0x 1917.0x 1042.1x 8333.0x 5208.0234++-+- >> 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）在编辑窗口输入：>> 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),'-',x1,y1,'-.r'),grid0.20.30.40.50.60.70.80.91 1.1 1.2图2 P(x)和S(x)函数插值图．（1）用这9个点作8次多项式插值()xL8S x．（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 pointsor try centering and scaling as described in HELP POL YFIT.> In polyfit at 81m =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')010203040506070图3 Lagran插值与Spline插值的比较图从图中结果可以看出：在区间[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-1temp=a(k-1)/b(k-1);b(k)=b(k)-temp*c(k-1);u(k)=u(k)-temp*u(k-1);endm(n)=u(n-1)/b(n-1);for k=n-2:-1:1m(k+1)=(u(k)-c(k)*m(k+2))/b(k);endm(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-1S(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);elsedisp();return;endor i=1:nfor i=1:nif(x(i)<=x0)&&(x(i+1)>=x0)index=i;break;endendfif(x(i)<=x0)&&(x(i+1)>=x0)index=i;return;endendA=diag(2*ones(1,n));A(1,2)=1;A(n,n-1)=1;u=zeros(n-2,1);lamda=zeros(n-1,1);for i=2:n-1u(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);endc(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;endend;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);endd(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)/d(i-1,1))*b(i-1,1);endx(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(:,1)=y';for j=2:nfor k=j:nd(k,j)=(d(k,j-1)-d(k-1,j-1))/(x(k)-x(k-j+1)); endendc=d(n,n);for k=(n-1):-1:1c=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+1v=1;for j=1:n+1if k~=jv=conv(v,poly(x(j)))/(x(k)-x(j));endendl(k,:)=v;endc=y*l;。