佛山科学技术学院实验报告课程名称数值分析实验项目插值法与数据拟合专业班级机械工程姓名余红杰学号2111505010指导教师陈剑成绩日期月日一、实验目的1、学会Lagrange插值、牛顿插值和三次样条插值等基本插值方法;2、讨论插值的Runge现象3、学会Matlab提供的插值函数的使用方法,会用这些函数解决实际问题。二、实验原理1、拉格朗日插值多项式2、牛顿插值多项式3、三次样条插值三、实验步骤1、用MATLAB编写独立的拉格朗日插值多项式函数2、用MATLAB编写独立的牛顿插值多项式函数3、用MATLAB编写独立的三次样条函数(边界条件为第一、二种情形)4、已知函数在下列各点的值为:ix0.20.40.60.81.0()ifx0.980.920.810.640.38根据步骤1,2,3编好的程序,试分别用4次拉格朗日多项式4()Lx、牛顿插值多项式4()Px以及三次样条函数()Sx(自然边界条件)对数据进行插值,并用图给出{(,),0.20.08,0,1,2,,10iiixyxii},4()Lx、4()Px和()Sx。5、在区间[-1,1]上分别取10,20n用两组等距节点对龙格函数21(),(11)125fxxx作多项式插值,对不同n值,分别画出插值函数及()fx的图形。6、下列数据点的插值x01491625364964y012345678可以得到平方根函数的近似,在区间[0,64]上作图。(1)用这9个点作8次多项式插值8()Lx。(2)用三次样条(第一边界条件)程序求()Sx。7、对于给函数21()125fxx在区间[-1,1]上取10.2(0,1,,10)ixii,试求3次曲线拟合,试画出拟合曲线并打印出方程,与第5题的结果比较。四、实验过程与结果:1、Lagrange插值多项式源代码:functionya=lag(x,y,xa)%x所有已知插值点%y插值点对应函数值%xa所求点,自变量%ya所求点插值估计量ya=0;mu=1;%初始化%循环方式求L系数,并求和:fori=1:length(y)forj=1:length(x)ifi~=jmu=mu*(xa-x(j))/(x(i)-x(j));elsecontinueendendya=ya+y(i)*mu;mu=1;end2、Newton源代码:functionya=newton(x,y,xa)%x所有已知插值点%y插值点对应函数值%xa所求点,自变量%ya所求点插值估计量%建立系数零矩阵D及初始化:D=zeros(length(x)-1);ya=y(1);xi=1;%求出矩阵D,该矩阵第一行为牛顿插值多项式系数:fori=1:(length(x)-1)D(i,1)=(y(i+1)-y(i))/(x(i+1)-x(i));endforj=2:(length(x)-1)fori=1:(length(x)-j)D(i,j)=(D(i+1,j-1)-D(i,j-1))/(x(i+j)-x(i));endend%xi为单个多项式(x-x(1))(x-x(2))...的值fori=1:(length(x)-1)forj=1:ixi=xi*(xa-x(j));endya=ya+D(1,i)*xi;xi=1;end3、三次样条插值多项式(1)(第一边界条件)源代码:functiony=yt1(x0,y0,f_0,f_n,x)_____________(1)%第一类边界条件下三次样条插值;%xi所求点;%yi所求点函数值;%x已知插值点;%y已知插值点函数值;%f_0左端点一次导数值;%f_n右端点一次导数值;n=length(x0);z=length(y0);h=zeros(n-1,1);k=zeros(n-2,1);l=zeros(n-2,1);S=2*eye(n);fori=1:n-1h(i)=x0(i+1)-x0(i);endfori=1:n-2k(i)=h(i+1)/(h(i+1)+h(i));l(i)=1-k(i);end%对于第一种边界条件:k=[1;k];_______________________(2)l=[l;1];_______________________(3)%构建系数矩阵S:fori=1:n-1S(i,i+1)=k(i);S(i+1,i)=l(i);end%建立均差表:F=zeros(n-1,2);fori=1:n-1F(i,1)=(y0(i+1)-y0(i))/(x0(i+1)-x0(i));endD=zeros(n-2,1);fori=1:n-2F(i,2)=(F(i+1,1)-F(i,1))/(x0(i+2)-x0(i));D(i,1)=6*F(i,2);end%构建函数D:d0=6*(F(1,2)-f_0)/h(1);___________(4)dn=6*(f_n-F(n-1,2))/h(n-1);___________(5)D=[d0;D;dn];______________(6)m=S\D;%寻找x所在位置,并求出对应插值:fori=1:length(x)forj=1:n-1if(x(i)=x0(j+1))&(x(i)=x0(j))y(i)=(m(j)*(x0(j+1)-x(i))^3)/(6*h(j))+...(m(j+1)*(x(i)-x0(j))^3)/(6*h(j))+...(y0(j)-(m(j)*h(j)^2)/6)*(x0(j+1)-x(i))/h(j)+...(y0(j+1)-(m(j+1)*h(j)^2)/6)*(x(i)-x0(j))/h(j);break;elsecontinue;endendend(2)(自然边界条件)源代码:仅仅需要对上面部分标注的位置做如下修改:__(1):functiony=yt2(x0,y0,x)__(2):k=[0;k]__(3):l=[l;0]__(4)+(5):删除—(6):D=[0:D:0]4、——————————————PS:另建了一个f方程文件,后面有一题也有用到。functiony=f(x0)y=1./(1+25.*x0.^2);___________________________clc;clear;x1=[0.2,0.4,0.6,0.8,1.0];y1=[0.98,0.92,0.81,0.64,0.38];plot(x1,y1,'.');holdonxo=[0.2:0.08:1];y=lag(x1,y1,xo);plot(xo,y,'o')holdon;y=newton(x1,y1,xo);plot(xo,y,'r');holdon;y=yt2(x1,y1,xo);plot(xo,y,'*')h=legend('原始','拉格','牛顿','自样',4);5、clc,clear;x1=linspace(-1,1,10);x2=linspace(-1,1,20);xo=[-1:0.02:1];yo=f(xo);plot(xo,yo);holdon;y=f(x1);y=newton(x1,y,xo);plot(xo,y,'k');holdon;y=f(x2);y=newton(x2,y,xo);plot(xo,y,'r')h=legend('原始','10插','20插',3);6、clc,clear;x1=[01491625364964];y1=[012345678];xo=[0:1:64];y=lag(x1,y1,xo)plot(xo,y,'k');holdony=yt2(x1,y1,xo)plot(xo,y,'r')h=legend('拉格','自然样条',2);7、clc,clear;x1=linspace(-1,1,11);xo=[-1:0.02:1];y=f(x1);yo=f(xo);plot(xo,yo,'r');holdonp=polyfit(x1,y,3);y1=polyval(p,xo);plot(xo,y1)h=legend('原图','三次曲线拟合',2);p%该曲线的三次多项式系数依次显示五、讨论分析及感想个人感觉就是在数据点比较少,要求比较低的情况下,使用拉格朗日或是牛顿插值就足够了。但是当数据点比较多的时候,使用样条曲线就更好。当数据更多时,就可以使用曲线拟合方法来求近似值。工程数学的基础知识并不是很苦难,但是很实用,是有必要好好学下的。Matlab也是比较强大的工具,比较符合人的思维逻辑,上手很快。看书上的程序例子和实际编写还是有区别的,看得懂不一定编写的好,主要还是思维方式的锻炼吧~格式方法之类,只要了解了其基本功能,然后就是一系列的组合,多训练就熟悉了,比较有趣的课程吧,关键是有数据,有输出图像,看的清晰明白。而且出现的错误,都有清晰的指导,修改起来也很快捷。至于具体的分析,在编写三次样条的时候发现,使用矩阵的思想分析世界是有很大优势的,以后自己会多加训练。