多项式插值[学术参考]

1、多项式插值典型算法及其应用一、任务简介21.1任务题目21.2任务目的21.3任务案例2二、总体设计3三、 算法原理43.1分段线性插值43.2保形插值43.3三次样条插值53.4最邻近点插值93.5分片线性插值93.6双线性插值9四、 算法介绍114.1分段线性插值114.2保形插值124.3三次样条插值134.4最邻近点插值144.5分片线性插值154.6双线性插值16五、软件功能16六、运行结果176.1主界面176.2二级界面176.3程序结果18七、个人总结22八、参考资料23附录23一、任务简介1.1任务题目多项式插值典型算法及其应用1.2任务目的基于MATLAB图形界面对测试案例

2、利用分段线性插值、保型插值、三次样条插值、最邻近点插值、分片线性插值、双线性插值实现。1.3任务案例案例一、数控机床加工零件待加工零件的外形根据工艺要求由一组数据（x,y）给出（在平面情况下），用数控机床加工时刀具必须沿这些数据点前进，并且由于刀具每次只能沿x方向或y方向走非常小的一步，所以需要将已知数据加密，得到加工所要求的步长很小的（x,y）坐标。图1是待加工零件的轮廓线，表1给出了轮廓线上x每间隔0.2（长度单位）的加工坐标x,y（顺时针方向为序，由轮廓线的左右对称性，表中只给出右半部的数据），假设需要得到x或y坐标每改变0.05时的坐标，试完成加工所需的加密数据，画出曲线。图1、零件的

3、轮廓线（x间隔0.2）表1、x间隔0.2的加工坐标（x,y）（图1右半部的数据）00.20.40.60.811.254.714.313.683.052.52.051.41.61.822.22.42.61.691.41.1810.860.740.642.833.23.43.63.840.570.50.440.40.360.320.294.24.44.64.854.84.60.260.240.20.150-1.4-1.964.44.243.83.63.43.2-2.37-2.71-3-3.25-3.47-3.67-3.8432.82.62.42.221.8-4-4.14-4.27-4.39-4.4

4、9-4.58-4.661.61.41.210.80.60.4-4.74-4.8-4.85-4.9-4.94-4.96-4.980.20-4.99-5案例二：山区地貌在某山区测得一些地点的高度如下表，平面区域，试作出该山区的地貌图。36003200280024002000160012001480 1500 1550 1510 1430 1300 1200 9801500 1550 1600 1550 1600 1600 1600 15501500 1200 1100 1550 1600 1550 1380 10701500 1200 1100 1350 1450 1200 1150 101013

5、90 1500 1500 1400 900 1100 1060 9501320 1450 1420 1400 1300 700 900 8501130 1250 1280 1230 1040 900 500 700Y/x1200 1600 2000 2400 2800 3200 3600 4000二、总体设计多项式插值案例二案例一双线性插值分片线性插值最邻近点插值三次样条插值保形插值分段线性插值三、 算法原理3.1分段线性插值给定n+1 个结点a=x0x1xn =b上的函数值y0, y1,yn ，求一折线函数Ih(x)满足：则称Ih(x)为分段线性插值函数。由定义可知Ih(x)在每一个小区间x

6、k,xk+1上可表示为若用插值基函数表示，则在整个区间a,b上Ih(x)为3.2保形插值给定a, b 上的一种分划及f(x) 在分划点上的函数值, 及 ，(k=0,1,n),则作一个分段三次Hermite 插值多项式 P(x) , 满足条件 P(x)C1a, b； = , , k=0,1,n； P(x) 在每个小区间上是三次多项式。分段三次Hermit插值函数的分段表达式为：3.3三次样条插值并且关于这个节点集的三次样条函数s(x)满足插值条件：则称这个三次样条函数s(x)为三次样条插值函数。三次样条插值函数的边界条件:(1) 插值条件：(2) 连续性条件：(3) 一阶导数连续条件：(4) 二

7、阶导数连续条件：因为s(x)在每个小区间上是一个次小于三次的多项式，故有四个未知系数；因为s(x)有n分段，从而共有4n个未知系数！但插值条件与样条条件仅给出4n-2个条件，无法定出4n个未知系数，还差2个条件！这2个条件我们用边界条件给出！第一边界条件：由区间端点处的一阶导数给出即第二边界条件：由区间端点处的二阶导数给出即特殊情况为自然边界条件：由区间端点处的二阶导数恒为0给出即第三类又称周期边界条件：由区间端点处的函数值或导数值满足周期条件给出为了确定三次样条插值函数的表达式 S(x)，我们采用待定系数法来求解，考虑到带一阶导数的分段三次Hermite插值多项式我们采用待定一阶导数的方法即

8、设因为分段三次Hermite插值多项式已经至少是一阶连续可导了，为了让它成为三次样条函数只需确定节点处的一阶导数使这些节点处的二阶导数连续即可。 由于在内部节点处二阶导数连续条件： 整理化简后得： 第一类三次样条插值问题方程组由于已知： 基本方程组化为n-1阶方程组化为矩阵形式这是一个严格对角占优的三对角方程组，用追赶法可以求解。第二类三次样条插值问题的方程组由于已知：联合基本方程组得一个n+1阶三对角方程组，化成矩阵形式为：仍然是严格对角占优第三类样条插值问题的方程组联合基本方程得一个广义三对角或周期三对角方程组，这个方程组的系数矩阵仍然是严格对角占优阵。3.4最邻近点插值 最邻近点插值取插

9、值点的4个邻点中距离最近的邻点值作为该点的值。最邻近插值一般不连续。具有连续性的最简单的插值是分片线性插值。3.5分片线性插值四个插值点（矩形的四个顶点）处的函数值：第一片(下三角形区域)： 插值函数为：第二片(上三角形区域)：插值函数为：由于(x, y)当然应该是在插值节点所形成的矩形区域内。显然，分片线性插值函数是连续的。3.6双线性插值双线性插值是一片一片的空间二次曲面构成。插值函数的形式如下：假如我们想得到未知函数f在点P=(x,y)的值，假设我们已知函数f在Q11=x1,y1,Q12=x1,y2,Q21=(x2,y1)及Q22=(x2,y2)四个点的值。首先在x方向进行线性插值，得到

10、fR1x2-xx2-x1fQ11+x-x1x2-x1fQ21 R1=x,y1,fR2x2-xx2-x1fQ12+x-x1x2-x1fQ22 R2=(x,y2)然后在y方向进行线性插值，得到fPy2-yy2-y1fR1+y-y1y2-y1f(R2)可得到所要的结果f(x,y)fx,yfQ11x2-x1y2-y1x2-xy2-y+fQ21x2-xy2-y1x-x1y2-y+fQ12x2-x1y2-y1x2-xy-y1+fQ22x2-x1y2-y1x-x1(y-y1)四、 算法介绍4.1分段线性插值输入已知点x,y判断x与y维数相等将x的区间，划分为若干个小区间在每个小区间上构造插值多项式判别插值点

11、所属区间计算所有插值点的函数值画出图形输出开始结束是否是开始输入已知点x,y及插值点x划分区间在每个小区间上构造多项式fx判断插值点所处区间计算所有插值x的函数值在此处键入公式。画出图形结束端点导数值用均差代替x,y,y的导数值个数是否相同否4.2保形插值 4.3三次样条插值开始根据区间a,b，利用边界条件整理结合得到的式子求线性方程组输入数据用追赶法求解判定插值点x在区间输出所求函数值图形结束4.4最邻近点插值开始输入已知点坐标x,y,z,x0,y0根据x与y的等分，划分若干方形区域判断插值点所在区域将该方形区域再划分成四个小方形区域判断插值点所在的小方形区域将插值点代入对应的小方形区域内的

12、插值函数得出相应z值输出所有z值画出图像结束4.5分片线性插值开始输入已知数据根据x的等分与y的等分，划分若干行区域判断插值点所在区域判断是否为方形区域的下三角将x0,y0代入下三角插值函数，得出相应z0将x0,y0代入上三角插值函数，得出相应z0输出所有z0，画出图形结束4.6双线性插值开始输入已知数据根据x与y的等分，划分若干方形区域判断x0,y0所在区域将x0,y0代入对应区域插值函数求z0结束输出所有z0画出图形五、软件功能可以针对已知点的二维坐标、三维坐标分别进行一维插值、二维插值，使计算机模拟出实际的图形。六、运行结果6.1主界面 6.2二级界面6.3程序结果分段线性插值：分段线性

13、插值所得结果与原轮廓线一致，因为原轮廓线也是由两点之间连线画出的。分段插值优点：只要小区间长度足够小，便可保证插值函数充分靠近原函数，即整体逼近效果好且计算简单快捷。分段插值缺点：图像光滑性差保形插值：保形插值优点：比分段线性插值效果明显改善。保形插值缺点：但是这种插值要求给出节点上的导数值，所要提供的信息太多，光滑度也不高。三次样条插值：三次样条插值优点：输出结果最为平滑，不但与被插值函数很接近，而且导数值也很接近，这样逼近效果是其他插值法所难以达到的。三次样条插值缺点：花费时间最多最邻近点插值：最邻近点插值优点：执行速度快；缺点：不连续，输出结果为直角转折分片线性插值：分片线性插值优点：简

14、单且连续双线性插值：双线性插值优点：效果最好；缺点：运行耗时长七、个人总结回顾起此次课程设计，至今我仍感慨颇多，的确，自从拿到题目到完成整个编程，从理论到实践，在不到一星期的日子里，可以学到很多很多的东西，同时不仅可以巩固了以前所学过的知识，而且学到了很多在书本上所没有学到过的知识。虽然因此开启了整日往返于机房、食堂、宿舍的“三点一线”的生活，虽然看上去着实不免枯燥乏味，但现在接近尾声之际又回想起那一次次运行成功的激动时分也的确让我们怀念不已。就要结束数值分析课程设计，也意味着我们整日忙于双手敲打键盘，双眼凝视屏幕编程的日子就要告一段落了，还真有些舍不得。在这一周的程序设计学习中，我认识到了许

15、多知识方面的不足，在编写程序时，我不能快速的编出程序，很多代码我都很生疏，很多知识上的漏洞暴露无遗，通过一周的学习，我慢慢弥补我的漏洞，克服我的缺点，在我认真学习下，我在慢慢的进步，慢慢的完善我自己的知识。使我对MATLAB有了更进一步的认识和了解，要想学好它要重在实践，要通过不断的上机操作才能更好地学习它，通过实践，我也发现我的好多不足之处，有对MATLAB的一些标准库函数不太了解，还有对函数调用的正确使用不够熟悉，还有对MATLAB中经常出现的错误也不了解，通过实践，使我在这几个方面的认识有所提高。虽然我们组只有两个人，但我们齐心协力，刻苦努力，圆满完成本次课设。在以后的学习生活中，我会以

16、这次课设为经验和教训，更加努力更加认真，争取做得更好。这次课程设计使我懂得了理论与实际相结合是很重要的，只有理论知识是远远不够的，只有把所学的理论知识与实践相结合起来，从理论中得出结论，才能提高自己的实际动手能力和独立思考的能力。在设计的过程中可以说得是困难重重，同时在设计的过程中发现了自己的不足之处，对一些前面学过的知识理解得不够深刻，掌握得不够牢固，这次课程设计相当于我们把前面所学过的知识又重新温故了一遍。将以前的知识串了起来，使自己的编程思路更为清晰流畅，运用起来更是有一气呵成的自然而然之感。一次又一次的出错、改错、调试、运行都是对知识的巩固练习，尤其是函数知识的运用，在编写这种大程序时

17、更是方便自如。所以我通过这次课程设计又掌握了新的知识，丰富了自己的储备量。我体会到了团队合作的重要性，提升了自己的协作能力；我养成了及时注释、一小段一小段查错、准确编写正规代码的习惯，形成了建立文件的思想；通过这一个星期的磨合，自己变得对待事物更有耐心了，不管喜不喜欢，也能坚持到底，认真对待完成了。只有通过自己亲自动手做，才能真正体会到汗水的意义所在，汗水预示着结果也见证着收获。“实践出真知”，它同样也是检验真理的唯一标准，所以，不要再一味的局限于书本知识，真正的去实践才会拥有更深的体会！这次课程设计让我受益匪浅！八、参考资料【】张德丰，MATLAB数值分析与应用，北京：国防工业出版社，20

18、年 月。【】王能超，数值分析简明教程，北京：高等教育出版社，20 年月。附录分段线性插值：clc,clear;%x1=0:0.2:5;y1= 5 4.71 4.31 3.68 3.05 2.5. 2.05 1.69 1.4 1.18 1 0.86. 0.74 0.64 0.57 0.5 0.44 0.4. 0.36 0.32 0.29 0.26 0.24. 0.2 0.15 0 y2=-5 -4.99 -4.98 -4.96 -4.94 -4.9. -4.85 -4.8 -4.74 -4.66 -4.58 -4.49. -4.39 -4.27 -4.14 -4 -3.84 -3.67 . -3

19、.47 -3.25 -3 -2.71 -2.37 -1.96 . -1.4 0; u1=0:0.05:5;n=length(u1);for i=1:nv1(i)=fenduan2(x1,y1,u1(i)endfor i=1:nv2(i)=fenduan2(x1,y2,u1(i)endplot(x1,y1,r,x1,y2,r),axis(-5 5 -5 5)hold onplot(-x1,y1,r,-x1,y2,r),axis(-5 5 -5 5)gtext(,FontSize,8)plot(u1,v1,b,u1,v2,b),axis(-5 5 -5 5)hold onplot(-u1,v1,b

20、,-u1,v2,b),axis(-5 5 -5 5)gtext(,FontSize,8)function f=fenduan2(x,y,x0) for k=1:length(x)-1 if(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); endend保形插值：clc,clear;%x1=0:0.2:5;y1= 5 4.71 4.31 3.68 3.05 2.5. 2.05 1.69 1.4 1.18 1 0.86. 0.74 0.64 0.57 0.5 0.44 0.4.

21、0.36 0.32 0.29 0.26 0.24. 0.2 0.15 0 y2=-5 -4.99 -4.98 -4.96 -4.94 -4.9. -4.85 -4.8 -4.74 -4.66 -4.58 -4.49. -4.39 -4.27 -4.14 -4 -3.84 -3.67 . -3.47 -3.25 -3 -2.71 -2.37 -1.96 . -1.4 0; u1=0:0.05:5;n=length(u1);for i=1:nv1(i)=baoxing2(x1,y1,u1(i)endfor i=1:nv2(i)=baoxing2(x1,y2,u1(i)endplot(x1,y1,r

22、,x1,y2,r),axis(-5 5 -5 5)hold onplot(-x1,y1,r,-x1,y2,r),axis(-5 5 -5 5)gtext(,FontSize,8)plot(u1,v1,b,u1,v2,b),axis(-5 5 -5 5)hold onplot(-u1,v1,b,-u1,v2,b),axis(-5 5 -5 5)gtext(,FontSize,8)function yi=baoxing2(x,y,xi) ydot=gradient(y,x); n=length(x);m1=length(y);m2=length(ydot);% x,yy. if n=m1|n=m2

23、|m1=m2 error(The length of X,Y and Ydot must be equal!); return; end for k=1:n-1 if x(k)=xi&xi=x(k+1) yi=y(k)*(1+2*(xi-x(k)/(x(k+1)-x(k). *(xi-x(k+1)2/(x(k)-x(k+1)2+. y(k+1)*(1+2*(xi-x(k+1)/(x(k)-x(k+1). *(xi-x(k)2/(x(k+1)-x(k)2+. ydot(k)*(xi-x(k)*(xi-x(k+1)2/(x(k)-x(k+1)2+. ydot(k+1)*(xi-x(k+1)*(xi

24、-x(k)2/(x(k+1)-x(k)2; return; end end三次样条插值：clc,clear;%x1=0:0.2:5;y1= 5 4.71 4.31 3.68 3.05 2.5. 2.05 1.69 1.4 1.18 1 0.86. 0.74 0.64 0.57 0.5 0.44 0.4. 0.36 0.32 0.29 0.26 0.24. 0.2 0.15 0 ;y2=-5 -4.99 -4.98 -4.96 -4.94 -4.9. -4.85 -4.8 -4.74 -4.66 -4.58 -4.49. -4.39 -4.27 -4.14 -4 -3.84 -3.67 . -3

25、.47 -3.25 -3 -2.71 -2.37 -1.96 . -1.4 0; u1=0:0.05:5;for i=1:101v1(i)=sanci2(x1,y1,u1(i)endfor i=1:101v2(i)=sanci2(x1,y2,u1(i)endplot(x1,y1,r,x1,y2,r),axis(-5 5 -5 5)hold onplot(-x1,y1,r,-x1,y2,r),axis(-5 5 -5 5)gtext(,FontSize,8)plot(u1,v1,b,u1,v2,b),axis(-5 5 -5 5)hold onplot(-u1,v1,b,-u1,v2,b),ax

26、is(-5 5 -5 5)gtext(,FontSize,8)function f0=sanci2(x,y,x0)syms t; n=length(x); for i=1:n if(x(i)=x0) index=i; break; endend %x0A=diag(2*ones(1,n); %m 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)

27、; c(i)=3*lamda(i)*(y(i)-y(i-1)/(x(i)-x(i-1)+3*u(i-1)*(y(i+1)-y(i)/(x(i+1)-x(i); A(i,i+1)=u(i-1); A(i,i-1)=lamda(i); %cendc(1)=3*(y(2)-y(1)/(x(2)-x(1)-(x(2)-x(1)*0/2;c(n)=3*(y(n)-y(n-1)/(x(n)-x(n-1)-(x(n)-x(n-1)*0/2;m=followup(A,c); %h=x(index+1)-x(index); %f=y(index)*(2*(t-x(index)+h)*(t-x(index+1)2

28、/h/h/h+y(index+1)*(2*(x(index+1)-t)+h)*(t-x(index)2/h/h/h+m(index)*(t-x(index)*(x(index+1)-t)2/h/h-m(index+1)*(x(index+1)-t)*(t-x(index)2/h/h; %x0f0=subs(f,t,x0); %x0 function x=followup(A,b)n=rank(A);for(i=1:n) if(A(i,i)=0) disp(Error:0); return; endendd=ones(n,1);a=ones(n-1,1);c=ones(n-1);for(i=1:

29、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) x(i,1)=(b(i,1)-c(i,1)*x(i+1,1)/d(i,1);end最邻近点插值：x=1200 1600 2000 2400 2800 3200 3600 4000;y=

30、1200 1600 2000 2400 2800 3200 3600 ;z=1130 1250 1280 1230 1040 900 500 700; 1320 1450 1420 1400 1300 700 900 850; 1390 1500 1500 1400 900 1100 1060 950; 1500 1200 1100 1350 1450 1200 1150 1010; 1500 1200 1100 1550 1600 1550 1380 1070; 1500 1550 1600 1550 1600 1600 1600 1550; 1480 1500 1550 1510 1430

31、 1300 1200 980;x0=1200:30:4000;y0=1200:30:3600;for i=1:length(y0) for j=1:length(x0) z0(i,j)=zuilin2(x,y,z,x0(j),y0(i); endend x0,y0=meshgrid(x0,y0);meshc(x0,y0,z0)function z0=zuilin2(x,y,z,x0,y0) for i=1:length(y)-1 for j=1:length(x)-1 if(y(i)=y0&y0=(y(i)+100)&x(j)=x0&x0=(x(j)+100) z0=z(i,j); end i

32、f(y(i)=y0&y0=(y(i)+100)&(x(j)+100)=x0&x0=(x(j+1) z0=z(i,j+1); end if(y(i)+100)=y0&y0=y(i+1)&x(j)=x0&x0=(x(j)+100) z0=z(i+1,j); end if(y(i)+100)=y0&y0=y(i+1)&(x(j)+100)=x0&x0=x(j+1) z0=z(i+1,j+1); end endendend分片线性插值：x1=1200 1600 2000 2400 2800 3200 3600 4000;y1=1200 1600 2000 2400 2800 3200 3600 ;z1

33、=1130 1250 1280 1230 1040 900 500 700; 1320 1450 1420 1400 1300 700 900 850; 1390 1500 1500 1400 900 1100 1060 950; 1500 1200 1100 1350 1450 1200 1150 1010; 1500 1200 1100 1550 1600 1550 1380 1070; 1500 1550 1600 1550 1600 1600 1600 1550; 1480 1500 1550 1510 1430 1300 1200 980;x0=1200:30:4000;y0=1200:30:3600; for i=1:length(y0) for j=1:length(x0) z0(i,j)=fenpian2(x1,y1,z1,x0(j),y0(i); endend x0,y0=meshgrid(x0,y0);meshc(x0,y0,z0)function z0=fenpian2(x,y,z,x0,y0) z0=;m=length(y);n=length(x);for i=1:m-1 for j=1:n-1 if(y(i