二元函数插值及其程序设计--毕业论文.doc

目 录摘要2关键词21前言21.1二元函数插值及其发展过程21.2本文所要达到的目的32二元函数插值32.1一元 Lagrange 插值的构造方法32.2二元函数插值的基本思想42.3二元函数插值的几种方法52.3.1分片双一次插值52.3.2分片不完全双二次插值72.3.3矩形域上分片双三次埃尔米特插值92.4二元函数插值程序设计112.4.1MATLAB 中插值描述及程序设计113总述15致谢16参考文献16英文摘要17二元函数插值及其程序设计摘要 本篇文章主要对二元函数插值进行了叙述。针对一元函数插值思想主要是拉格朗日插值，我们将其中构造基函数的方法推广到二元函数，讨论了二元函数的插值问题。其中，主要讨论矩形区域上的插值、分片低次插值，将矩形域上分片插值问题分作分片双一次插值，分片不完全的双二次插值。并且针对插值做了MATLAB的程序设计，简单分析了插值问题的解决办法。关键词 二元函数插值；拉格朗日；MATLAB；分片双一次插值；分片不完全双二次插值；矩形区域；1前言1.1二元函数插值及其发展过程二元函数插值在生活中有着广泛的应用。例如在计算几何与计算辅助几何设计中有着重要的作用。在许多实际问题及科学研究中，因素之间往往存在着函数关系，然而这种关系很难有明显的解析表达式，通常只是由观察或测试得到一些离散数值。有时，即使给出了解析表达式，却由于表达式过于复杂，不仅使用不便，而且不易于进行计算和理论分析。从几何角度来说，就是要由给定的这组数据点去描绘曲线的近似图形。解决这类问题的方法有两种：插值法和曲线拟合法。二元逼近是一元函数逼近理论的发展，是在逼近工具和被逼近对象方面的二元推广。由于现代科学技术的发展的需要，二元函数逼近理论的研究日益受到数学、计算机数学、物理及工程等领域的专家和科学工作者的重视，已成为当今逼近理论和计算数学的研究热点之一。 在许多实际问题中需要建立模拟曲面，例如飞机、汽车、轮船和雕塑零件的外形设计以及一些描述科学现象的曲面的拟合等等。在这些外形设计和曲面拟合中经常应用局部多元插值方法。多元插值是一个活跃的研究领域，至今，已有相当多的多元插值公式，并且还在与日俱增。在这里，我们主要讨论二元插值公式。插值公式的产生和发展是由实际问题决定的，因此它们各有特色，但是应该说可供应用的公式还很少。基于矩阵网格上的矩形曲面片是容易理解的，因此得到使用者们欢迎。只要有可能，大家尽量用矩形曲面片。当然，在应用矩形曲面片时，需要注意相容性条件。一般来说，曲面是很复杂的，仅通过数值表示是难于想象的，因此需要作图或实际的三维模型。 1.2本文所要达到的目的 本文的目的是将一元函数插值思想如拉格朗日插值中构造基函数的方法推广到二元函数，讨论二元函数的插值问题。主要讨论矩形区域上的插值、分片低次插值等。 2二元函数插值2.1一元 Lagrange 插值的构造方法我们首先来看在一元函数中的Lagrange插值多项式，对给定了插值点来求如式的插值多项式的方法有多种。当节点是n+1时，可以先构造函数，它们的次数不超过n,且满足 当时，或者当j=i时，此为（2.1.1）然后以对应点处的函数值为系数作线性组合，即得所要求的插值多项式。下面推导的表达式。 由式（2.1.1），多项式 有 n 个根 ，且，故它必定是如下形式 （2.1.2）这些函数称为 Lagrange 插值基函数。利用它们立即得出插值问题的解 （2.1.3）事实上，因为每个插值 基函数 li (x)(i = 0,1,L, n) 都是 n 次多项 式。由式（2.1.1）又得 即 满足插值条件。式（2.1.3）称为 n 次 Lagrange 插值多项式。为了以后便于区别，常用 由前面讨论可得结论， n + 1 个节点的 n 次 Lagrange 插值多项式存在唯一。2.2二元函数插值的基本思想 二元函数分片光滑逼近是随着电子计算机的广泛应用而活跃起来的一个研究领域，在应用上也很重要。这类问题一般提法是，给定了被逼近曲面或函数 u = u (x, y ) ，或者是给定了u(x,y)的一组离散的近似值，要求构造一个比较简单的函数U(x,y)去逼近u(x,y)或离散值。若，则称为插值逼近或简称插值；通常由于 总有观测误差，因此并不要求，只要近似满足就行，近似通过给定点的曲面逼近法，我们称为曲面拟合法。我们主要介绍曲面插值法。 逼近二元函数最简单的函数类，自然是二元多项式。但实际问题往往给定的点很多，如果用一个整片多项式去逼近，则必然使得多项式次数很高，效果并不好。因此类似于一元函数的分段多项式逼近或样条函数逼近，这里就采用分片二元多项式逼近，并使不同曲面片之间光滑地联接起来。这种分片逼近的方法应十分值得注意。 例如，在用有限元方法解弹性力学问题或其它数学物理问题时，首先要对定解区域 进行几何剖分，也就是将 划分成一定的网格（或称单元）。当 是二维区域时，最基本的单元是三角形和四边形，也可以是曲边三角形和曲边四边形。区域剖分后紧接着的问题是在这些小单元上选取一近似函数U(x, y)去代替数学物理问题中的解u(x, y)，如果选取U(x, y)使它与u(x, y)在小单元的某些点的值相等，那么U(x, y)就称为u(x, y)在小单元上的插值函数。由于计算机辅助几何设计的发展，又提出了另一类的曲面逼近问题。美国的康斯提出了描述曲面的一种数学方法。一张曲面可以用若干个小的曲面片拼起来，适当选择曲面片的方程，使它联接起来保持一定的光滑性；换句话说，康斯曲面是由它的边界条件决定的，这又是另一类的插值问题。插值函数类的选择，最简单的方法是二元双k次式，即选择函数类G: （2.2.1）为插值函数。在研究插值函数的具体构造时，我们强调模仿一元函数的插值函数构造法。例如对于n个互异节点的二元插值函数U(x, y)，我们只要找到n个属于G的函数类，使它满足，当时（i,j=1,2,，n），或者，当i=j时（i,j=1,2,，n），此式为（2.2.2）。那么函数 （2.2.3）就是G中满足n个插值条件，i=1,2,，n (2.2.4)的插值函数。这是模仿一元函数拉格朗日的插值方法。对于某些插值问题，我们可先固定x（或y），将u(x, y)看成是y的函数，用一元插值的方法得到插值函数，简记之为。 然后将函数对x进行插值，得到的插值函数，记之为。对于某类特定的问题，就是满足插值条件的插值函数。这种方法称为乘积型插值法。 待定系数法是寻求一元函数插值法的重要方法，我们在构造二元插值函数时也要灵活运用这一数学方法。2.3二元函数插值的几种方法2.3.1分片双一次插值 在 Oxy 平面给定一个矩形域,其顶点坐标分别为（参见图1）。 Y O X 图 1 如果在角点上给定了函数值（或简记为）i = 1,2,3,4。现在我们来寻求它的插值函数U(x, y)。由于插值条件有四个，就将插值函数G 取为a + bx + cy + dxy（2.3.1）上式称为双一次式，其满足上述条件的插值问题称为矩形域上分片双一次插值。我们来寻找函数，它属于函数类G，又满足，或者，此式为（2.3.2）为此，我们将矩形各边的方程写出来：1+ y = 0 ：1 x = 0 ：1 y = 0 ：1+ x = 0现在，考虑函数的构造，从图1看出的对顶点是,过的两条邻边是，而就在由,组成的折线上,足见，若取为，（2.3.3）（1 x = 0，1 y = 0分别为，的方程）那么它属于函数类G，且满足，i=2,3,4。为了使，只要选取常数满足,所以。类似的方法可求得，列出如下， （2.3.4）或将其统一写成 （2.3.5）于是所求得分片双线性插值的表达式为（2.3.6）。（2.3.6）也可通过下列方式获得，先对u(x, y)作关于x的线性插值，记为,然后对作关于y的线性插值,记为，就是（2.3.6）的右端，现验明如下 ： 由关于x 的线性插值，再由关于y 的线性插值有 而这就是（2.3.6）的右端。2.3.2分片不完全双二次插值 给定矩形域上的四个顶点的函数值及四个边中点的函数值，要求作插值。由于条件有八个，选择函数类 (2.3.7)作插值函数，对于二元双二次式来讲，（2.3.7）式缺少项，故称 （2.3.7）为不完全的双二次插值。我们运用拉格朗日形式来构造插值函数，令, （2.3.8）这里 均属函数类（2.3.7）并且满足： B3 Y A4 A3 B2 O X B4 A1 A2 B1 图2 记，坐标为(见图2)矩形四边方程为：1+ y = 0 ：1 x = 0 ：1 y = 0 ：1+ x = 0矩形相邻二边中点联线的方程为：1+ x + y = 0 ：1 x + y = 0 ：1 x y = 0 ：1+ x y = 0现在我们来寻求函数，由于在线段在线段，上，及在直线上，可见应取为包括因子(1 x)(1 y)(1+ x + y)，这只要注意，及的方程就可得到，为了使，只要取类似的分析可得。最后我们将其统一写成，i=1,2,3,4 (2.3.9)注意到方程，就可将 取成,归纳之有,i=1,2,3,4 （2.3.10） 将（2.3.9），（2.3.10）代入（2.3.8）即为所求的插值多项式。2.3.3矩形域上分片双三次埃尔米特插值假设矩形单元如图 3，顶点坐标分别为 Y A4 A3 O X A1 A2 图3插值条件如下： i=1,2,3,4 (2.3.11)这个问题称为矩形域上的埃尔米特插值。插值条件有16 个，故可选择插值函数类为完全三次式 (2.3.12)为了寻求插值函数U(x, y)，先对y方向进行埃尔米特插值，记之为，然后对x方向进行埃尔米特插值，记之为，具体步骤如下： （2.3.13）其中 （2.3.14）对函数 , 进行x方向进行埃尔米特插值有（2.3.15）（2.3.15）中的。现在再引入记号 （2.3.16）由函数的性质可推知函数属于函数类（2.3.12）且有下列重要性质： （2.3.17）根据引入的记号（2.3.15）式可写成（2.3.18）由性质（2.3.17），容易看出上式右端就是满足插值条件（2.3.11）的完全三次式，或即有 （2.3.19）2.4二元函数插值程序设计 由上述分析可知，二元函数插值计算过程比较比较复杂，我们可以考虑设计程序用计算机实现：2.4.1MATLAB 中插值描述及程序设计 用函数 interp2 来拟合二维网格（X，Y）上的数据Z，语法是YI = INTERP2(X,Y, Z，XI, YI，方法)， 其中（X， Y，Z） 是已给的数据点， （XI，YI）是插值点坐标，其中方法主要有linear -线性插值，默认pchip -逐段三次Hermite 插值spline -逐段三次样条函数插值其中最后一种插值的曲面比较平滑例1：测得平板表面3*5网格点处的温度分别为： 82 81 80 82 84 79 63 61 65 81 84 84 82 85 86 试作出平板表面的温度分布曲面z=f(x,y)的图形。程序如下：clear;close;x=0:4;y=2:4;z=82 81 80 82 84;79 63 61 65 81;84 84 82 85 86;subplot(2,2,1);mesh(x,y,z);title(RAW DATA);xi=0:0.1:4;yi=2:0.1:4;zspline=interp2(x,y,z,xi,yi,spline);subplot(2,2,2);mesh(xi,yi,zspline);title(SPLINE);若数据是不规则的，即z 的数据不完全，不能构成一个矩阵，从而不能直接用interp2 插值。【ZI = griddata(x, y, z, xi, yi) x,y,z 均为向量（不必单调）表示数据，xi,yi 为网格向量，griddata 采用三角形线性插值ZI = griddata(x, y, z, xi, yi, cubic) 采用三角形三次插值】如果数据不完全如下残缺不全： yx01234234*7984*848061*8265*84*86则程序如下：x=2,3,4,0,2,3,0,1,4;y=2,2,2,3,3,3,4,4,4;z=80,82,84,79,61,65,84,84,86;subplot(2,2,3);stem3(x,y,z);title(RAW DATA);ZI=griddata(x,y,z,xi,yi,cubic); subplot(2,2,4);mesh(xi,yi,ZI);title(GRIDDATA)图形如下： 例2：设有数据x=1,2,3,4,5,6，y=1,2,3,4，在由x,y构成的网格上，数据为：12,10,11,11,13,1516,22,28,35,27,2018,21,26,32,28,2520,25,30,33,32,20求通过这些点的插值曲面。解：程序为：x=1:6;y=1:4;t=12,10,11,11,13,15 16,22,28,35,27,20 18,21,26,32,28,25; 20,25,30,33,32,20subplot(1,2,1)mesh(x,y,t)x1=1:0.1:6;y1=1:0.1:4;x2,y2=meshgrid(x1,y1);t1=interp2(x,y,t,x2,y2,cubic);subplot(1,2,2)mesh(x1,y1,t1); 针对这些基本方法及分片双一次插值，分片不完全的双二次插值和分片双三次埃尔米特插值的公式，设计了以下插值算法基本思想：针对分片双一次插值的设计Function U=nalagr(x,y,z,xx,yy)x=a1,a2,a3,a4;y=b1,b2,b3,b4;for i=1:4U(i)=(x(i),y(i);q(i)=(1+x(i)x)(1+y(i)y)/4;f(i)=U(i)*q(i);endU=f(1)+f(2)+f(3)+f(4);针对分片不完全的双二次插值的设计Function U=fchazhi(fh,x,y)t=a1,a2,a3,a4;p=b1,b2,b3,b4;A=c1,c2,c3,c4;B=d1,d2,d3,d4;for i=1:4U(i)=feval(fh,t(i),p(i);V(i)=feval(fh,A(i),B(i);q(i)=-(1+t(i)x)(1+p(i)y)(1-t(i)x-p(i)y)/4;s(i)=(1-p(i)x)(1+B(i)x)(1+t(i)y)(1+A(i)x)/2;f(i)=U(i)*q(i)+ V(i)* s(i);endU=f(1)+f(2)+f(3)+f(4);针对分片双三次埃尔米特插值的设计Function U=fchazhi(fh,x,y)t=a1,a2,a3,a4;p=b1,b2,b3,b4;A=difffh,x;B=difffh,y;C=diff(A,y); (1)=(1-x)2*(1+2*x)*(1-y)2*(1+2*y); (2)=x2*(3-2*x)*(1-y)2*(1+2*y); (3)=x2*(3-2*x)*(3-2*y); (4)=(1-x)2*(1+2*x)*y2*(3-2*y); (1)=x*(1-x)2*(1-y)2*(1+2*y); (2)=x2*(x-1)*(1-y)2*(1+2*y); (3)=x2*(x-1)*y2*(3-2*y); (4)=x*(1-x)2*y2*(3-2*y); (1)=(1-x)2*(1+2*x)*y*(1-y)2; (2)=x2*(3-2*x)*y*(1-y)2; (3)=x2*(3-2*x)*y2*(y-1); (4)=(1-x)2*(1+2*x)*y2*(y-1); (1)=x*(1-x)2*y*(1-y)2; (2)=x2*(x-1) *y*(1-y)2; (3)=x2*(x-1)*y2*(y-1); (4)=x*(1-x)2*y2*(y-1);for i=1:4D(i)=feval(fh,t(i),p(i);E(i)=feval(A, t(i),p(i);F(i)=feval(B, t(i),p(i);G(i)=feval(C, t(i),p(i);f(i)= D(i)* (i)+ E(i)* (i)+ F(i)* (i)+ G(i)* (i);endU=f(1)+f(2)+f(3)+f(4);3总述本文主要讲了二元函数插值的几种方法及其程序设计，以便解决实际中的问题，但在实行的过程中需要从一组数据去预测函数的表达式。从几何角度来说，这个问题就是要由给定的数据点去描绘曲线的近似图形。插值方法是处理这类问题的一种数值方法，不过由于插值曲线要求严格通过所给的每一个数据点，这会保留所给数据的误差，如果个别数据误差很大，那么插值效果显然不好。但是从足够多的数据当中，我们能发现规律即设法构造一条曲线（所谓拟合曲线）来反映所给数据点总的趋势，以消除局部的波动。二元函数插值就是解决的利器。最后，由于本人能力有限，论文中还有许多不完善的地方，但在我以后的学习过程中，我会努力的完善自己的理论知识、严格严谨的做学问！在这次的毕设中，让我掌握了matlab这个软件的应用，并且了解了二元函数插值的系统知识，进一步完善了自己的知识系统。更重要的是看到了自己的许多不足之处，让我懂得了面对困难要努力克服它，要知难而上！致谢毕业论文到此结束了，在此，首先感谢我的指导老师徐俊文老师对我整个毕业设计过程中的指导和帮助，使我能够顺利的完成我的学士学位论文。他严肃的科学态度，严谨的治学精神，精益求精的工作作风，深深地感染和激励着我。徐老师不仅在学业上给我以精心指导，同时还在思想、生活上给我以无微不至的关怀，在此谨向徐老师致以诚挚的谢意和崇高的敬意。参考文献1刘春凤，何亚丽.应用数值分析.冶金工业出版社，2005.1568 .2梁学章，李强.多元逼近.国防工业出版社，2005.2122 .3李岳生，黄友谦.数值逼近.人民教育出版社，1978.360367 .4（德）梅拉德斯 （Meinardus,G.）. 函数逼近.高等教育出版社 .5李庆扬，关治，白峰山.数值计算原理. 清华大学出版社 .6曹立凡，史万明.数值分析.北京工业学院出版社，1986.403404 .7姜健飞，胡良剑，唐俭.数值分析及其 MATLAB 实验. 科学出版社，2004.103 104 .英文摘要Dual function interpolation and its programmingWang guo-dong 20091105156Institute of mathematical sciences, information and computation sc