《插值法部分习题.doc》由会员分享，可在线阅读，更多相关《插值法部分习题.doc（7页珍藏版）》请在金锄头文库上搜索。
1、20.给定数据表如下：0.250.300.390.450.530.50000.54770.64250.67080.7280试求三次样条插值S（x），并满足条件：（1），（2）解：（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
2、 0.3151 -0.4029 0.7452 0.6708 （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-2493213042602153
3、1/90071992547409920)*(39/100-t)2-10000/81*(140377307168768943/450359962737049600-3599418132532537/4503599627370496*t)*(t-3/10)21、已知函数在下列各点的值0.2 0.4 0.6 0.8 1.00.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=newpol
4、y(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）i0 1 11 10 0.2 0.28 1.08 1.00.98 0.9596 0.2403 0.38在编辑窗口输入： x=0.2,0.28,0.4,0.6,0.8
5、,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
6、,-.r),grid3下列数据点的插值x01491625364964y012345678可以得到平方根函数的近似，在区间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
7、 scaling as described in HELP POLYFIT. 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
8、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);
9、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);endm(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);endm(1)=3*(d(1)-dx0)/h(1)-m(2)/2;m(n+1)=3*(dxn-d
10、(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);elsedisp(); return;endor i=1:n for i=1:
11、n if(x(i)=x0) index=i; break; endendf if(x(i)=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);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
12、+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
13、+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)
