第二章 插值與擬合

第二章插值與擬合2.1插值與擬合的基本概念2.1.1插值與插值函數(shù)已知由g3)(可能未知或非常復(fù)雜)產(chǎn)生的一批離散數(shù)據(jù)(土,七)，i=°」，，n，且n+1個互異插值節(jié)點(diǎn)a=x<x< <x<x=b，在插值區(qū)間內(nèi)尋找一個相對簡單的函0 1 n-1 n ???數(shù)f(x)，使其滿足下列插值條件：了氏)Z="…部再利用已求得的f(x)計(jì)算任一非插值節(jié)點(diǎn)x*的近似值y*=f(x*)，這就是插值。其中f(x)稱為插值函數(shù)，g(x)稱為被插函數(shù)。下面介紹幾種常用的而且有現(xiàn)成MATLAB命令的插值方法的數(shù)學(xué)原理。分段線性插值將兩個相鄰節(jié)點(diǎn)用直線連起來，如此形成的一條折線就是分段線性插值函數(shù)，記作I(x)，它滿足I(x)=y，且I(x)在每個小區(qū)間[x,x],j=0,1,,n，上是線性函數(shù)。n njj n jj+1I(x)可以表示為TOC\o"1-5"\h\zn …I(x)=Wyl(x)njjj=0x—x ..\o"CurrentDocument" 1,x.<x<x.,j=0舍去,j j—1+x)+x)=1x—x j+1-x—xjj+10x<x<x ,j=n舍去,其他I〃(x)有良好的收斂性，即對于xG[a,b]時，有l(wèi)imI〃(x)=g⑴。n—8用In(x)計(jì)算x點(diǎn)的插值時,只用到x左右的兩個節(jié)點(diǎn)，計(jì)算量與節(jié)點(diǎn)個數(shù)n無關(guān)。但是n越大，分段越多，插值誤差越小。MATLAB中有現(xiàn)成的分段線性插值命令：y=interp1(x0,y0,x)其中x0,y0為節(jié)點(diǎn)數(shù)組(同長度)，x為插值點(diǎn)數(shù)組，y為插值數(shù)組。三次樣條插值三次樣條函數(shù)記作S(x)(a<x<b)，要求它滿足以下條件：在每個小區(qū)間[x「x.](i=1,,n)上是三次多項(xiàng)式；在a<x<b上二階導(dǎo)數(shù)連續(xù)；S(x)=y.,(i=0,1,,n)?！蓷l件(1)，不妨將S(x)記為???S(x)={S(x),xg[x,x](i=1,,n)}，S(x)=ax3+bx2+cx+d，i i—1i i ii ii其中a,b,c,d為待定系數(shù)，共4n個。由條件(2)，iiii ???S(土)=Si+1(」Sf(x)=S'(x),i=1,2,,n-1ii i+1iS"(x)=S"(x).ii i+1i ,,,上式與條件(3)共有4n-2個方程，為確定S(x)的4n個待定參數(shù)，還需兩個條件。在實(shí)際應(yīng)用中通常有以下三種類型的端點(diǎn)條件作為附加條件。第一類給定兩端點(diǎn)的一階導(dǎo)數(shù)S'(x°),Sf(x^)；第二類給定兩端點(diǎn)的二階導(dǎo)數(shù)S”(；°),S"")，最常用的是所謂的自然邊界條件：S”(x)=S"(x)=0； "0 n第三類對于周期函數(shù)，即兩端點(diǎn)已經(jīng)滿足S(x)=S(x)時，令它們的一階導(dǎo)數(shù)及二0n階導(dǎo)數(shù)分別相等，即S'(x0)=S'(xn)，S"(x0)=S"(xn)，稱為周期條件。這樣，就構(gòu)成了4n元線性方程組，可以證明有蜂一解，即S(x)被唯一確定。對于三次樣條插值，MATLAB中有現(xiàn)成的命令y=interp1(x0,y0,x,’spline’)或者y=spline(x0,y0,x)其中x0,y0為節(jié)點(diǎn)數(shù)組(同長度)，x為插值點(diǎn)數(shù)組，y為插值數(shù)組，端點(diǎn)為自然邊界條件。Spline命令還可以處理上述第一類端點(diǎn)條件，只需將輸入數(shù)組y0改為yy0=[ay0b]，其中a，b分別為S'(x0),S'(x.)。高維插值 ° ”MATLAB還給出了高維插值函數(shù)interpN()，其中N可以為2,3，…，例如N=2時，為二維插值，調(diào)用格式為Zi=interp2(x,y,z,xi,yi,’method’)其中x,y,z為插值節(jié)點(diǎn)，zi為被插值點(diǎn)(xi,yi)處的插值結(jié)果。’method’表示采用的插值方法：’nearest’表示最鄰近插值；’linear’表示線性插值；’cubic’表示三次插值。所有插值方法都要求x,y是單調(diào)的網(wǎng)格。2.1.2 最小二乘擬合已知一批離散的數(shù)據(jù)(xi,七),i=0,1,,n，xi互不相同，尋求一個擬合函數(shù)y=f(x)，使f(x)與y的誤差平方和在最小二乘意義下最小。在最小二乘意義下確定的ii ???f(x)稱為最小二乘擬合函數(shù)。一元最小二乘法給定平面上的點(diǎn)(x,,七),i=0,1,,n，進(jìn)行曲線擬合有多種方法，最小二乘法是解決曲線擬合最常用的一種方法。最小二乘法的提法是：求f(x)，使8=X?2=X[f(x)—y.]2達(dá)到最小。i=1 i=1擬合時選用一定形式的擬合函數(shù)，擬合函數(shù)可由一些簡單的“基函數(shù)”(例如冪函數(shù)，三角函數(shù)等等)覽(x),氣(x),,氣(x)來線性表示：f(x)=c中(x)+c中(x)+ +c中(x)現(xiàn)在要確定系數(shù)c0,c1,,cm，使S達(dá)到嚶小。"為此，將f(x)的表達(dá)式代入5中，8就成TOC\o"1-5"\h\z為c,c,,c的函數(shù)，令5對c(i=1,2,,m)的偏導(dǎo)數(shù)等于零，于是得到m+1個方程組，0 1m ??? i由此求出c(i=1,2,,m)。通常取基函數(shù)為1,x,x2,x3,,xm，這時擬合函數(shù)f(x)為多項(xiàng)?.?i ???式函數(shù)。特別地，當(dāng)m=1時，f(x)=a+bx，稱為一元線性擬合函數(shù)。??? ???MATLAB中提供了多項(xiàng)式擬合的語句：a=polyfit(x,y,m)x,y為要擬合的數(shù)據(jù)，是長度自定義的數(shù)組，m為擬合多項(xiàng)式的次數(shù)。輸出參數(shù)為擬合多項(xiàng)式的系數(shù)。多項(xiàng)式為降冪形式：y=axm+axm-1+ +ax+a輸出系數(shù)次序a=[a「a2, ,a^]。多項(xiàng)式在x處的擬合值「可用下面程序計(jì)算：???y=polyval(a,x)如何選擇擬合函數(shù)已知一組數(shù)據(jù)(七,七),i=0,1,,n，選擇什么樣的函數(shù)f⑴呢？一是根據(jù)機(jī)理分析來確定函數(shù)形式，二是根據(jù)散點(diǎn)圖直觀判斷函數(shù)f(x)的形式。要比較兩個模型哪個擬合效果更佳，則比較兩個模型的殘差平方和。殘差平方和較小者更佳。設(shè)y為擬合函數(shù)的值，y為測量值，則殘差e=￡(y-y)2i i iii曲面擬合簡介實(shí)際問題中可能遇到曲面擬合問題，可以將一元最小二乘方法的有關(guān)概念和結(jié)論推廣到多元最小二乘方法。已知m個自變量(x「 ，xm)和一個因變量y的一組觀測值(x,x,,x,y),i=1,2,,n，要確定函數(shù)y=f(x, ,x)，使得TOC\o"1-5"\h\z1i2imii ??? 1m\o"CurrentDocument"minJ=Y[f(x,x..,,x)-y]2 ...1i2i miii=1一般地，首先通過機(jī)理分析或數(shù)據(jù)的直觀判斷，去確定函數(shù)y=f(x1, ,xm)的結(jié)構(gòu)，假定函數(shù)中含有未知參數(shù)a,a,,a，然后，通過最小二乘原理具體確定參數(shù)a,a,,a。1 2k ??- 1 2k2.1.3溫度預(yù)測問題??????在12小時內(nèi)，每隔1小時測量一次溫度。溫度依次為：5,8,9，15,25,29,31，30，22，25，27，24。(單位:°C)試分別用分段線性插值、三次樣條插值方法估計(jì)在3.2h，6.5h，7.1h，11.7h的溫度值，每隔1/10h估計(jì)一次溫度值并畫出其圖形。用多項(xiàng)式擬合，估計(jì)在3.2h，6.5h，7.1h，11.7h的溫度值。hours=1:12;temps=[5,8,9,15,25,29,31,30,22,25,27,24];t=interp1(hours,temps,[3.2,6.5,7.1,11.7]) %線性插值T=interp1(hours,temps,[3.2,6.5,7.1,11.7],'spline') %三次樣條插值計(jì)算結(jié)果為t=10.2000 30.0000 30.9000 24.9000T=9.6734 30.0427 31.1755 25.3820每隔1/10h估計(jì)一次溫度值并畫出其圖形：hours=1:12;temps=[5,8,9,15,25,29,31,30,22,25,27,24];h=1:0.1:12;t=interp1(hours,temps,h,'spline');plot(hours,temps,'+',h,t,hours,temps,'r:')xlabel('時間,),ylabel(溫度')

圖2.1溫度插值圖形用三次多項(xiàng)式擬合：hours=1:12;temps=[5,8,9,15,25,29,31,30,22,25,27,24];a=polyfit（hours,temps,3）temps1=polyval（a,hours）;plot（hours,temps,'ro',hours,temps1,'b.'）得到y(tǒng)=—0.0065x3—0.3283x2+7.1281x-4.4343，圖形如下：計(jì)算在3.2h,6.5h,7.1h,11.7h的溫度值為：14.8,26.24,27.30,23.61。2.1.4山區(qū)地貌圖在某山區(qū)（平面區(qū)域（0，2800）x（0，2400）內(nèi)，單位：m）測得一些地點(diǎn)的高度（單位：m）如表1.1所示，試作出該山區(qū)的地貌圖。表1.1某山區(qū)一些地點(diǎn)的高度z040080012001600200024002800240014301450147013201280120010809402000145014801500155015101430130012001600146015001550160015501600160016001200137015001200110015501600155013808001270150012001100135014501200115040012301390150015001400900110010600118013201450142014001300700900用二維三次插值方法，可以得到不同方向上的山體的高度（插值結(jié)果），然后在matlab中建立M-file程序如下：y=0:400:2400;x=0:400:2800;z=[118013201450142014001300700900;1230139015001500140090011001060;12701500120011001350145012001150;13701500120011001550160015501380;14601500155016001550160016001600;14501480150015501510143013001200;1430145014701320128012001080940;];[xi,yi]=meshgrid(0:20:2800,0:20:2400);zi=interp2(x,y,z,xi,yi,'cubic');surf(xi,yi,zi);xlabel('x軸'),ylabel('y軸'),zlabel('z軸')；axis([028000240002000]);title('山區(qū)地貌圖')；2.2行駛汽車車距問題美國的某些司機(jī)培訓(xùn)課程中有這樣的規(guī)則：正常駕駛條件下車速每增加16公里/小時，后面與前面一輛車的距離應(yīng)增加一個車身的長度。又云，實(shí)現(xiàn)這個規(guī)則的一種簡便辦法是所謂2秒準(zhǔn)則。即后車司機(jī)從前車經(jīng)過某一標(biāo)志開始默數(shù)2秒鐘后到達(dá)同一標(biāo)志，而不管車速如何。試判籪2秒準(zhǔn)則與上述規(guī)則是一樣的么？這個規(guī)則的合理性如何，是否有更好的規(guī)則。下面表2-1是測得的車速和剎車距離的一組數(shù)據(jù)表2-1車速和剎車距離的一組數(shù)據(jù)車速（km/h）20406080100120140剎車距離（m）6.517.833.657.183.4118.0153.51.問題分析制訂這樣的規(guī)則是為了在后車急剎車情況下不致撞上前車，即要確定汽車的剎車距離。剎車距離顯然與車速有關(guān)，先計(jì)算出汽車以16km/h的車速在2秒鐘行駛的距離，這個距離為：16公里/小時X（1小時/3600秒）X2秒=8.9米，而汽車平均長度4.6m，所以“2秒準(zhǔn)則”與上述規(guī)則并不一致。要判斷規(guī)則是否合理，則要對剎車距離進(jìn)行分析。剎車距離由反應(yīng)距離和制動距離兩部分組成，反應(yīng)距離指司機(jī)看到需要剎車的情況到汽車制動器開始起作用汽車行駛的距離，制動距離指制動器開始起作用到汽車完全停止行駛的距離。反應(yīng)距離由反應(yīng)時間和車速決定，反應(yīng)時間取決于司機(jī)個人狀況（靈敏、機(jī)警等）和制動系統(tǒng)的靈敏性，由于很難對反應(yīng)時間進(jìn)行區(qū)別，因此，通常認(rèn)為反應(yīng)時間為常數(shù)，而且在這段時間內(nèi)車速不變。制動距離與制動作用力、車重、車速以及路面狀況等因素有關(guān)。由能量守恒制動力所作的功被汽車動能的改變所抵消。設(shè)計(jì)制動器的一個合理原則是，最大制動力大體上與汽車質(zhì)量成正比，使汽車的減速度基本上是常數(shù)。路面狀況可認(rèn)為是固定的。模型假設(shè)基于以上分析，作如下假設(shè):剎車距離d等于反應(yīng)距離d1和制動距離d2之和。反應(yīng)距離di與車速v成正比，比例系數(shù)為反應(yīng)時間七。剎車時使用最大制動力F，F(xiàn)作的功等于汽車動能的改變，且F與車的質(zhì)量m成正比。模型的建立由假設(shè)2d=tv由假設(shè)3,在F作用下行駛距離d2作的功Fd2使車速從v變成0，動能的變化為mv2/2，有Fd2=mv2/2，又F憂m，按照牛頓第二定理F=ma可知，剎車時的減速度a為常數(shù)，于是d=kv其中k為比例系數(shù)，實(shí)際上k=1/2a。由假設(shè)1，剎車距離為d=tvkv1下面對參數(shù)t和k進(jìn)行估計(jì)。通常有經(jīng)驗(yàn)估計(jì)和數(shù)據(jù)擬合兩種方法，這里我們采用數(shù)據(jù)1擬合的方法。如果沒有首先推導(dǎo)出模型形式的話，也可直接用多項(xiàng)式擬合。由上述推導(dǎo)得出，實(shí)際剎車距離的擬合多項(xiàng)式為d=k1V2+k2V。其M-file程序如下：v=[20:20:140]/3.6;v2=v.A2;x=[v;v2]';d=[6.5,17.8,33.6,57.1,83.4,118.0,153.5]';a=x\ddd=x*addd=[6.5,17.8,33.6,57.1,83.4,118.0,153.5];b=polyfit(v,ddd,2)y=polyval(b,v)plot(v,ddd,'ro',v,dd,'b')運(yùn)行后得到k1=0.0853,k2=0.6522，即d=0.6522v+0.0853v2如果不知道函數(shù)的形式可以作出實(shí)際距離的2次多項(xiàng)式擬合擬合，d=k1v2+k2v+k3其M-file程序如下：v=[20:20:140]/3.6;d=[6.5,17.8,33.6,57.1,83.4,118.0,153.5];

k=polyfit(v,d,2)%插值d1=polyval(k,v)%擬合plot(v,d1,'b',v,d,'r*');圖2.2函數(shù)d=kiv2+k2v擬合圖 二次多項(xiàng)式d=kiv2+k2v+k3擬合圖圖2.2給出了實(shí)際剎車距離與擬合曲線的計(jì)算剎車距離的比較，對應(yīng)數(shù)據(jù)列在表2-2中。模型應(yīng)用由d=0.6522v+0.0853v2，根據(jù)假設(shè)及模型的建立知道，反應(yīng)時間匕=0.6522，如果設(shè)制動時間為，，如果已知速度則由物理學(xué)勻變速運(yùn)動的運(yùn)動方程可得到相應(yīng)速度的制動距2離，d2=vt2-+at;，其中a=或=5.86。再考察2秒準(zhǔn)則的合理性，根據(jù)車距和汽車車速，即可計(jì)算出汽車行駛一段車距路程的時間，而車距就取為擬合的剎車距離（當(dāng)然如果有最大剎車距離數(shù)據(jù)，則利用最大剎車距離計(jì)算會更好）。這個時間我們稱為“車距時間”，部分結(jié)果列在表2-2中。根據(jù)計(jì)算結(jié)果，可以將“車距時間”2秒規(guī)則改為：車速0~19（km/h）,“車距時間”為1秒；車速20~59（km/h）,“車距時間”為2秒；車速60~99（km/h）,“車距時間”為3秒；車速100~140（km/h）,“車距時間”為4秒。表2-2車速與剎車距離和計(jì)算剎車距離及車距時間車速（km/h）20406080100120140實(shí)際剎車距離（m）6.517.833.657.183.4118.0153.5計(jì)算剎車距離（m）6.2617.7834.5656.6283.93116.52154.37車距時間（s）1.11.62.12.53.03.53.72.3國土面積的數(shù)值計(jì)算1.問題的提出已知?dú)W洲一個國家的的地圖，為了算出它的國土面積，首先對地圖作如下測量：以由西向東方向?yàn)閄軸，由南向北方向?yàn)閥軸，選擇方便的原點(diǎn)，并將從最西邊界點(diǎn)到最東邊界點(diǎn)

在X軸上的區(qū)間適當(dāng)?shù)姆譃槿舾啥危诿總€分點(diǎn)的y方向測出南邊界點(diǎn)和北邊界點(diǎn)的y坐標(biāo)約和y2，這樣就得到了表2-3的數(shù)據(jù)（單位：mm）。1根據(jù)地圖的比例我們知道18mm相當(dāng)于40km，試由測量數(shù)據(jù)計(jì)算該國國土的近似面積，與它的精確值41288km2比較。表2-3 某國國土地圖邊界測量值（單位：mm）X7.010.513.017.534.040.544.548.056.0y444547505038303034y4459707293100110110110X61.068.576.580.591.096.0101.0104.0106.5y363441454643373328y117118116118118121124121121X111.5118.0123.5136.5142.0146.0150.0157.0158.0y326555545250666668y1211221168381828685682.模型的假設(shè)我們作如下假設(shè)（1） 假設(shè)測量的地圖和數(shù)據(jù)準(zhǔn)確，由最西邊界點(diǎn)與最東邊界點(diǎn)分為上下兩條連續(xù)的邊界曲線，邊界內(nèi)的所有土地均為該國國土。（2） 假設(shè)我們從最西邊界點(diǎn)到最東邊界點(diǎn)，變量xe［a,b］，劃分［a,b］為n小段［七,氣］,并由此將國土分為n小塊，設(shè)每一小塊均為X型區(qū)域。即作垂直于X軸的直線穿過該區(qū)域，直線與邊界曲線最多只有兩個交點(diǎn)。3.數(shù)值積分方法計(jì)算國土面積根據(jù)測量數(shù)據(jù)我們利用Matlab軟件對上下邊界進(jìn)行三次多項(xiàng)式插值，得到圖2.3。數(shù)值積分法的基本思想是將上邊界點(diǎn)與下邊界點(diǎn)分別利用插值函數(shù)求出兩條曲線，則曲線所圍面積即為國土面積（地圖上的國土面積），然后根據(jù)比例縮放關(guān)系求出國土面積的近似解。在求國土面積時，利用求平面圖形面積的數(shù)值積分方法一一將該面積分成若干個小長方形，分別求出長方形的面積后相加即為該面積的近似解。設(shè)上邊界函數(shù)為f（X），下邊界函數(shù)為f（X），由定積分定義可知曲線所圍區(qū)域