-拟合与插值课件

1 第四章第四章 拟合与插值拟合与插值2一、数据拟合：一、数据拟合：问题的提法：问题的提法：已知曲线上已知曲线上n个点个点(xi,yi),(i=1,2,n),xi互不相同，互不相同，寻求一个函数寻求一个函数y=F(x),使使y=F(x)在某种准则下与所有在某种准则下与所有数据点最为接近，即曲线拟合得最好。数据点最为接近，即曲线拟合得最好。基本思路是基本思路是:设设是是m+1个线性无关的连续函数，个线性无关的连续函数，称为基本拟合函数，取称为基本拟合函数，取 为待定系数，为待定系数，m n，拟合准则是拟合准则是:（最小二乘法）（最小二乘法）3按此方法得到的函数按此方法得到的函数称为数据集称为数据集(xi,yi),i=1,2,n的最小二乘拟合函数。的最小二乘拟合函数。衡量拟合情况优劣的量之一衡量拟合情况优劣的量之一 是：是：其中其中R2的值可作为拟合质量的指标。的值可作为拟合质量的指标。若它接近若它接近1，表示，表示F(x)能较好的拟合该数据集。能较好的拟合该数据集。4Matlab中的多项式拟合函数：中的多项式拟合函数：polyfit格式：格式：a=polyfit(x,y,n)功能：对一组点功能：对一组点(x,y)进行进行n次多项式拟合。次多项式拟合。x,y为要拟合的数据，为要拟合的数据，n为拟合多项式的次数，为拟合多项式的次数，a为拟合多项式为拟合多项式 系数，系数，若取：若取：即用即用m次多项式拟合给定数据，称多项式拟合。次多项式拟合给定数据，称多项式拟合。拟合多项式在拟合多项式在x处的值可用处的值可用polyval(a,x)计算。计算。可结合可结合polytool(x,y,n)处理多项式拟合问题。处理多项式拟合问题。5例例4.1 某通信公司在一次施工中，需要在一段水面宽某通信公司在一次施工中，需要在一段水面宽20米的河底沿直线走向铺设一条河底光缆，米的河底沿直线走向铺设一条河底光缆，A点表示点表示光缆入水处，光缆入水处，B点表示光缆出水处，点表示光缆出水处，C、D点表示河底点表示河底的光缆的两个触底点，在铺设前对河底地形做了测的光缆的两个触底点，在铺设前对河底地形做了测量，测得一组等分点位置量，测得一组等分点位置xi处的水深处的水深hi，x0位于点位于点A处，处，x20位于点位于点B处，数据如表（单位：米）处，数据如表（单位：米）试估算通过这条河沟所需光缆长度的近似值。试估算通过这条河沟所需光缆长度的近似值。ABCD6x0 x1x2x3x4x5x6x7x8x9x109.28.918.027.958.119.0510.1211.0812.2513.2613.45x11x12x13x14x15x16x17x18x19x2012.6211.2810.239.257.918.058.869.8510.8110.96河底深度测量数据（单位：米）河底深度测量数据（单位：米）ABCD7得出对河底光缆长得出对河底光缆长度的估计度的估计s，所需光缆为所需光缆为s+y(0)+y(20)对数据求出合适次数对数据求出合适次数的拟合多项式的拟合多项式p(x)，由弧长公式：由弧长公式：x=0:1:20;y=-9.2,-8.91,-8.02,-7.95,-8.11,-9.05,-10.12,-11.08,-12.25,-13.26,-13.45,-12.62,-11.28,-10.23,-9.25,-7.91,-8.05,-8.86,-9.85,-10.81,-10.96;polytool(x,y,3)n=input(拟合次数拟合次数n=)format longp=polyfit(x,y,n)dp=polyder(p)x1=0:0.01:20;y1=sqrt(1+polyval(dp,x1).2);s=trapz(x1,y1)-y(1)-y(end)6次多项式拟合较为次多项式拟合较为合适。合适。光缆总长：光缆总长：46.28m8二、一维插值二、一维插值 已知函数已知函数y=f(x)在区间在区间a,b上的上的n+1个不同点个不同点的函数值为的函数值为若存在一个简单函数若存在一个简单函数F(x),使使称称F(x)为为f(x)在区间在区间a,b上上的插值函数，的插值函数，称称(xi,yi)为插值节点为插值节点。当函数的表示式未知或非常复杂时，可通过插当函数的表示式未知或非常复杂时，可通过插值的方法研究函数。值的方法研究函数。若若F(x)为多项式，称为多项式插值为多项式，称为多项式插值(或代数插值或代数插值)；常用的代数插值方法有：拉格朗日插值，牛顿插值。常用的代数插值方法有：拉格朗日插值，牛顿插值。9 n次代数次代数插值：已知函数插值：已知函数f(x)在在n+1个点个点x0,x1,xn处处的函数值为的函数值为 y0,y1,yn，求一个求一个n次多项式函数次多项式函数Pn(x)，使其满足：使其满足：Pn(xi)=yi,(i=0,1,n).若若Pn(x)按下述方式构造，称为按下述方式构造，称为拉格朗日拉格朗日插值插值，称为称为拉格朗日插值基函数拉格朗日插值基函数其中其中Li(x)为为n次多项式：次多项式：10特别地特别地:（1 1）已知两个节点时，得线性插值多项式）已知两个节点时，得线性插值多项式:（2 2）已知三个节点时，得抛物插值多项式）已知三个节点时，得抛物插值多项式:已知已知n+1+1个节点时，可得个节点时，可得n次拉格朗日插值多项式。次拉格朗日插值多项式。11在区间在区间-1,1上上的代数插值的代数插值当插值次数过高时，常常会出现龙格现象，当插值次数过高时，常常会出现龙格现象，因此在数据建模中不建议使用过高次的代数插值。因此在数据建模中不建议使用过高次的代数插值。12若若F(x)为分段函数，称为分段插值；为分段函数，称为分段插值；插值中使用较多的是分段线性插值和三次样条插值。插值中使用较多的是分段线性插值和三次样条插值。1.分段线性插值分段线性插值 设函数设函数y=f(x)的插值节点为的插值节点为(xi,yi)，(i=0,1,2,n)将每两个相邻的节点用直线连起来，如此形成的一条将每两个相邻的节点用直线连起来，如此形成的一条折线就是分段线性插值函数记为：折线就是分段线性插值函数记为：In(x),它满足它满足In(xi)=yi,且在每个小区间上是一次函数。且在每个小区间上是一次函数。分段线性插值具有收敛性：分段线性插值具有收敛性：但在节点处但在节点处In(x)不可导，光滑性较差。不可导，光滑性较差。13 在在数数学学上上，光光滑滑程程度度的的定定量量描描述述是是：函函数数(曲曲线线)的的k阶导数连续，则称阶导数连续，则称该曲线具有该曲线具有k阶光滑性。阶光滑性。是否次数较低且是否次数较低且光滑性阶数较高光滑性阶数较高的分段多项式的分段多项式函数？函数？三次样条插值三次样条插值:且满足且满足:1)为三次多项式。为三次多项式。(自然边界条件自然边界条件)则则:且插值函数且插值函数S(x)具有具有2阶光滑性。阶光滑性。14Matlab中的一维插值函数中的一维插值函数interp1格式格式1：yi=interp1(x,y,xi)功能：已知函数功能：已知函数y=f(x)的的一组插值节点一组插值节点(x,y)，求出其，求出其插值函数插值函数F(x)在点在点xi处的函数值处的函数值yi。格式格式2：yi=interp1(y,xi)默认默认x=1:n,n为向量为向量y的元素个数。的元素个数。格式格式3：yi=interp1(x,y,xi,method)Method：(1)nearest 线性最近项插值线性最近项插值 (2)linear 线性插值线性插值 (3)spline 三次样条插值三次样条插值 (4)cubic 立方插值立方插值 缺省时为分段线性插值。缺省时为分段线性插值。15 x=0,3,5,7,9,11,12,13,14,15;y=0,1.2,1.65,2.1,2.15,2.0,1.85,1.65,1.55,1.25;xi=0:0.1:15;yi=interp1(x,y,xi,spline);yi1=interp1(x,y,xi,linear);yi2=interp1(x,y,xi,cubic);plot(x,y,*,xi,yi,r-,xi,yi1,b-,xi,yi2,g-)legend(节点节点,三次样条插值三次样条插值,线性插值线性插值,立方插值立方插值)例例4.2 通过实验测得某函数的一组数据如下：通过实验测得某函数的一组数据如下：x:0,3,5,7,9,11,12,13,14,15 y:0,1.2,1.65,2.1,2.15,2.0,1.85,1.65,1.55,1.25试作出其插值函数的图形。试作出其插值函数的图形。1617三次样条插值是解决一维插值问题最常用的方法，三次样条插值是解决一维插值问题最常用的方法，但是按但是按yi=interp1(x,y,xi,spline)只能得到希望求得只能得到希望求得插值点的函数值，而不能得到插值函数的表达式。插值点的函数值，而不能得到插值函数的表达式。若要得到插值函数的表达式，可用若要得到插值函数的表达式，可用csape函数函数。x=0,3,5,7,9,11,12,13,14,15;y=0,1.2,1.65,2.1,2.15,2.0,1.85,1.65,1.55,1.25;xi=0:0.1:15;S=csape(x,y)P=S.coefs yi=ppval(S,xi);plot(x,y,o,xi,yi,r);18得：得：P=0.0039 -0.0716 0.5800 0 0.0110 -0.0368 0.2548 1.2000 -0.0181 0.0290 0.2392 1.6500 0.0112 -0.0793 0.1386 2.1000 -0.0019 -0.0119 -0.0437 2.1500 -0.0129 -0.0232 -0.1139 2.0000 0.0608 -0.0618 -0.1990 1.8500 -0.0802 0.1205 -0.1403 1.6500 -0.0399 -0.1201 -0.1399 1.550019三、三、二维插值二维插值1.网格节点数据插值函数：网格节点数据插值函数：interp2格式：格式：z=interp2(x0,y0,z0,x,y,method)x0,y0,z0:插值节点坐标，要求插值节点坐标，要求x0,y0单调；单调；x,y是被插值点的横坐标与纵坐标是被插值点的横坐标与纵坐标(x,y不能超过不能超过x0,y0的的范围范围)，z是被插值点的函数值。是被插值点的函数值。Method：(1)nearest 最邻近插值最邻近插值 (2)linear 双线性插值双线性插值 (3)cubic双三次插值双三次插值默认为双线性插值。默认为双线性插值。20例例4.3 要在一山区修建公路，首先测得一些点的高要在一山区修建公路，首先测得一些点的高程（见附件程（见附件,表中数据为坐标点的高程，单位：米，表中数据为坐标点的高程，单位：米，y轴正向为北轴正向为北)绘制该地地貌图。绘制该地地貌图。480048001350135013701370139013901400140014101410960960940940880880800800690690570570430430290290210210150150440044001370137013901390141014101430143014401440114011401110111010501050950950820820690690540540380380300300210210400040001380138014101410143014301450145014701470132013201280128012001200108010809409407807806206204504503703703503503600360014201420143014301450145014801480150015001550155015101510143014301300130012001200980980850850750750550550500500320032001430143014501450146014601500150015501550160016001550155016001600160016001600160015501550150015001500150015501550150015002800280095095011901190137013701500150012001200110011001550155016001600155015501380138010701070900900105010501150115012001200240024009109101090109012701270150015001200120011001100135013501450145012001200115011501010101088088010001000105010501100110020002000880880106010601230123013901390150015001500150014001400900900110011001060106095095087087090090093093095095016001600830830980980118011801320132014501450142014201400140013001300700700900900850850840840380380780780750750120012007407408808801080108011301130125012501280128012301230104010409009005005007007007807807507506506505505508008006506507607608808809709701020102010501050120012008308308008007007003003005005005505504804803503504004005105106206207307308008008508508708708508507807807207206506505005002002003003003503503203200 0370370470470550550600600670670690690670670620620580580450450400400300300100100150150250250y/xy/x0 040040080080012001200160016002000200024002400280028003200320036003600400040004400440048004800520052005600560021先将数据保存为先将数据保存为data403数据文件数据文件 clear,clc load data403 x=0:400:5600;y=4800:-400:0;X,Y=meshgrid(x,y);surf(X,Y,A)pause figure(2)xi=linspace(0,5600,80);yi=linspace(0,4800,80);Xi,Yi=meshgrid(xi,yi);Zi=interp2(X,Y,A,Xi,Yi);surf(Xi,Yi,Zi)pause 22 figure(3)t=0:100:1600;c,h=contourf(Xi,Yi,Zi,t);clabel(c,h)colormap cool colorbar232425262.散点数据插值函数：散点数据插值函数：griddata已知已知n个插值节点个插值节点(xi,yi,zi),(i=1,2,n),求在点求在点(x,y)处的插处的插值值z，matlab提供函数提供函数griddata。格式：格式：cz=griddata(x,y,z,cx,cy,method)其中其中x,y,z 均为均为n 维向量，指明所给数据点（插值节点）维向量，指明所给数据点（插值节点）的横坐标、纵坐标和竖坐标。的横坐标、纵坐标和竖坐标。cx,cy是给定被插值点的是给定被插值点的横坐标和纵坐标，横坐标和纵坐标，cz为相应点的竖坐标。为相应点的竖坐标。若若cx,cy是向量，则给定以它们所确定网格点的横坐标和是向量，则给定以它们所确定网格点的横坐标和纵坐标，这时要求纵坐标，这时要求cx,cy一个为行向量一个为列向量。一个为行向量一个为列向量。编程时也可先用编程时也可先用meshgrid将将cx,cy定义成网格矩阵。定义成网格矩阵。27例例 4.4 在某海域测得一些点在某海域测得一些点(x,y)处的水深处的水深z（英尺）（英尺）如下表，船的吃水深度为如下表，船的吃水深度为5英尺，在矩形区域英尺，在矩形区域(75,200)(-50,150)里那些地方船要避免进入。里那些地方船要避免进入。x129140103.588185.5195105157.5107.57781162162117.5y7.5141.52314722.5137.585.5-6.5-81356.5-66.584-33.5z4868688998894928 clear,clcx=129,140,103.5,88,185.5,195,105,157.5,107.5,77,81,162,162,117.5;y=7.5,141.5,23,147,22.5,137.5,85.5,-6.5,-81,3,56.5,-66.5,84,-33.5;z=-4,8,6,8,6,8,8,9,9,8,8,9,4,9;x1,y1=meshgrid(75:5:200,-50:5:150);z1=griddata(x,y,z,x1,y1,v4);surf(x1,y1,z1)figure(2)c,h=contourf(x1,y1,z1);clabel(c,h)293031四、模型中的参数识别四、模型中的参数识别 由试验数据，确定模型中的未知参数。由试验数据，确定模型中的未知参数。可用：可用：fit函数，函数，lsqcurvefit函数，函数，lsqnonlin函数函数例例4.5 在某次阻尼振荡实验中测得在某次阻尼振荡实验中测得18组数据点，组数据点，试确定其振动方程。试确定其振动方程。t：0 0.4 1.2 2 2.8 3.6 4.4 5.2 6 y:1 0.85 0.29 -0.27-0.53 -0.4 -0.12 0.17 0.28 t：7.2 8 9.2 10.4 11.6 12.4 13.6 14.4 15 y:0.15 -0.03 -0.15 -0.07 0.06 0.08 0.032 -0.015 -0.02 由物理背景知：由物理背景知：32 clc,clearx=0,0.4,1.2,2,2.8,3.6,4.4,5.2,6,7.2,8,9.2,10.4,11.6,12.4,13.6,14.4,15;y=1,0.85,0.29,-0.27,-0.53,-0.4,-0.12,0.17,0.28,0.15,-0.03,-0.15,-0.07,0.059,0.08,0.032,-0.015,-0.02;f=fittype(a*cos(k*x)*exp(w*x),independent,x,coeff,a,k,w);cfun=fit(x,y,f)xi=0:0.1:20;yi=cfun(xi);plot(x,y,r*,xi,yi,b-)33 运行结果：运行结果：General model:cfun(x)=a*cos(k*x)*exp(w*x)Coefficients(with 95%confidence bounds):a=0.9988 (0.9837,1.014)k=1.001 (0.9959,1.006)w=-0.2067 (-0.2132,-0.2003)3435练习题练习题1.海水温度随着深度的变化而变化，海面温度较高，随海水温度随着深度的变化而变化，海面温度较高，随着深度的增加温度越来越低，这样也就影响了海水的对着深度的增加温度越来越低，这样也就影响了海水的对流和混合，使得深层海水中的氧气越来越少，根据这一流和混合，使得深层海水中的氧气越来越少，根据这一规律可对海水鱼层作一划分，现通过实验观测得一组海规律可对海水鱼层作一划分，现通过实验观测得一组海水深度水深度h与温度与温度t的数据如下：的数据如下：t23.522.420.919.115.411.59.58.88.48.1h01.52.54.68.212.516.520.222.526.2(1)找出深度与温度之间的函数找出深度与温度之间的函数t=t(h)。(2)找出温度变化最快的深度位置。找出温度变化最快的深度位置。362.在对山区农作物种植面积进行评估时，需对某一区域在对山区农作物种植面积进行评估时，需对某一区域的地表面积进行估计，现测得一个矩形区域内一些地面的地表面积进行估计，现测得一个矩形区域内一些地面点的高程点的高程(见附表见附表)，假设所观测的矩形区域内均是陆地，假设所观测的矩形区域内均是陆地且整个地表可近似看做一张光滑曲面。且整个地表可近似看做一张光滑曲面。要求：要求：（1）通过二维插值绘出该观察区域地表的曲面图）通过二维插值绘出该观察区域地表的曲面图（2）估算出观察区域地表面积的近似值）估算出观察区域地表面积的近似值（3）计算从矩形区域一边上任意一点到其对边上任意）计算从矩形区域一边上任意一点到其对边上任意一点沿直线走向的地表曲线的近似长度（自己在矩形的一点沿直线走向的地表曲线的近似长度（自己在矩形的对边上取两个定点求出）。对边上取两个定点求出）。37 xy01002003004005006007008009001000110012000381.5397.3401.5396.1379.2361.5339.6320.7307.8299.6293.5285.5273.0100398.0418.5423.2416401380.5356.3336.4325.3315.5310.4306.3298.5200411.2428.3436.5433.5420.5398.5374.5353.4338.3330.2325.4323.1317.0300420.2439.6448.8443.1435.5415.2392.3370.2354.1346.5342.3338.2332.5400423.8444.6456.7459.2447.6429.2407.3385.5369.4360.3356.5352.5345.2500422.1441.5459.5463.2454.6439.5419.5400.5384.5374.5369.5364.5356.4600416.2439.2454.5466.5457.6446.5429.5412.4398.5386.3381.2375.5362.1700406.3429.5442.3454.6450.2448.1435.2422.2409.3399.5391.2380.2370.5800391.5414.6431.4443.9447.8440.2436.5426.5416.5406.3397.4382.2363.5900373.5395.6413.5426.6433.5434.5430.5424.5416.5407.6390.0388.6370.01000352.3374.5391.6405.5414.3418.5416.8413.6407.5398.5389.5376.5350.5 附表：区域附表：区域0,1200;0,1000内观测点处的高程表内观测点处的高程表（单位：米）（单位：米）