插值法与数据拟合

佛山科学技术学院实验报告课程名称数值分析实验项目插值法专业班级姓名学号指导教师成绩日期5月12日一、实验目的1、学会Lagrange插值、牛顿插值和分段线性插值等基本插值方法；2、讨论插值的Runge现象，掌握分段线性插值方法3、学会Matlab提供的插值函数的使用方法，会用这些函数解决实际问题。二、实验原理1、拉格朗日插值多项式2、牛顿插值多项式3、分段线性插值三、实验步骤1、用MATLAB编写独立的拉格朗日插值多项式函数；2、用MATLAB编写独立的牛顿插值多项式函数；3、利用编写好的函数计算本章P66例1、P77例1的结果并比较；4、已知函数在下列各点的值为：ix0.20.40.60.81.0()ifx0.980.920.810.640.38试用4次牛顿插值多项式4()Px对数据进行插值，根据｛(,),0.20.08,0,1,2,,10iiixyxii｝，画出图形。5、在区间[-1,1]上分别取10,20n用两组等距节点对龙格函数21(),(11)125fxxx作多项式插值，对不同n值，分别画出插值函数及()fx的图形。6、下列数据点的插值x01491625364964y012345678(1)可以得到平方根函数的近似，在区间[0,64]上作图。（2）用这9个点作8次多项式插值8()Lx。四、实验结果1、用MATLAB编写独立的拉格朗日插值多项式函数Lagrange插值多项式源代码functionyi=Lagrange(x,y,xi)%Lagrange插值多项式，其中%x---向量，全部的插值节点%y---向量，插值节点处的函数值%xi---标量，自变量x%yi---xi处的函数估计值n=length(x);m=length(y);ifn~=merror('ThelengthsofXandYmustbeequal');return;endp=zeros(1,n);%对向量p赋初值0fork=1:nt=ones(1,n);forj=1:nifj~=kifabs(x(k)-x(j))epserror('theDATAiserror!');return;endt(j)=(xi-x(j))/(x(k)-x(j));endendp(k)=prod(t);endyi=sum(y.*p);2、用MATLAB编写独立的牛顿插值多项式函数functionyi=New_Int(x,y,xi)%Newton基本插值公式，其中%x---向量，全部的插值节点，按行输入%y---向量，插值节点处的函数值，按行输入%xi---标量，自变量x%yi---xi处的函数估计值n=length(x);m=length(y);ifn~=merror('ThelengthsofXandYmustbeequal');return;end%计算均差表YY=zeros(n);Y(:,1)=y';%Y(:,1)表示矩阵中第一列的元素fork=1:n-1fori=1:n-kifabs(x(i+k)-x(i))epserror('theDATAiserror!');return;endY(i,k+1)=(Y(i+1,k)-Y(i,k))/(x(i+k)-x(i));endend%计算Newton插值公式N(xi)yi=0;fori=1:nz=1;fork=1:i-1z=z*(xi-x(k));endyi=yi+Y(1,i)*z;end3、利用编写好的函数计算本章P66例1、P77例1的结果并比较。P66例1x=[144169225];y=[121315];yi=Lagrange(x,y,175)yi=13.230158730158733P77例1x=[0.40,0.55,0.65,0.80];y=[0.41075,0.57815,0.69675,0.88811];yi=New_int(x,y,0.596)yi=0.6319144055040004、已知函数在下列各点的值为：ix0.20.40.60.81.0()ifx0.980.920.810.640.38试用4次牛顿插值多项式4()Px对数据进行插值，根据｛(,),0.20.08,0,1,2,,10iiixyxii｝，画出图形。解：X=[0.2:0.2:1.0];y=[0.98,0.92,0.81,0.64,0.38];xx=[0.2:0.08:1.0];m=length(xx);z=zeros(1,m);fori=1:mz(i)=Lagrange(x,y,xx(i));endholdonplot(x,y,'o');plot(xx,z,'r*');holdoff得到如下图形：图一练习4的图形5、在区间[-1,1]上分别取10,20n用两组等距节点对龙格函数21(),(11)125fxxx作多项式插值，对不同n值，分别画出插值函数及()fx的图形。解：a=-1;b=1;n=100;h=(b-a)/n;x=a:h:b;y=1./(1+25.*x.^2);plot(x,y,'k')其函数原图形分别如下所示：-1-0.500.5100.10.20.30.40.50.60.70.80.91xy图二龙格函数的图形用龙格函数的Lagrange（）插值函数画图源程序当n=10时，有：functionRunge(10)%Runge现象%n---等距离节点a=-1;b=1;h=(b-a)/n;x=[a:h:b];y=1./(1+25.*x.^2);xx=[a:0.01:b];yy=1./(1+25.*xx.^2);m=length(xx);z=zeros(1,m);fori=1:mz(i)=Lagrange(x,y,xx(i));endholdonplot(x,y,'o');plot(xx,z,'r-');holdoff当n=20时，有：functionRunge(10)%Runge现象%n---等距离节点a=-1;b=1;h=(b-a)/n;x=[a:h:b];y=1./(1+25.*x.^2);xx=[a:0.01:b];yy=1./(1+25.*xx.^2);m=length(xx);z=zeros(1,m);fori=1:mz(i)=Lagrange(x,y,xx(i));endholdonplot(x,y,'o');plot(xx,z,'r-');holdoff其图形分别如下所示：图三Runge(10)图四Runge(20)6、下列数据点的插值x01491625364964y012345678(1)可以得到平方根函数的近似，在区间[0,64]上作图。（2）用这9个点作8次多项式插值8()Lx。functionL8%平方根函数xy的8次多项式插值)(8xLx=[01491625364964];y=sqrt(x);m=length(x);z=zeros(1,m);fori=1:mz(i)=Lagrange(x,y,x(i));endholdonplot(x,y,'o');plot(x,y,'o');xlabel('x');ylabel('y');plot(x,z,'r-');holdoff得到其图形如下所示：图五xy的图形x=[01491625364964];y=[012345678];xx=[0:0.5:64];yy=sqrt(xx);m=length(xx);z=zeros(1,m);fori=1:mz(i)=Lagrange(x,y,xx(i));endplot(x,y,'o',xx,yy,'.',xx,z,'r-');xy的8次Lagrange插值函数图形图六xy的8次Lagrange插值函数图形五、讨论分析从下面这段代码可以看出，利用mathlab自带的增长型变量来作图，是一种在c语言学习中没有发现的方法。而且，代码中还体现了一种比较作图和数据共用的思维方法。例如：x=[01491625364964];y=[012345678];xx=[0:0.5:64];yy=sqrt(xx);m=length(xx);z=zeros(1,m);fori=1:mz(i)=Lagrange(x,y,xx(i));endplot(x,y,'o',xx,yy,'.',xx,z,'r-');利用以上简单的程序，matble就可以做出图形，方便观察和研究。我们知道，利用计算机则可以方便地增加曲线上样点的个数，随着样点的个数增加，让学生观察函数的图像性质，从而抽象出函数图像的形状。当然，作图的话我们也可以利用z＋z超级画板。下面我们完成这个过程：（1）在坐标系的属性对话框中，选择“画坐标网格”和“显示刻度”选项。（2）在作图区空白处单击鼠标右键，在快捷菜单中单击【函数或参数方程曲线…】命令，如下图所示弹出函数作图属性对话框，在“xy”对应的编辑框中输入：√x，然后在下面的曲线属性中，设置“曲线的点数”为：9，选择曲线以“折线段”方式显示，设置变量x的参数范围为：0到64，选择在折线段上“画点”，设置点的大小为：2；最后单击【确定】按钮完成,同样可以画出函数图像，通过改变曲线的点数，函数曲线越来光滑；六、改进实验建议三次样条插值直接用spline函数做。边界条件加在y的首尾，第一个表示y'(x0)，最后一个表示y'(xt)。如果不加边界条件，默认是not-a-knot边界条件（注意不是自然边界条件）自然边界条件的插值要用csape函数才能得到。如果用interp1，则只能使用spline函数的默认边界条件，即not-a-knot条件。下面是例子:x=0:3:9;y=x.*cos(x);xx=linspace(0,9);plot(x,y,'o');%样本点holdon;plot(xx,interp1(x,y,xx,'spline'),'r');%interp1只能使用默认边界条件plot(xx,spline(x,[0y0],xx),'r:');%spline可以使用第一类边界条件，这里y'(0)=y'(9)=0pp=csape(x,y,'second');plot(xx,fnval(pp,xx))%第二类边界条件要用csape做，这里自然边界条件legend('样本点','默认边界条件','一阶导数为0','自然边界条件','location','south')holdoff;另外，我们可以设：y=a0+a1*x+a2*x^2+…+a(n-1)*x^(n-1)，则其n次拉格朗日函数插值程序：functionLagrangesNs()%用于求过n点的拉格朗日n-1次插值多项式options={'Nameofdatafile'};title='Lagranges_points';lineNo=2;def={'Lagranges.dat'};outval=inputdlg(options,title,lineNo,def);ifisempty(outval)==1,return,endfilename=outval{1};data=load(filename);x=data(:,1);y=data(:,2);lagrangesN(x,y);endfunctionlagrangesN(x,y)%画出已知n个点的位置plot(x,y,'*');holdon%n次拉格朗日多项式为y=a0+a1*x+a2*x^2+…+a(n-1)*x^(n-1)%其中a0a1a2…a(n-1)为待求系数n=length(x);X=Vandermonde(x,1);A=X\y;%绘制插值函数图象x1=linspace(0,max(x));x2=Vandermonde(x1',n);y1=x2*A;plot(x1',y1);holdon%显示公式func=['y=',num2str(A(1))];fori=2:n;b=['+',num2str(A(i)),'*x^',num2str(i-1)];func=[func,b];endtext(0.8,0.8,func);end%创建一个Vandermonde行列