插值法与数据拟合法

第页第七讲插值方法及数据拟合§7.1引言在工程和科学实验中，常常需要从一组实验观测数据(xi,yi)(i=1,2,…,n)揭示自变量x及因变量y之间的关系，一般可以用一个近似的函数关系式y=f(x)来表示。函数f(x)的产生办法因观测数据及要求的不同而异，通常可采用两种方法：插值及数据拟合。§7.1.1插值方法引例1已经测得在北纬32.3海洋不同深度处的温度如下表：表7.1.1深度x(m)46671495014221634水温y(C)7.044.283.402.542.13根据这些数据，我们希望能合理地估计出其它深度（如500米、600米、1000米…）处的水温。解决这个问题，可以通过构造一个及给定数据相适应的函数来解决，这是一个被称为插值的问题。2．插值问题的基本提法对于给定的函数表xx0x1…xny=f(x)y0y1…yn其中f(x)在区间[a,b]上连续，x0，x1，…，xn为[a,b]上n1个互不相同的点，要求在一个性质优良、便于计算的函数类{P(x)}中，选出一个使P(xi)=yi，i=0,1,…,n(7.1.1)成立的函数P(x)作为f(x)的近似，这就是最基本的插值问题（见图7.1.1）。为便于叙述，通常称区间[a,b]为插值区间，称点x0，x1，…，xn为插值节点，称函数类{P(x)}为插值函数类，称式(7.1.1)为插值条件，称函数P(x)为插值函数，称f(x)为被插函数。求插值函数P(x)的方法称为插值法。§7.1.2数据拟合引例2在某化学反应中，已知生成物的浓度及时间有关。今测得一组数据如下：表7.1.2时间t（分）12345678浓度y1034.006.408.008.809.229.509.709.86时间t（分）9116浓度y10310.0010.2010.3210.3210.5010.5510.5810.60根据这些数据，我们希望寻找一个y=f(t)的近似表达式（如建立浓度y及时间t之间的经验公式等）。从几何上看，就是希望根据给定的一组点（1,4.00），…,（16,10.60），求函数y=f(t)的图象的一条拟合曲线。数据拟合问题的基本提法对于给定的函数表xx0x1…xny=f(x)y0y1…yn其中f(x)在区间[a,b]上连续，x0，x1，…，xn为[a,b]上n1个互不相同的点，要求找一个简单合理的函数近似表达式(x)，使(x)及f(x)在某种准则下最为接近，这就是最基本的数据拟合问题（见图7.1.2）。通常，我们称(x)为给定数据点的拟合函数。图7.1.1插值问题示意图图7.1.2数据拟合问题示意图§7.1.3插值方法及数据拟合的基本理论依据插值方法及数据拟合的基本理论依据，就是数学分析中的Weierstrass定理：设函数f(x)在区间[a,b]上连续，则对>0，存在多项式P(x)，使得即：有界区间上的连续函数被多项式一致逼近。§7.1.4实际应用中两种方法的选择在实际应用中，究竟选择哪种方法比较恰当？总的原则是根据实际问题的特点来决定采用哪一种方法。具体说来，可从以下两方面来考虑：1．如果给定的数据是少量的且被认为是严格精确的，那么宜选择插值方法。采用插值方法可以保证插值函数及被插函数在插值节点处完全相等。2．如果给定的数据是大量的测试或统计的结果，并不是必须严格遵守的，而是起定性地控制作用的，那么宜选用数据拟合的方法。这是因为，一方面测试或统计数据本身往往带有测量误差，如果要求所得的函数及所给数据完全吻合，就会使所求函数保留着原有的测量误差；另一方面，测试或统计数据通常很多，如果采用插值方法，不仅计算麻烦，而且逼近效果往往较差。§7.2一维数据的基本插值方法简介插值函数类的取法很多，可以是代数多项式，也可以是三角多项式或有理函数；可以是[a,b]上任意光滑函数，也可以是分段光滑函数。在此介绍最基本、最常用的两种插值方法：分段多项式插值及三次样条插值，及其Matlab实现。§7.2.1一维数据的分段多项式插值对于给定的一维数据xx0x1…xny=f(x)y0y1…yn分段多项式插值就是求一个分段（共n段）多项式P(x)，使其满足P(xi)=yi（i=0,1,…,n）或更高的要求。一般地，分段多项式插值中的多项式都是低次多项式（不超过三次）。分段线性插值y分段线性插值函数P1(x)是一个分段一次多项式（分段线f(x)性函数）。在几何上就是用折线代替曲线，如图7.2.1，故分段线性插值亦称为折线插值。其插值公式为P(x)，x[xi,xi+1](7.2.1)0x0x1xn-1xnx2．分段二次插值图7.2.1分段线性插值示意图分段二次插值函数P2(x)是一个分段二次多项式。在几何上就是分段抛物线代替曲线y=f(x)，故分段二次插值又称为分段抛物插值。其插值公式为，x[xi-1,xi+1](7.2.2)3．三次Hermite插值三次Hermite插值问题的基本提法一：已知一维数据xx0x1y=f(x)y0y1y=f(x)m0m1求一个三次多项式P3(x)，使之满足P3(xi)=yi，P3(xi)=mi，i=0,1(7.2.3)构造三次插值基函数0(x)，1(x)，0(x)，1(x)，使之满足(7.2.4)利用这四个插值基函数，取三次多项式P3(x)为P3(x)=0(x)y01(x)y10(x)m01(x)m1(7.2.5)将插值条件(7.2.3)式代入，可推得：(7.2.6)(7.2.5)、(7.2.6)两式构成了三次Hermite插值基本提法一的插值公式。三次Hermite插值问题的基本提法二：已知一维数据xx0x1x2y=f(x)y0y1y2y=f(x)m1求一个三次多项式P3(x)，使之满足P3(xi)=yi，i=0,1,2，P3(x1)=mi(7.2.7)构造三次插值基函数0(x)，1(x)，2(x)，1(x)，使之满足(7.2.8)利用这四个插值基函数，取三次多项式P3(x)为P3(x)=0(x)y01(x)y12(x)y21(x)m1(7.2.9)将插值条件(7.2.7)式代入，可推得：(7.2.10)(7.2.9)、(7.2.10)两式构成了三次Hermite插值基本提法二的插值公式。§7.2.2一维数据的三次样条插值上述介绍的分段多项式插值，其优点为计算简单、稳定性好、收敛性有保证，且易于在计算机上实现。但它也明显存在着缺陷。它只能保证在每个小区间段[xi,xi+1]内光滑，在各小区间连接点xi处连续，却不能保证整条曲线的光滑、光顺性，难以满足某些工程的要求。对于象高速飞机的机翼形线，船体放样等型值线往往要求有二阶光滑度，即有二阶连续导数。而由60年代开始，首先起源及航空、造船业等工程设计的实际需要而发展起来的样条插值，既保留了分段多项式插值的各种优点，又提高了插值函数的光滑度。在此，仅介绍应用最广且具有二阶连续导数的三次样条插值方法。三次样条插值问题的基本提法对于给定的一维数据xx0x1…Xny=f(x)y0y1…Yn求一个三次多项式S(x)满足条件（1）S(xi)=yi，i=0,1,…,n；（2）S(x)具有二阶连续导数，特别在节点xi上应满足连续性要求，即对i=0,1,…,n有三次样条插值函数给定区间[a,b]的一个划分：a=x0<x1<…<xn=b，设函数y=f(x)在节点xi上的值为yi=f(xi)，i=0,1,…,n。如果S(x)于[a,b]有二阶连续导数，且在每个小区间[xi,xi+1]上是三次多项式，则称S(x)是节点x0，x1，…，xn上的三次样条函数。如果S(x)在节点xi上还满足插值条件S(xi)=yi，i=0,1,…,n，(7.2.11)则称S(x)为三次样条插值函数。对应于划分的三次样条插值函数的表达式为(7.2.12)其中，。边界条件在式(7.2.12)给出的三次多项式中，共含有n3个待定系数。而由插值条件(7.2.11)式，可列出n1个方程，方程组中未知数的个数比方程个数多2，还需附加2个条件才能进行求解。通常可在区间端点x0=a和xn=b处各附加一个条件（称为边界条件或边值条件）去确定S(x)。边界条件类型很多，较基本而又常见的有三类：第一边值条件，即给出边界点的一阶导数值S(x0)=y0，S(xn)=yn(7.2.13)第二边值条件，即给出边界点的二阶导数值S(x0)=y0，S(xn)=yn(7.2.14)特别地，当S(x0)=S(xn)=0时，称为自然边界条件。满足自然边界条件的三次样条插值函数称为自然样条插值函数。第三边值条件（混合边值条件）(7.2.15)其中1、2、1、2、1、2为定数。当1、2为零时，则为第一边值条件，当1、2为零时，则为第二边值条件。§7.2.3一维数据插值的Matlab实现1．一维数据插值的Matlab实现Matlab中一维数据的插值函数为：interp1()。其调用格式为yi=interp1(x,y,xi,'methos'),其中x，y——为插值节点，均为向量；xi——任取的被插值点可以是一个数值，也可以是一个向量；yi——为被插值点xi处的插值结果；'methos'——为采用的插值方法：'nearest'：表示最临近插值，'linear'：表示线性插值，'cubic'：表示三次插值，'spline'：表示三次样条插值。注意：（1）上述'methos'中所有的插值方法都要求x是单调的，并且xi不能超过x的取值范围；（2）三次样条插值函数的调用格式有两种，yi=interp1(x,y,xi,'spline')和yi=spline(x,y,xi)，它们是等价的。2．例：§7.1.1引例1的求解解法一直接借助于Matlab命令，分别采用最临近插值方法、线性插值方法、三次插值方法和三次样条插值方法进行插值。Matlab程序如下：x=[466,714,950,1422,1634];y=[7.04,4.28,3.40,2.54,2.13];x1=linspace(400,1660,100);y1=interp1(x,y,x1,'nearest');y2=interp1(x,y,x1,'linear');y3=interp1(x,y,x1,'cubic');y4=interp1(x,y,x1,'spline');subplot(2,2,1),plot(x,y,'o',x1,y1),title('最邻近插值曲线')subplot(2,2,2),plot(x,y,'o',x1,y2),title('线性插值曲线')subplot(2,2,3),plot(x,y,'o',x1,y3),title('三次插值曲线')subplot(2,2,4),plot(x,y,'o',x1,y4),title('三次样条插值曲线')x2=500:100:1600w1=interp1(x,y,x2,'nearest')w2=interp1(x,y,x2,'linear')w3=interp1(x,y,x2,'cubic')w4=interp1(x,y,x2,'spline')要想求得任一深度值xi处的水温，只需在上述程序运行完毕之后，继续运行如下命令即可：xi=500:100:1600yi1=interp1(x,y,xi,'nearest')yi2=interp1(x,y,xi,'linear')yi3=interp1(x,y,xi,'cubic')yi4=interp1(x,y,xi,'spline')下面是运行上述程序后得到的四种不同插值方法在500米到1600米不同深度处的水温值（见表7.2.1），以及四种不同的插值函数的曲线图（见图7.2.2）。从图中可以看出，三次插值和三次样条插值得到的插值曲线的光滑性更要好些。表7.2.1深度(m)59001000最邻近插值的水温(C)7.04004.28004.28004.28003.40003.4000线性插值的水温(C)6.66165.54874.43583.95933.58643.3089三次插值的水温(C)6.48875.21444.37033.84333.52033.2905三次样条插值的水温(C)6.48875.21444.37033.84333.52033.2905深度(m)110012001600最邻近插值的水温(C)3.40002.54002.54002.54002.54002.1300线性插值的水温(C)3.12672.94452.76232.58012.38922.1958三次插值的水温(C)3.09252.91422.74602.57792.40032.2032三次样条插值的水温(C)3.09252.91422.74602.57792.40032.2032图7.2.2四种不同插值函数的曲线图如果要想求得任一深度值xi处的水温，只需在上述程序运行完毕之后，继续运行如下命令即可：xi=500:100:1600yi1=interp1(x,y,xi,'nearest')yi2=interp1(x,y,xi,'linear')yi3=interp1(x,y,xi,'cubic')yi4=interp1(x,y,xi,'spline')解法二求三次样条插值函数选取第一边值条件，并用差商和分别代替y0和yn，则所求三次样条插值函数为(7.2.16)且满足（见表7.1.1）：S3(xi)=yi，(i=0,1,…,4)，S3(x0)=0.0111，S3(x4)=0.0019(7.2.17)其中，。将已知数据代入(7.2.17)式，得到一个7阶线性方程组。解之得（用alpha=linsolve(a,b)命令）：0=3.8854，1=0.0363，2=1.7650104，3=3.2126107，1=5.4789107，2=2.2710107，3=4.5398109，即0=3.8854，1=0.0363，2=0.0002，3=1=2=3=0（用alpha=a\b命令也可得到这个结果）。并计算得在500米到1600米不同深度处的水温值如下：表7.2.2深度(m)59001000水温(C)6.64925.44114.39523.77463.48823.3157深度(m)110012001600水温(C)3.14372.96572.77992.58452.37912.1763Matlab程序如下：x=[466,714,950,1422,1634];y=[7.04,4.28,3.40,2.54,2.13];a=[1,x(1),(x(1)^2)/2,(x(1)^3)/6,0,0,0;...1,x(2),(x(2)^2)/2,(x(2)^3)/6,0,0,0;...1,x(3),(x(3)^2)/2,(x(3)^3)/6,((x(3)-x(2))^3)/6,0,0;...1,x(4),(x(4)^2)/2,(x(4)^3)/6,((x(4)-x(2))^3)/6,((x(4)-x(3))^3)/6,0;...1,x(5),(x(5)^2)/2,(x(5)^3)/6,((x(5)-x(2))^3)/6,((x(5)-x(3))^3)/6,((x(5)-x(4))^3)/6;...0,1,x(1),(x(1)^2)/2,0,0,0;...0,1,x(5),(x(5)^2)/2,((x(5)-x(2))^2)/2,((x(5)-x(3))^2)/2,((x(5)-x(4))^2)/2];b=[7.04;4.28;3.40;2.54;2.13;-0.0111;-0.0019];alpha=a\bx1=500:100:1600;fori=1:length(x1)ifx1(i)>=466&x1(i)<=714a1=[1,x1(i),(x1(i)^2)/2,(x1(i)^3)/6,0,0,0];y1(i)=a1*alpha;elseifx1(i)>714&x1(i)<=950a3=[1,x1(i),(x1(i)^2)/2,(x1(i)^3)/6,((x1(i)-x(2))^3)/6,0,0];y1(i)=a3*alpha;elseifx1(i)>950&x1(i)<=1422a4=[1,x1(i),(x1(i)^2)/2,(x1(i)^3)/6,((x1(i)-x(2))^3)/6,((x1(i)-x(3))^3)/6,0];y1(i)=a4*alpha;elsea5=[1,x1(i),(x1(i)^2)/2,(x1(i)^3)/6,((x1(i)-x(2))^3)/6,((x1(i)-x(3))^3)/6,((x1(i)-x(4))^3)/6];y1(i)=a5*alpha;endendy1§7.3二维数据的基本插值方法简介对于二维数据的插值，首先要考虑两个问题：一是二维区域是任意区域还是规则区域，二是给定的数据是有规律分布的还是散乱的、随机分布的。第一个问题比较容易处理。目前的插值方法基本上是基于规则区域的，对于不规则区域，只需将其，划分为规则区域或扩充为规则区域来讨论即可。对于第二个问题，当给定的数据是有规律分布时，方法较多也较成熟；而给定的数据是散乱的、随机分布时，没有固定的方法，但一般的处理思想是：从给定的数据出发，依据一定的规律恢复出规则分布点上的数据，转化为数据分布有规律的情形来处理。二维数据插值的方法也有很多。在此，针对给定数据有规律分布和散乱分布两种情形，简单介绍双三次样条插值方法和改进的Shepard方法（反距离平方法）的基本概念和基本思想，及其Matlab实现。§7.3.1双三次样条插值双三次样条插值方法，是用来解决规则区域上给定数据有规律分布的插值问题的常用方法。设R：[a,b][c,d]是xy平面上的一个矩形区域。在x轴和y轴上分别取定分割x：a=x0<x1<…<xn<b，y：c=y0<y1<…<ym<d。由此导出R上的一个矩形网格分割：xy，将R分成mn个子矩形Rij：[xi1,xi][yi1,yi]。它的两条邻边长分别是hi=xixi1(i=1,2,…,n)，gj=yjyj1(j=1,2,…,m)。直线x=xi(i=1,2,…,n)和y=yj(j=1,2,…,m)称为分割的两族网格线。网格线的交点(xi,yj)(i=1,2,…,n，j=1,2,…,m)称为分割的节点，总共有(m1)(n1)个节点，见图7.3.1。ydRyiRijyi1c0axi1xib图7.3.1矩形分割设函数S(x,y)定义在矩形域R上且满足下列条件：（1）在每个子矩形Rij(i=1,2,…,n，j=1,2,…,m)上，S(x,y)是一个关于x和y的都是三次的多项式函数，即；(7.3.1)（2）在整个R上，函数S(x,y)的偏导数都是连续的。则称S(x,y)为双三次样条函数。如果给定一数组{sij}(i=1,2,…,n，j=1,2,…,m)，双三次样条函数S(x,y)还满足插值条件S(xi,yi)=sij(i=1,2,…,n，j=1,2,…,m)，就称S(x,y)为插值双三次样条函数。实际上，双三次样条函数是由两个一维三次样条函数作直积产生的。对任意固定的y0[c,d]，S(x,y0)是关于x的三次样条函数，同理，对任意固定的x0[a,b]，S(x0,y)是关于y的三次样条函数。从而，根据一维三次样条函数的算法可以设计出S(x,y)的具体算法。§7.3.2改进的Shepard方法改进的Shepard方法，也称反距离加权平均法，这是解决规则区域上给定数据散乱、随机分布的插值问题的一个常用的方法。问题：设T=[a,b][c,d]上散乱分布N个点V1，V2，…，VN，其中Vk=(xk，yk)处给出数据fk，k=1,2,…,N。要寻求T上的二元函数F(x，y)，使F(xk，yk)=fk，k=1,2,…,N。一个典型的容易想到的方法是“反距离加权平均”方法，又称之为Shepard方法。这方法的基本思想是，在非给定数据的点处，定义其函数值由已知数据点及该点距离的近或远作加权平均决定。记，则二元函数（曲面）定义为，其中(7.3.2)如此定义的曲面是全局相关的，对曲面的任一点作数据计算都要涉及到全体数据，这在大量实测数据插值和拟合中是很慢的。此外，F(x，y)在每个插值点(xk，yk)附近产生一个小的“平台”，使曲面不具有光顺性。当fk=c（常数）时，F(x，y)=c。当(xk，yk)皆位于非常数的同一平面上时，F(x，y)并不能保证再生这个平面。但因为这种作法思想简单，人们对它作了种种改进，例如，取适当的常数R>0，令由于()是可微函数，使得如下定义的F(x，y)在性能上有所改善，其中。(7.3.3)按照上述的思想，可从给定的数据恢复出规则分布点上的数据，接下来就可应用双三次样条插值或其它的二维数据插值方法来处理。§7.3.3二维数据插值的Matlab实现1．规则区域上给定数据有规律分布的二维插值数据形式为：y1y2…ynx1x11z12…z1nx2z21z22…z2n……………xmzm1zm2…zmn插值函数为：interp2()。其调用格式为zi=interp2(x,y,z,xi,yi,'methos'),其中x，y，z——为插值节点，均为向量；zi——为被插值点(xi,yi)处的插值结果；'methos'——为采用的插值方法：'nearest'：表示最临近插值，'linear'：表示双线性插值，'cubic'：表示双三次插值，'spline'：表示双三次样条插值。注意：上述'methos'中所有的插值方法都要求x和y是单调的网格，x和y可以是等距的也可以是不等距的。2．规则区域上给定数据散乱或随机分布的二维插值数据形式为：(x1,y1)(x2,y2)…(xn,yn)z1z2…zn插值函数为：e01sef和e01sff，。通常两者配合使用，其调用格式为[fnodes,a,rnw,b,c]=e01sef(x,y,z)；[pf(I,j),ifail]=e01sff(x,y,z,rnw,fnodes,px(j),py(i))；其中x，y，z——为插值节点，均为向量；px(j)，py(i)——为被插值节点；pf(i,j)——为被插值点(px(j),py(i))处的插值结果；其它输出参数涉及插值算法，可以不用了解。e01sef的输出fnodes和rnw为确定插值的参数，它们是e01sff需要的输入参数，因此两函数需配合使用。3．例例1气旋变化情况的可视化表7.3.1是气象学家测量得到的气象资料，它们分别表示在南半球地区按不同纬度、不同月份的平均气旋数字。根据这些数据，绘制出气旋分布曲面图形。表7.3.10—1010—2020—3030—4040—5050—6060—7070—8080—901月2.418.720.822.137.348.225.65.30.32月1.621.418.520.128.836.624.25.303月2.416.218.220.527.835.525.55.404月3.29.216.625.137.24024.64.90.35月1.02.812.929.240.337.621.14.906月0.51.710.132.641.735.422.27.107月0.41.48.333.046.23520.25.30.18月0.22.411.231.039.934.721.27.30.28月0.55.812.528.625.935.722.670.310月0.89.221.132.040.339.528.58.6011月2.410.323.928.138.24025.36.30.112月3.61625.525.643.441.924.36.60.3解下面分别用最邻近插值、双线性插值、双三次插值和双三次样条插值，给出不同月份按纬度变化的气旋值（插值结果），并作出可视化图形如下。图7.3.2四种插值方法的可视化图形最邻近插值的Matlab程序为：x=1:12;y=5:10:85;z=[2.4,18.7,20.8,22.1,37.3,48.2,25.6,5.3,0.31.6,21.4,18.5,20.1,28.8,36.6,24.2,5.3,02.4,16.2,18.2,20.5,27.8,35.5,25.5,5.4,03.2,9.2,16.6,25.1,37.2,40,24.6,4.9,0.31.0,2.8,12.9,29.2,40.3,37.6,21.1,4.9,00.5,1.7,10.1,32.6,41.7,35.4,22.2,7.1,00.4,1.4,8.3,33.0,46.2,35,20.2,5.3,0.10.2,2.4,11.2,31.0,39.9,34.7,21.2,7.3,0.20.5,5.8,12.5,28.6,25.9,35.7,22.6,7,0.30.8,9.2,21.1,32.0,40.3,39.5,28.5,8.6,02.4,10.3,23.9,28.1,38.2,40,25.3,6.3,0.13.6,16,25.5,25.6,43.4,41.9,24.3,6.6,0.3];[xi,yi]=meshgrid(1:12,5:1:85);zi=interp2(x,y,z,xi,yi,'nearest');figuresurf(xi,yi,zi)xlabel('月份'),ylabel('纬度'),zlabel('气旋'),axis([012090050])title('南半球气旋可视化图形')双线性插值、双三次插值、双三次样条插值的Matlab程序为：分别将最邻近线性插值程序中的zi=interp2(x,y,z,xi,yi,'nearest')改写为zi=interp2(x,y,z,xi,yi,'linear')zi=interp2(x,y,z,xi,yi,'cubic')zi=interp2(x,y,z,xi,yi,'spline')例2水道测量数据（AMCM86A—题）在某海域测得一些点（x,y）处的水深z（单位：英尺）由表7.3.2给出，水深数据是在低潮时测得的。船的吃水深度为5英尺，问在矩形区域(75，200)(50，150)里的哪些地方船要避免进入。表7.3.2水道水深测量数据（单位：英尺）x129.0140.0108.588.0185.5195.0105.5y7.5141.528.0147.022.5137.585.5z4868688x157.5107.577.081.0162.0162.0117.5y6.581.03.056.566.584.038.5z9988949解（1）假设由题目给出的信息是很少的，除了14个位置的水深之外一无所知。显然，题目要求我们找出水深不到5英尺的区域。为了讨论方便，下面三个假设是合理的：=1\*GB3①所给数据是精确的；=2\*GB3②讨论区域的海底曲面是光滑的，更确切地说，可以认为曲面的一阶、二阶导数是连续的。因为我们可以认为讨论区域为浅水海域，由于长期的海水水流作用，形成的是以砾石或沙为主要组成部分的海底，不存在珊瑚礁、水底峡谷、山脊等不可意料的突变地形。=3\*GB3③水深是一个按区域来划分的变量，在某个位置的水深及其周围区域的水深是相互依赖的，但这种依赖作用随距离的增大而减小。就我们讨论的问题来说，每一个给定数据点影响周围的每一个未知点，一个给定数据点离未知点越近，作用就越大。（2）问题分析根据假设，海底曲面是连续光滑的，不存在珊瑚礁、水底峡谷、山脊等不可意料的突变地形，因而很自然的想法就是用某种光滑的拟合曲面去逼近已知的14个数据点或以14个已知的数据点为基础，利用二维插值补充一些点的水深，以求得水深不超过5米的区域。在此，我们采用二维插值方法，应用Matlab程序，作出矩形区域(75，200)(50，150)范围内的海底地形图、水深不超过5米的危险区域的平面图以及水深不超过5米的危险区域的海底地貌图，并求出水深不超过5米的危险海域范围。（3）问题求解采用改进的Shepard方法，利用Matlab软件作出作出矩形区域(75，200)(50，150)范围内的海底地形图、水深不超过5米的危险区域的平面图以及水深不超过5米的危险区域的海底地貌图（见图7.3.3），并求出水深不超过5米的危险海域范围为：[115，200][3，119]。（4）求解的Matlat程序如下：clear;x=[129,140,108.5,88,185.5,195,105.5,157.5,107.5,77,81,162,162,117.5];y=[7.5,141.5,28,147,22.5,137.5,85.5,-6.5,-81,3,56.5,-66.5,84,-38.5];subplot(2,2,1),plot(x,y,'+'),title('测量点分布图');z=[-4,-8,-6,-8,-6,-8,-8,-9,-9,-8,-8,-9,-4,-9];[fnodes,minnq,rnw,rnq,ifail]=e01sef(x,y,z);nx=100;px=linspace(75,200,nx);图图7.3.3ny=200;py=linspace(-50,150,ny);fori=1:nyforj=1:nx[pf(i,j),ifail]=e01sff(x,y,z,rnw,fnodes,px(j),py(i));endendsubplot(2,2,2),meshz(px,py,pf+5),title('(75，200)x(-50，150)范围内的海底地形图');[a,b]=find(pf>=-5);amin=min(a);amax=max(a);bmin=min(b);bmax=max(b);xmin=75+((200-75)/100)*bminxmax=75+((200-75)/100)*bmaxymin=-50+((150+50)/200)*aminymax=-50+((150+50)/200)*amaxfork=1:length(b)i0(k)=75+((200-75)/100)*b(k);endfork=1:length(a)j0(k)=-50+((150+50)/200)*a(k);endsubplot(2,2,3),plot(i0,j0,'+'),title('水深不超过5米的危险区域的平面图');[i1,j1]=find(pf<-5);fork=1:length(i1)pf(i1(k),j1(k))=-5;endsubplot(2,2,4),meshc(px,py,pf),title('水深不超过5米的危险区域的海底地貌图');rotate3d§7.4一维数据拟合的最小二乘法简介数据拟合问题的形式多种多样，解决的方法也有许多。在此，我们只简单介绍以最小二乘法为准则的一维数据拟合方法。§7.4.1数据拟合的最小二乘法简介对于给定的一组测量数据xx1x2…xny=f(x)y1y2…yn设y=(x)为一拟合函数，记i=(xi)yi(i=1,2,…,n)，(7.4.1)则称i为拟合函数(x)在xi点处的偏差或残量。为使(x)在整体上尽可能及给定数据的函数f(x)近似，我们常采用偏差的平方和达到最小，即(7.4.2)来保证每个偏差的绝对值|i|都很小。这一原则称为最小二乘原则，根据最小二乘原则确定拟合函数(x)的方法称为最小二乘法。线性最小二乘拟合我们知道，函数系{xk|k=0,1,…,m}的线性组合(x)=a0a1xa2x2…amxm为m次多项式。一般地，若函数系{k(x)|k=0,1,…,m}是线性无关的，则其线性组合(7.4.3)称为函数系{k(x)|k=0,1,…,m}的广义多项式。如三角多项式就是函数系{1,cosx,sinx,cos2x,sin2x,…,cosmx,sinmx}的广义多项式。设{k(x)|k=0,1,…,m}为一线性无关的函数系，取拟合函数为(7.4.3)式给出的广义多项式，使得(7.4.2)成立。由于(x)的待定系数a0,a1,a2,…,am全部以线性形式出现，故我们称之为线性最小二乘拟合。在式(7.4.2)中，目标函数S是关于参数a0,a1,a2,…,am的多元函数，由多元函数取得最小值的必要条件知，欲使S达到极小，须满足，（k=0,1,…,m）即亦即，（k=0,1,…,m）(7.4.4)式(7.4.4)是关于a0,a1,a2,…,am的线性方程组，称为正规方程组。从正规方程组(7.4.4)中解出a0,a1,a2,…,am，于是就得到了最小二乘拟合函数(x)。非线性最小二乘拟合如果拟合函数(x)=(x,a0,a1,a2,…,am)的待定参数a0,a1,a2,…,am不能全部以线性形式出现，如指数拟合函数等，这便是非线性最小二乘拟合问题。一般地，非线性最小二乘拟合问题是一个非线性函数的极小化问题，可用非线性优化方法求解。最小二乘拟合函数的选择最小二乘法中，拟合函数的选择是很重要的。可以通过对给定数据的分析来选择，也可以直接由实际问题给定。最常用的是多项式和样条函数，尤其是当不知道该选择什么样的拟合函数时，通常可以考虑选择样条函数。另外，对同一问题，也可选择不同的函数进行最小二乘拟合，比较各自误差的大小，从中选出误差较小的作为拟合函数。§7.4.2一维数据拟合的Matlab实现多项式函数拟合拟合函数的命令为：polyfit()，其调用格式为a=polyfit(xdata,ydata,m),其中m——表示多项式的最高阶数；xdata，ydata——为将要拟合的数据，它们都是以数组方式输入；a——输出参数，为拟合多项式的系数a=[a0,a1,a2,…,am]。多项式在x处的值y可用如下命令格式计算：y=polyval(a,x)。一般的曲线拟合拟合函数的命令为：curvefit()，或lsqcurvefit()，其调用格式为p=curvefit('Fun',p0,xdata,ydata)，或p=lsqcurvefit('Fun',p0,xdata,ydata)，其中Fun——表示函数Fun(p,xdta)的M文件；P0——为函数的初值。若要求点x处的函数值y，可用程序f=Fun(p,x)计算。例：§7.1.2引例2的求解解：（1）选择拟合曲线作出给定数据的散点图如下：图7.4.1给定数据的散点图通过对散点图的分析可以看出，数据点的分布为一条单调上升的曲线。具有这种特性的曲线很多，我们选取如下三种函数来拟合：=1\*GB3①多项式(x)=a0a1x…amxm，m为适当选取的正整数；=2\*GB3②；=3\*GB3③（a>0,b<0）。（2）拟合运算首先，分别用二、三、六次多项式拟合，计算得输出参数分别为p1=[0.0445，1.0711，4.3252]p2=[0.0060，0.1963，2.1346，2.5952]p3=[0.0000，0.0004，0.0103，0.1449，1.1395，4.9604，0.0498]拟合函数分别为(1)(x)=0.04451.0711x4.3252x2(2)(x)=0.00600.1963x2.1346x22.5952x3(3)(x)=0.0004x0.0103x20.1449x31.1395x44.9304x50.0498x6；其次，再用有理分式拟合，计算得输出参数分别为p=[0.0841，0.1392]拟合函数为最后，用指数函数拟合，计算得输出参数分别为p=[11.3578，1.0873]拟合函数为三种方式五个种函数的拟合曲线见图7.4.2—7.4.4。图7.4.2多项式函数拟合曲线图图7.4.3有理分式函数拟合曲线图图7.4.4指数函数拟合曲线图（3）误差分析和给定的16组数据比较，三种方式五个函数拟合的误差见下表：表7.4.1五个函数拟合的误差表偏差平方和平均偏差平方和最大偏差二次多项式拟合4.44160.27761.3518三次多项式拟合1.22920.07680.6067六次多项式拟合0.11690.00730.2684有理分式拟合0.57320.03580.4772指数函数拟合0.17770.01110.2544其中：偏差平方和为，平均偏差平方和为。从上表中求得的误差情况来看，好似六次多项式函数拟合的最好，指数函数拟合次之，然后分别是有理分式函数拟合、三次多项式函数拟合和二次多项式函数拟合。但是，就这个实际问题的本质来说，化学反应中生成物的浓度到一定时间后应基本稳定，即当t时，f(t)常数。而我们有并且三个多项式函数拟合曲线的趋势也可从图7.4.5看出。因而，我们有理由认为指数函数拟合的效果最好，有理分式函数拟合的效果次之，本问题不宜使用多项式拟合。（4）结论通过以上的计算和分析，我们得出如下结论：本问题可用指数函数或有理分式函数来拟合，其拟合函数分别为和拟合误差为表7.4.2偏差平方和平均偏差平方和最大偏差指数函数拟合0.17770.01110.4772有理分式拟合0.57320.03580.2544图7.4.5（5）Matlab程序=1\*GB3①多项式拟合的Matlab程序：t=1:16;y=[4,6.4,8,8.4,9.28,9.5,9.7,9.86,10,10.2,10.32,10.42,10.5,10.55,10.58,10.6];figure(1)%作出给定数据的散点图plot(t,y,'o')p1=polyfit(t,y,2)p2=polyfit(t,y,3)p3=polyfit(t,y,6)x1=linspace(0,16,160);y1=polyval(p1,x1);y2=polyval(p2,x1);y3=polyval(p3,x1);figure(2)%作出三条拟合曲线图plot(t,y,'.',x1,y1,'-',x1,y2,':',x1,y3,'-.')u1=polyval(p1,t);%以下计算误差u2=polyval(p2,t);u3=polyval(p3,t);d1=(u1-y).^2;d2=(u2-y).^2;d3=(u3-y).^2;delta11=sum(d1)delta21=sum(d2)delta31=sum(d3)delta12=sum(d1)/16delta22=sum(d2)/16delta32=sum(d3)/16delta13=sqrt(abs(max(d1)))delta23=sqrt(abs(max(d2)))delta33=sqrt(abs(max(d3)))x2=linspace(0,25,60);%画三条拟合曲线的趋势图w1=polyval(p1,x2);w2=polyval(p2,x2);w3=polyval(p3,x2);figure(3)plot(t,y,'.',x2,w1,'-',x2,w2,':',x2,w3,'-.')title('三条拟合曲线的趋势图')=2\*GB3②有理分式拟合的Matlab程序：tdata=1:16;ydata=[4,6.4,8,8.4,9.28,9.5,9.7,9.86,10,10.2,10.32,10.42,10.5,10.55,10.58,10.6];figure(1)plot(tdata,ydata,'o')p0=[0.1,0.1];cde:\matlab\model\%input'p=curvefit('e111fun',p0,tdata,ydata)fori=1:16phi(i)=p(1)*exp(p(2)/tdata(i));endfigure(2)plot(tdata,phi)figure(3)plot(tdata,ydata,'.',tdata,phi)d=(phi-ydata).^2;delta1=sum(d)delta2=sum(d)/16delta3=abs(max(d))=3\*GB3③指数函数拟合的Matlab程序：tdata=1:16;ydata=[4,6.4,8,8.4,9.28,9.5,9.7,9.86,10,10.2,10.32,10.42,10.5,10.55,10.58,10.6];figure(1)plot(tdata,ydata,'o')p0=[0.1,0.1];cde:\matlab\model\%input'p=curvefit('e110fun',p0,tdata,ydata)fori=1:16phi(i)=tdata(i)/(p(1)*tdata(i)+p(2));endfigure(2)plot(tdata,phi)figure(3)plot(tdata,ydata,'.',tdata,phi)title('有理分式拟合')d=(phi-ydata).^2;delta1=sum(d)delta2=sum(d)/16delta3=abs(max(d))=4\*GB3④调用的M文件程序：e110fun.m文件functionphi=e90fun(p,tdata)fori=1:16phi(i)=tdata(i)/(p(1)*tdata(i)+p(2));endphi;e111fun.m文件functionphi=e91fun(p,tdata)fori=1:16phi(i)=p(1)*exp(p(2)/tdata(i));endphi;参考文献：易大义等，数值方法，浙江科学技术出版社，1984。徐涛，数值计算方法，吉林科学技术出版社，2019。寿纪麟，数学建模——方法及范例，西安交通大学出版社，1993。叶其孝主编，大学生数学建模竞赛辅导教材，湖南教育出版社，1993。李尚志主编，数学建模竞赛教程，江苏教育出版社，1996。傅鹏等编著，数学实验，科学出版社，2000。谢云荪主编，数学实验，科学出版社，2000。应用Matlab求解水道测量数据问题摘要：水道测量数据问题是一个给定数据散乱、随机分布的二维离散数据的插值问题。本文以水道测量数据问题为例，应用Matlab软件提供的求解三维网格点数据的函数，对求解决给定数据散乱、随机分布的二维离散数据插值问题，给出了一个简便易行的方法。问题重述水道测量数据问题是1986年美国大学生数学建模竞赛的A题，由加州海军研究生院数学系的RichardFranke提供。问题如下：在某海域测得一些点（x,y）处的水深z（单位：英尺）由表1给出，水深数据是在低潮时测得的。船的吃水深度为5英尺，问在矩形区域(75，200)(50，150)里的哪些地方船要避免进入。表1水道水深测量数据（单位：英尺）x129.0140.0108.588.0185.5195.0105.5y7.5141.528.0147.022.5137.585.5z4868688x157.5107.577.0