数值分析实验报告最新版

实验项目名称 插值法 实 验 室 数学实验室 所属课程名称 数值逼近 实 验 类 型 算法设计 实 验 日 期

班 级 学 号 姓 名 成 绩

可编辑word,供参考版！

实验概述： 【实验目的及要求】 本次实验的目的是熟练《数值分析》第二章“插值法”的相关内容，掌握三种插值方法：牛顿多项式插值，三次样条插值，拉格朗日插值，并比较三种插值方法的优劣。 本次试验要求编写牛顿多项式插值，三次样条插值，拉格朗日插值的程序编码，并在MATLAB软件中去实现。 【实验原理】 《数值分析》第二章“插值法”的相关内容，包括：牛顿多项式插值，三次样条插值，拉格朗日插值的相应算法和相关性质。 【实验环境】（使用的软硬件） 软件： MATLAB 2012a 硬件： 电脑型号：联想 Lenovo 昭阳E46A笔记本电脑 操作系统：Windows 8 专业版 处理器：Intel（R）Core（TM）i3 CPU M 350 @2.27GHz 2.27GHz 实验内容： 【实验方案设计】 第一步，将书上关于三种插值方法的内容转化成程序语言，用MATLAB实现；第二步，分别用牛顿多项式插值，三次样条插值，拉格朗日插值求解不同的问题。 【实验过程】（实验步骤、记录、数据、分析） 实验的主要步骤是：首先分析问题，根据分析设计MATLAB程序，利用程序算出问题答案，分析所得答案结果，再得出最后结论。 实验一： 已知函数在下列各点的值为 xi 0.2 0.4 0.6 .0.8 1.0 f（xi） 0.98 0.92 0.81 0. 0.38 试用4次牛顿插值多项式P4（x）及三次样条函数S（x）（自然边界条件）对数据进行插值。用图给出{（xi，yi），xi=0.2+0.08i，i=0，1, 11, 10}，P4（x）及S（x）。 (1)首先我们先求牛顿插值多项式，此处要用4次牛顿插值多项式处理数据。 已知n次牛顿插值多项式如下： Pn=f(x0)+f[x0,x1](x-x0)+ f[x0,x1,x2](x-x0) (x-x1)+···+ f[x0,x1，···xn](x-x0) ···(x-xn-1) 可编辑word,供参考版！

我们要知道牛顿插值多项式的系数，即均差表中得部分均差。 在MATLAB的Editor中输入程序代码，计算牛顿插值中多项式系数的程序如下： function varargout=newtonliu(varargin) clear,clc x=[0.2 0.4 0.6 0.8 1.0]; fx=[0.98 0.92 0.81 0. 0.38]; newtonchzh(x,fx); function newtonchzh(x,fx) %由此函数可得差分表 n=length(x); fprintf('*****************差分表*****************************\\n'); FF=ones(n,n); FF(:,1)=fx'; for i=2:n for j=i:n FF(j,i)=(FF(j,i-1)-FF(j-1,i-1))/(x(j)-x(j-i+1)); end end for i=1:n fprintf('%4.2f',x(i)); for j=1:i fprintf('%10.5f',FF(i,j)); end fprintf('\\n'); end 由MATLAB计算得： xi f(xi) 一阶差商 0.20 0.980 0.40 0.920 -0.30000 0.60 0.810 -0.55000 0.80 0.0 -0.85000 1.00 0.380 -1.30000 所以有四次插值牛顿多项式为： 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 二阶差商 -0.62500 -0.75000 -1.12500 三阶差商 -0.20833 -0.62500 四阶差商 -0.52083 可编辑word,供参考版！

(x-0.2)(x-0.4)(x-0.6)(x-0.8） (2)接下来我们求三次样条插值函数。 用三次样条插值函数由上题分析知,要求各点的M值： 0000M0020.500M-3.750020.50000100.50020.5000M2-4.5000 000.50020.500M3-6.7500000020M4三次样条插值函数计算的程序如下： function tgsanci(n,s,t) %n代表元素数，s,t代表端点的一阶导。 x=[0.2 0.4 0.6 0.8 1.0]; y=[0.98 0.92 0.81 0. 0.38]; n=5 for j=1:1:n-1 h(j)=x(j+1)-x(j); end for j=2:1:n-1 r(j)=h(j)/(h(j)+h(j-1)); end for j=1:1:n-1 u(j)=1-r(j); end for j=1:1:n-1 f(j)=(y(j+1)-y(j))/h(j); end for j=2:1:n-1 d(j)=6*(f(j)-f(j-1))/(h(j-1)+h(j)); end d(1)=0 d(n)=0 a=zeros(n,n); for j=1:1:n a(j,j)=2; end 可编辑word,供参考版！

r(1)=0; u(n)=0; for j=1:1:n-1 a(j+1,j)=u(j+1); a(j,j+1)=r(j); end b=inv(a) m=b*d' p=zeros(n-1,4); %p矩阵为S(x)函数的系数矩阵 for j=1:1:n-1 p(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); end p 得到m=(0 -1.6071 -1.0714 -3.1071 0) 即M0=0 ;M1= -1.6071;M2= -1.0714; M3= -3.1071; M4=0 则根据三次样条函数定义，可得： 3 （00.4x）1.3393（ x0.2）4.900（0.4x）4.6536(x0.2)，x[0.2,0.4]33 -1.3393（0.6x）0.29（ x-0.4） 4.6536（0.6x）4.085（7x0.4），x[0.4,0.6] S(x)= 33（0.8x）-2.53（x-0.6）4.0857（0.8x）3.3036（x0.6），x[0.6,0.8]-0.2933-2.53（1.0x）-（0x0.8）3.3036（ 1.0x）1.9(x0.8)，x[0.8,1.0]T 接着，在Command Window里输入画图的程序代码， 下面是画牛顿插值以及三次样条插值图形的程序： x=[0.2 0.4 0.6 0.8 1.0]; y=[0.98 0.92 0.81 0. 0.38]; plot(x,y) hold on for i=1:1:5 y(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) end k=[0 1 10 11] 可编辑word,供参考版！

x0=0.2+0.08*k for i=1:1:4 y0(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) end plot( x0,y0,'o',x0,y0 ) hold on y1=spline(x,y,x0) plot(x0,y1,'o') hold on s=csape(x,y,'variational') fnplt(s,'r') hold on gtext('三次样条自然边界') gtext('原图像') gtext('4次牛顿插值') 运行上述程序可知：给出的{（xi，yi），xi=0.2+0.08i，i=0，1, 11, 10}点，S（x）及P4（x）图形如下所示： 10.94次牛顿插值0.80.70.60.50.40.30.20.2三次样条自然边界原图像0.30.40.50.60.70.80.911.11.2 实验二： 在区间[-1,1]上分别取n10,20用两组等距节点对龙格函数f(x)1作多项式插值及三次样条插125x2可编辑word,供参考版！

值，对每个n值，分别画出插值函数即f(x)的图形。 我们先求多项式插值： 在MATLAB的Editor中建立一个多项式的M-file,输入如下的命令（如牛顿插值公式）： 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 当n=10时，我们在Command Window中输入以下的命令： clear,clf,hold on; 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,'.'); grid on; xp=-1:0.2:1; z=1./(1+25*xp.^2); plot(xp,z,'r') 得到插值函数和f（x）图形： 可编辑word,供参考版！

当n=20时，我们在Command Window中输入以下的命令： clear,clf,hold on; 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,'.'); grid on; xp=-1:0.1:1; z=1./(1+25*xp.^2); plot(xp,z,'r') 得到插值函数和f（x）图形： 可编辑word,供参考版！

下面再求三次样条插值函数，在MATLAB的Editor中建立一个多项式的M-file, 输入下列程序代码： 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)); B(N-1)=B(N-1)-H(N)/2; U(N-1)=U(N-1)-3*(-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 可编辑word,供参考版！

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 当n=10时，我们在Command Window中输入以下的命令： clear,clc X=-1:0.2:1; Y=1./(25*X.^2+1); dx0= 0.07394970414201;dxn= -0.07394970414201; 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时，我们在Command Window中输入以下的命令： 可编辑word,供参考版！

clear,clc X=-1:0.1:1; Y=1./(25*X.^2+1); dx0= 0.07394970414201;dxn= -0.07394970414201; 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,'.') 结果如图： 实验三： 下列数据点的插值 x 0 1 4 9 16 25 36 49 y 0 1 2 3 4 5 6 7 8 可以得到平方根函数的近似，在区间[0,]上作图。 （1）用这9各点作8次多项式插值L8(x). （2）用三次样条（自然边界条件）程序求S（x）。从结果看在[0,]上，那个插值更精确；在区间[0,1]上，两种哪个更精确？ L8(x)可由公式Ln(x)=yk0nkkl(xj)得出。 三次样条可以利用自然边界条件。写成矩阵： 可编辑word,供参考版！

21000其中j= 020122003000M0d0Md001120M2d2 23M3d342M4d4hj-1hj-1hj，i=hjhj-1hj，dj=6f[xj-1,xj,xj+1]，n=0=0 d0=dn=0 l0(x)=(x-1)(x4)(x9)(x16)(x25)(x36)(x49)(x) (01)(04)(09)(016)(025)(036)(049)(0)(x-0)(x4)(x9)(x16)(x25)(x36)(x49)(x) (10)(14)(19)(116)(125)(136)(149)(1)(x-0)(x1)(x9)(x16)(x25)(x36)(x49)(x)l2(x)= (40)(41)(49)(416)(425)(436)(449)(4)(x-0)(x1)(x4)(x16)(x25)(x36)(x49)(x)l3(x)= (90)(91)(94)(916)(925)(936)(949)(9)(x-0)(x1)(x4)(x9)(x25)(x36)(x49)(x)l4(x)= (160)(161)(164)(169)(1625)(1636)(1649)(16)(x-0)(x-1)(x4)(x9)(x16)(x36)(x49)(x)l5(x)= (250)(251)(254)(259)(2516)(2536)(2549)(25)(x-0)(x1)(x4)(x9)(x16)(x25)(x49)(x)l6(x)= (360)(361)(364)(369)(3616)(3625)(3649)(36)（x-0）(x-1)(x4)(x9)(x16)(x25)(x36)(x)l7(x)= （490）(491)(494)(499)(4916)(4925（)4936)(49)(x-0)(x-1)(x4)(x9)(x16)(x25)(x36)(x49)l8(x)= (0)(1)(4)(9)(16)(25)(36)(49)l1(x)= L8(x)= l1(x)+2 l2(x)+3 l3(x)+4 l4(x)+5 l5(x)+6 l6(x)+7 l7(x)+8 l8(x) 求三次样条插值函数由MATLAB计算： 可得矩阵形式的线性方程组为： 00000000M002.00000.25002.00000.75001.0000000M100.37502.00000.625000000M20.1M00.41672.00000.583300000.028630 0000.43752.00000.5625000M40.0119M00000.45002.00000.5500000.006150M600000.45832.00000.1700.0035M0000000.46342.00000.53570.00227000000002.00000M8在MATLAB中的Editor中输入程序代码， 以下是三次样条函数的程序代码： function tgsanci(n,s,t) %n代表元素数，s,t代表端点的一阶导。 可编辑word,供参考版！

y=[0 1 2 3 4 5 6 7 8]; x=[0 1 4 9 16 25 36 49 ]; n=9 for j=1:1:n-1 h(j)=x(j+1)-x(j); end for j=2:1:n-1 r(j)=h(j)/(h(j)+h(j-1)); end for j=1:1:n-1 u(j)=1-r(j); end for j=1:1:n-1 f(j)=(y(j+1)-y(j))/h(j); end for j=2:1:n-1 d(j)=6*(f(j)-f(j-1))/(h(j-1)+h(j)); end d(1)=0 d(n)=0 a=zeros(n,n); for j=1:1:n a(j,j)=2; end r(1)=0; u(n)=0; for j=1:1:n-1 a(j+1,j)=u(j+1); a(j,j+1)=r(j); end b=inv(a) m=b*d' t=a p=zeros(n-1,4); %p矩阵为S(x)函数的系数矩阵 可编辑word,供参考版！

for j=1:1:n-1 p(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); end p 解得: M0=0;M1=-0.5209;M2=0.0558;M3=-0.0261;M4=0.0006;M5=-0.0029;M6=-0.0008;M7=--0.0009; M8=0，则三次样条函数： 3（01x）0.0868(x0)3（01x）1.0868(x-0)，x[0,1]33（4x）-0.0031（x-1）0.5938（4x）0.6388（x1），x[1,4]-0.02330.0019（9x）0.0009（x-4）0.3535（9x）0.6271（x4），x[4,9]33（16x）（0x9）0.4590（16x）-0.5708(x-9)，x[9,16]-0.0006S(x)=  033025x）-0.0001（x16）-0.4436（25x）0.5600(x-16)，x[16,25]（33（036x）-（0x25）0.4599（36x）0.70(x-25)，x[25,36]33（048x）-（0x36）0.4633（48x）0.04(x-36)，x[36,48]30x）0(x48)30.46（x）0.5333(x-48)，x[48,]（下面进行画图，在Command Window中输入画图的程序代码： %画图形比较那个插值更精确的函数： x0=[0 1 4 9 16 25 36 49 ]; y0=[0 1 2 3 4 5 6 7 8]; x=0:;y=lagr1(x0,y0,x); plot(x0,y0,'o') hold on plot(x,y,'r'); hold on; pp=csape(x0,y0,'variational') fnplt(pp,'g'); hold on; plot(x0,y0,':b');hold on %axis([0 2 0 1]); %看[0 1]区间的图形时加上这条指令 gtext('三次样条插值') gtext('原图像') gtext('拉格朗日插值') 可编辑word,供参考版！

%下面是求拉格朗日插值的函数 function y=lagr1(x0,y0,x) n=length(x0); m=length(x); for i=1:m z=x(i); s=0.0; for k=1:n p=1.0; for j=1:n if j~=k p=p*(z-x0(j))/(x0(k)-x0(j)); end end s=p*y0(k)+s; end y(i)=s; end 拉格朗日插值函数与三次样条插值函数如图中所示，绿色实线条为三次样条插值曲线，蓝色虚线条为x=y2的曲线，另外一条红色线条为拉格朗日插值曲线。图3-1为[0 1]的曲线，图3-2为[0 ]区间上的曲线。 图3—1 可编辑word,供参考版！

10.90.80.7拉格朗日插值0.60.50.40.30.20.1000.20.40.60.811.21.41.61.82三次样条插值原图像1008060拉格朗日插值4020三次样条插值0原图像-20010203040506070 图3—2 由图3-1可以看出，红色的线条与蓝色虚线条几乎重合，所以可知拉格朗日插值函数的曲线更接近开平方根的函数曲线，在[0，1]朗格朗日插值更精确。而在区间[0，]上从图3-2中可以看出蓝色虚线条和绿色线条是几乎重合的，而红色线条在[30，70]之间有很大的振荡，所在在区间[0，]三次样条插值更精确写。 【结论】（结果） 单个多项式高次插值效果并不理想，有龙格现象，偏差大，没有使用价值。而分段低次插可编辑word,供参考版！

值则精确度较高，拟合效果较好，而三次样条插值具有良好的收敛性与稳定性，与分段低次插值相比较光滑度更高，而且提供的信息也相对少一些。我们可以看到，在以上的三道实验题里，我们可以从图形中看出，三次样条的拟合程度是三种插值函数里最好的。 【小结】 通过此次实验，我对牛顿多项式插值，三次样条插值，拉格朗日插值有了更进一步的了解，知道了三次样条的拟合程度在高次的情况下更高，在理论上和应用上都有重要意义，在利用计算机编程软件进行高次插值的时候，我们可以多考虑利用三次样条进行插值。 指导教师评语及成绩： 成绩： 指导教师签名： 批阅日期：

【此文档部分内容来源于网络，如有侵权请告知删除，本文档可自行编辑和修改内容，感谢您的支持！】

可编辑word,供参考版！