matlab计算拉格朗日牛顿及分段线性插值的程序

《工程常用算法》综合实践作业二完成日期：2013年4月14日班级学号姓名主要工作说明自评成绩07182010071826崔洪亮算式与程序的编写1807182010071815侯闰上流程图的编辑，程序的审查07182010071809赵化川报告的整理汇总一.作业题目：三次样条插值与分段插值已知飞机下轮廓线数据如下：x035791112131415y01.21.72.02.12.01.81.21.01.6飞机下轮廓线形状大致如下图所示：要求分别用拉格朗日插值法、Newton插值法、分段线性插值法和三次样条插值法计算x每改变0.5时y的值，即x取0.5,1,1.5,…,14.5时对应的y值。比较采用不同方法的计算工作量、计算结果和优缺点。二.程序流程图及图形1.拉格朗日插值法开始x,y,x0Length(x)==lEngth(y)?n=length(x)i=1:n,l=1。j=1:i-1&j=i+1:nl=l.*(x0-x(j)/x(i)-x(j)f=f+l*y(i)结束否是机翼下轮廓线2.牛顿插值法开始x,y,xiLength(x)==length(y)?n=length(x)Y=zeros（n),Y（：1）=y,f=0a=1:n-1,b=1:n-a,Y(b,a+1)=（Y（b+1，a)-Y(b,a))/(x（b+a)-x(b))i=1:n,z=1结束j=1:i-1,z=z.*(xi-x(j))f=f+Y(1,i)*z否是3.分段线性插值法开始x，y，x0length（x）==length(y)?k=1:n-1x(k)=x0&x0《=x(k+1)temp=x(k)-x(k+1)f=(x0-x(k+1))/temp*y(k)+(x0-x(k))/(-temp)*y(k+1)结束否否是是三.matlab程序及简要的注释（m文件）1.拉格朗日插值法2.牛顿插值法functionf=newdun(x,y,xi)%x为已知数据点的x坐标向量%y为已知数据点的y坐标向量functionf=lang(x,y,x0)%x为已知数据点的x坐标向量%y为已知数据点的y坐标向量fora=1:n-13.分段线性插值法%xi为插值点的x坐标%f为求得的均差牛顿插值多项式if(length(x)==length(y))n=length(x);elsedisp('ThelengthsofXangYmustbenotequal!');return;end%检错Y=zeros(n);Y(:,1)=y;fora=1:n-1forb=1:n-aY(b,a+1)=(Y(b+1,a)-Y(b,a))/(x(b+a)-x(b));%计算均差函数（差商）endendf=0;fori=1:nz=1;forj=1:i-1z=z.*(xi-x(j));%计算多项式函数endf=f+Y(1,i)*z;%计算牛顿插值函数end%x0为插值点的x坐标%f为求得的拉格朗日插值多项式if(length(x)==length(y))n=length(x);elsedisp('ThelengthsofXangYmustbenotequal!')return;end%检错f=0;fori=1:nl=1;forj=1:i-1l=l.*(x0-x(j))/(x(i)-x(j));end;forj=i+1:nl=l.*(x0-x(j))/(x(i)-x(j));%计算拉格朗日基函数end;f=f+l*y(i);%计算拉格朗日插值函数endReturnfunctionf=fd2(x,y,x0)%x为已知数据点的x坐标向量%y为已知数据点的y坐标向量%x0为插值点的x坐标%f为求得的分段线性插值多项式if(length(x)==length(y))n=length(x);elsedisp('ThelengthsofXandYmustbenotequal!')return;end%检错fork=1:n-1if(x(k)=x0&x0=x(k+1))temp=x(k)-x(k+1);f=(x0-x(k+1))/temp*y(k)+(x0-x(k))/(-temp)*y(k+1);end;end四.程序运行结果1.拉格朗日插值法2.牛顿插值法3.分段线性插值法五、对不同实现方法的运行结果进行比较及总结x=[035791112131415];y=[01.21.72.02.12.01.81.21.01.6];xi=0.5:0.5:14.5;yi=lang(x,y,xi)yi=Columns1through8-15.4117-15.9238-10.9898-5.4272-1.22531.20002.17652.2666Columns9through161.98941.70001.57031.62491.79952.00002.14772.2040Columns17through242.17522.10002.02691.99041.99282.00001.95371.8000Columns25through291.52721.20000.96561.00001.3480plot(x,y,'b:',xi,yi)yi=newdun(x,y,xi)yi=Columns1through8-15.4117-15.9238-10.9898-5.4272-1.22531.20002.17652.2666Columns9through161.98941.70001.57031.62491.79952.00002.14772.2040Columns17through242.17522.10002.02691.99041.99282.00001.95371.8000Columns25through291.52721.20000.96561.00001.3480plot(x,y,xi,yi,'g+')yi=fd2(x,y,xi)yi=Columns1through80.20000.40000.60000.80001.00001.20001.32501.4500Columns9through161.57501.70001.77501.85001.92502.00002.02502.0500Columns17through242.07502.10002.07502.05002.02502.00001.90001.8000Columns25through291.50001.20001.10001.00001.3000plot(x,y,'b:',xi,yi,'g+')拉格朗日插值法的优点是表达式简单明确，形式对称，它的缺点是如果要想增加插值节点，整个公式必须都发生改变，且容易发生龙格现象。牛顿插值法却很好的改善了这一点，从而变得更加灵活方便。此外两者的拟合程度也很相似。分段插值的缺点是不能保证曲线在连接点处的光滑性。分段插值可以步进地插值计算，同时也带来了内在的高度稳定性和较好的收敛性。与前两种相比还具有良好的拟合性。总结：我们组只编写了前三种插值方法，对于三次样条插值参数太多，故没能编写出其程序，但通过matlab中的interp1函数可以实现。在此次作业我们也学到了很多知识，感到很充实。补充三次样条插值图形如下：4.三次样条插值法六.计算公式及计算方法（手写）yi=interp1(x,y,xi,'spline')yi=Columns1through80.24210.46650.67390.86491.04011.20001.34541.4767Columns9through161.59471.70001.79301.87401.94302.00002.04502.0773Columns17through242.09592.10002.08932.06682.03582.00001.94691.8000Columns25through291.50991.20001.01341.00001.1866plot(x,y,xi,yi)