用MATLAB進(jìn)行數(shù)據(jù)插值資料

我們經(jīng)常會遇到大量的數(shù)據(jù)需要處理，而處理數(shù)據(jù)的關(guān)鍵就在于這些算法，例如數(shù)據(jù)擬合、參數(shù)估計(jì)、插值等數(shù)據(jù)處理算法。此類問題在MATLAB中有很多現(xiàn)成的函數(shù)可以調(diào)用，熟悉MATLAB，這些方法都能游刃有余的用好。一、概述第1頁/共107頁第一頁，共108頁。數(shù)據(jù)擬合在很多賽題中有應(yīng)用，與圖形處理有關(guān)的問題很多與插值和擬合有關(guān)系，例如98年美國賽A題，生物組織切片的三維插值處理，94年A題逢山開路，山體海拔高度的插值計(jì)算，2003年吵的沸沸揚(yáng)揚(yáng)的“非典”問題也要用到數(shù)據(jù)擬合算法，觀察數(shù)據(jù)的走向進(jìn)行處理，2005年的雨量預(yù)報(bào)的評價(jià)的插值計(jì)算。2001年的公交車調(diào)度擬合問題，2003年的飲酒駕車擬合問題。第2頁/共107頁第二頁，共108頁。插值問題——雨量預(yù)報(bào)的評價(jià)預(yù)測點(diǎn)和實(shí)測點(diǎn)的圖形插值后的圖形第3頁/共107頁第三頁，共108頁。擬合問題——飲酒駕車喝兩瓶酒的擬合曲線喝1-5瓶酒的擬合曲線第4頁/共107頁第四頁，共108頁。在實(shí)際中，常常要處理由實(shí)驗(yàn)或測量所得到的一些離散數(shù)據(jù)。插值與擬合方法就是要通過這些數(shù)據(jù)去確定某一類已知函數(shù)的參數(shù)或?qū)で竽硞€(gè)近似函數(shù)，使所得到的近似函數(shù)與已知數(shù)據(jù)有較高的擬合精度。如果要求這個(gè)近似函數(shù)（曲線或曲面）經(jīng)過所已知的所有數(shù)據(jù)點(diǎn)，則稱此類問題為插值問題。（不需要函數(shù)表達(dá)式）二、基本概念第5頁/共107頁第五頁，共108頁。如果不要求近似函數(shù)通過所有數(shù)據(jù)點(diǎn)，而是要求它能較好地反映數(shù)據(jù)變化規(guī)律的近似函數(shù)的方法稱為數(shù)據(jù)擬合。（必須有函數(shù)表達(dá)式）近似函數(shù)不一定（曲線或曲面）通過所有的數(shù)據(jù)點(diǎn)。第6頁/共107頁第六頁，共108頁。1、聯(lián)系都是根據(jù)實(shí)際中一組已知數(shù)據(jù)來構(gòu)造一個(gè)能夠反映數(shù)據(jù)變化規(guī)律的近似函數(shù)的方法。2、區(qū)別插值問題不一定得到近似函數(shù)的表達(dá)形式，僅通過插值方法找到未知點(diǎn)對應(yīng)的值。數(shù)據(jù)擬合要求得到一個(gè)具體的近似函數(shù)的表達(dá)式。

三、插值與擬合的區(qū)別和聯(lián)系第7頁/共107頁第七頁，共108頁。

插值當(dāng)數(shù)據(jù)量不夠，需要補(bǔ)充，且認(rèn)定已有數(shù)據(jù)可信時(shí),通常利用函數(shù)插值方法。實(shí)際問題當(dāng)中碰到的函數(shù)f(x)是各種各樣的，有的表達(dá)式很復(fù)雜，有的甚至給不出數(shù)學(xué)的式子，只提供了一些離散數(shù)據(jù)，警如，某些點(diǎn)上的函數(shù)值和導(dǎo)數(shù)值。第8頁/共107頁第八頁，共108頁。拉格朗日插值分段線性插值三次樣條插值一維插值一、插值的定義二、插值的方法三、用Matlab解插值問題返回第9頁/共107頁第九頁，共108頁。返回二維插值一、二維插值定義二、網(wǎng)格節(jié)點(diǎn)插值法三、用Matlab解插值問題最鄰近插值分片線性插值雙線性插值網(wǎng)格節(jié)點(diǎn)數(shù)據(jù)的插值散點(diǎn)數(shù)據(jù)的插值第10頁/共107頁第十頁，共108頁。一維插值的定義已知n+1個(gè)節(jié)點(diǎn)其中互不相同，不妨設(shè)求任一插值點(diǎn)處的插值節(jié)點(diǎn)可視為由產(chǎn)生,，表達(dá)式復(fù)雜,，或無封閉形式,，或未知.。第11頁/共107頁第十一頁，共108頁。

構(gòu)造一個(gè)(相對簡單的)函數(shù)通過全部節(jié)點(diǎn),即再用計(jì)算插值，即返回第12頁/共107頁第十二頁，共108頁。

稱為拉格朗日插值基函數(shù)。

已知函數(shù)f(x)在n+1個(gè)點(diǎn)x0,x1,…,xn處的函數(shù)值為y0,y1,…,yn

。求一n次多項(xiàng)式函數(shù)Pn(x)，使其滿足：

Pn(xi)=yi,i=0,1,…,n.

解決此問題的拉格朗日插值多項(xiàng)式公式如下其中Li(x)為n次多項(xiàng)式：拉格朗日(Lagrange)插值第13頁/共107頁第十三頁，共108頁。拉格朗日(Lagrange)插值特別地:兩點(diǎn)一次(線性)插值多項(xiàng)式:三點(diǎn)二次(拋物)插值多項(xiàng)式:第14頁/共107頁第十四頁，共108頁。

拉格朗日多項(xiàng)式插值的這種振蕩現(xiàn)象叫Runge現(xiàn)象

采用拉格朗日多項(xiàng)式插值：選取不同插值節(jié)點(diǎn)個(gè)數(shù)n+1，其中n為插值多項(xiàng)式的次數(shù)，當(dāng)n分別取2,4,6,8,10時(shí)，繪出插值結(jié)果圖形.例返回ToMatlablch(larg1)第15頁/共107頁第十五頁，共108頁。第16頁/共107頁第十六頁，共108頁。分段線性插值n越大，誤差越小.xjxj-1xj+1x0xnxoy第17頁/共107頁第十七頁，共108頁。ToMATLABxch11，xch12，xch13，xch14返回例用分段線性插值法求插值,并觀察插值誤差.1.在[-6,6]中平均選取5個(gè)點(diǎn)作插值(xch11)4.在[-6,6]中平均選取41個(gè)點(diǎn)作插值(xch14)2.在[-6,6]中平均選取11個(gè)點(diǎn)作插值(xch12)3.在[-6,6]中平均選取21個(gè)點(diǎn)作插值(xch13)第18頁/共107頁第十八頁，共108頁。第19頁/共107頁第十九頁，共108頁。比分段線性插值更光滑。xyxi-1xiab

在數(shù)學(xué)上，光滑程度的定量描述是：函數(shù)(曲線)的k階導(dǎo)數(shù)存在且連續(xù)，則稱該曲線具有k階光滑性。光滑性的階次越高，則越光滑。是否存在較低次的分段多項(xiàng)式達(dá)到較高階光滑性的方法？三次樣條插值就是一個(gè)很好的例子。三次樣條插值第20頁/共107頁第二十頁，共108頁。

三次樣條插值g(x)為被插值函數(shù)。第21頁/共107頁第二十一頁，共108頁。例用三次樣條插值選取11個(gè)基點(diǎn)計(jì)算插值(ych)返回ToMATLABych第22頁/共107頁第二十二頁，共108頁。用MATLAB作插值計(jì)算一維插值函數(shù)：yi=interp1(x，y，xi，'method')插值方法被插值點(diǎn)插值節(jié)點(diǎn)xi處的插值結(jié)果‘nearest’

：最鄰近插值‘linear’

：線性插值；‘spline’

：三次樣條插值；‘cubic’

：立方插值。缺省時(shí)：分段線性插值。注意：所有的插值方法都要求x是單調(diào)的，并且xi不能夠超過x的范圍。第23頁/共107頁第二十三頁，共108頁。（1）nearest方法速度最快，占用內(nèi)存最小，但一般來說誤差最大，插值結(jié)果最不光滑；(2)‘liner’分段線性插值：插值點(diǎn)處函數(shù)值由連接其最鄰近的兩側(cè)點(diǎn)的線性函數(shù)預(yù)測，MATLAB中interp1的默認(rèn)方法（3）spline三次樣條插值是所有插值方法中運(yùn)行耗時(shí)最長的，其插值函數(shù)以及插值函數(shù)的一階、二階導(dǎo)函數(shù)都連續(xù)，因此是最光滑的插值方法，占用內(nèi)存上比cubic方法小，但當(dāng)已知數(shù)據(jù)點(diǎn)不均勻分布時(shí)可能出現(xiàn)異常結(jié)果。（4）cubic三次多項(xiàng)式插值法中插值函數(shù)及其一階導(dǎo)數(shù)都是連續(xù)的，因此其插值結(jié)果也比較光滑，運(yùn)算速度比spline方法略快，但占用內(nèi)存最多。在實(shí)際的使用中，應(yīng)根據(jù)實(shí)際需求和運(yùn)算條件選擇合適的算法。第24頁/共107頁第二十四頁，共108頁。例用其他一維插值方法對以下7個(gè)離散數(shù)據(jù)點(diǎn)

(1,3.5)、(2,2.1)、(3,1.3)、(4.0.8)、(5,2.9)、(6,4.2)、(7,5.7)進(jìn)行一維插值方法。解：在MATLAB命令窗口中輸入以下命令：

>>x=[1234567];>>y=[3.52.11.30.82.94.25.7];>>xx=1:0.5:7;>>y1=interp1(x,y,xx,'nearest');>>y2=interp1(x,y,xx,'spline');>>y3=interp1(x,y,xx,'cubic');>>plot(x,y,'o',xx,y1,'-',xx,y2,'-.',xx,y3,':')第25頁/共107頁第二十五頁，共108頁。第26頁/共107頁第二十六頁，共108頁。

例：在1-12的11小時(shí)內(nèi)，每隔1小時(shí)測量一次溫度，測得的溫度依次為：5，8，9，15，25，29，31，30，22，25，27，24。試估計(jì)每隔1/10小時(shí)的溫度值。ToMATLAB(temp)hours=1:12;temps=[589152529313022252724];h=1:0.1:12;t=interp1(hours,temps,h,'spline');(直接輸出數(shù)據(jù)將是很多的)plot(hours,temps,'+‘)holdonplot(h,t,hours,temps,'r:')%作圖xlabel('Hour'),ylabel('DegreesCelsius’)第27頁/共107頁第二十七頁，共108頁。第28頁/共107頁第二十八頁，共108頁。xy機(jī)翼下輪廓線例已知飛機(jī)下輪廓線上數(shù)據(jù)如下，求x每改變0.1時(shí)的y值。ToMATLAB(plane)返回第29頁/共107頁第二十九頁，共108頁。第30頁/共107頁第三十頁，共108頁。二維插值的定義xyO第一種（網(wǎng)格節(jié)點(diǎn)）：第31頁/共107頁第三十一頁，共108頁。

已知mn個(gè)節(jié)點(diǎn)其中互不相同，不妨設(shè)

構(gòu)造一個(gè)二元函數(shù)通過全部已知節(jié)點(diǎn),即再用計(jì)算插值，即第32頁/共107頁第三十二頁，共108頁。第二種（散亂節(jié)點(diǎn)）：yx0第33頁/共107頁第三十三頁，共108頁。已知n個(gè)節(jié)點(diǎn)其中互不相同，

構(gòu)造一個(gè)二元函數(shù)通過全部已知節(jié)點(diǎn),即再用計(jì)算插值，即返回第34頁/共107頁第三十四頁，共108頁。

注意：最鄰近插值一般不連續(xù)。具有連續(xù)性的最簡單的插值是分片線性插值。最鄰近插值xy(x1,y1)(x1,y2)(x2,y1)(x2,y2)O

二維或高維情形的最鄰近插值，與被插值點(diǎn)最鄰近的節(jié)點(diǎn)的函數(shù)值即為所求。返回第35頁/共107頁第三十五頁，共108頁。

將四個(gè)插值點(diǎn)（矩形的四個(gè)頂點(diǎn)）處的函數(shù)值依次簡記為：

分片線性插值xy(xi,yj)(xi,yj+1)(xi+1,yj)(xi+1,yj+1)Of(xi,yj)=f1，f(xi+1,yj)=f2，f(xi+1,yj+1)=f3，f(xi,yj+1)=f4第36頁/共107頁第三十六頁，共108頁。插值函數(shù)為：第二片(上三角形區(qū)域)：(x,y)滿足插值函數(shù)為：注意：(x,y)當(dāng)然應(yīng)該是在插值節(jié)點(diǎn)所形成的矩形區(qū)域內(nèi)。顯然，分片線性插值函數(shù)是連續(xù)的；分兩片的函數(shù)表達(dá)式如下：第一片(下三角形區(qū)域)：(x,y)滿足返回第37頁/共107頁第三十七頁，共108頁。

雙線性插值是一片一片的空間二次曲面構(gòu)成。雙線性插值函數(shù)的形式如下：其中有四個(gè)待定系數(shù)，利用該函數(shù)在矩形的四個(gè)頂點(diǎn)（插值節(jié)點(diǎn)）的函數(shù)值，得到四個(gè)代數(shù)方程，正好確定四個(gè)系數(shù)。雙線性插值xy(x1,y1)(x1,y2)(x2,y1)(x2,y2)O返回第38頁/共107頁第三十八頁，共108頁。

要求x0,y0單調(diào)；x，y可取為矩陣，或x取行向量，y取為列向量，x,y的值分別不能超出x0,y0的范圍。z=interp2(x0,y0,z0,x,y,’method’)被插值點(diǎn)插值方法用MATLAB作網(wǎng)格節(jié)點(diǎn)數(shù)據(jù)的插值插值節(jié)點(diǎn)被插值點(diǎn)的函數(shù)值‘nearest’

最鄰近插值‘linear’

雙線性插值‘cubic’

雙三次插值缺省時(shí),雙線性插值第39頁/共107頁第三十九頁，共108頁。例：測得平板表面3*5網(wǎng)格點(diǎn)處的溫度分別為：828180828479636165818484828586試作出平板表面的溫度分布曲面z=f(x,y)的圖形。輸入以下命令：x=1:5;y=1:3;temps=[8281808284;7963616581;8484828586];mesh(x,y,temps)1.先在三維坐標(biāo)畫出原始數(shù)據(jù)，畫出粗糙的溫度分布曲圖.2．以平滑數(shù)據(jù),在x、y方向上每隔0.2個(gè)單位的地方進(jìn)行插值.第40頁/共107頁第四十頁，共108頁。第41頁/共107頁第四十一頁，共108頁。再輸入以下命令:xi=1:0.2:5;yi=1:0.2:3;zi=interp2(x,y,temps,xi',yi,'cubic');mesh(xi,yi,zi)畫出插值后的溫度分布曲面圖.ToMATLAB(wendu)第42頁/共107頁第四十二頁，共108頁。第43頁/共107頁第四十三頁，共108頁。例>>[x,y]=meshgrid(-3:.6:3,-2:.4:2);>>z=(x.^2-2*x).*exp…(-x.^2-y.^2-x.*y);>>surf(x,y,z),第44頁/共107頁第四十四頁，共108頁。選較密的插值點(diǎn)，用默認(rèn)的線性插值算法進(jìn)行插值>>[x1,y1]=meshgrid(-3:.2:3,-2:.2:2);>>z0=interp2(x,y,z,x1,y1);>>surf(x1,y1,z0)第45頁/共107頁第四十五頁，共108頁。立方和樣條插值：>>z1=interp2(x,y,z,x1,y1,'cubic');>>z2=interp2(x,y,z,x1,y1,'spline');>>surf(x1,y1,z1),axis([-3,3,-2,2,-0.7,1.5])>>figure;surf(x1,y1,z2),axis([-3,3,-2,2,-0.7,1.5])第46頁/共107頁第四十六頁，共108頁。算法誤差比較>>z=(x1.^2-2*x1).*exp(-x1.^2-y1.^2-x1.*y1);>>surf(x1,y1,abs(z-z1))>>figure;surf(x1,y1,abs(z-z2))>>figure;surf(x1,y1,abs(z-z0))第47頁/共107頁第四十七頁，共108頁。山區(qū)地形地貌圖已知某處山區(qū)地形選點(diǎn)測量坐標(biāo)數(shù)據(jù)為：x=0

0.5

1

1.5

2

2.5

3

3.5

4

4.5

5y=0

0.5

1

1.5

2

2.5

3

3.5

4

4.5

5

5.5

6海拔高度數(shù)據(jù)為：z=8990878592919693908782

9296989995918986848284

9698959290888584838185

8081828995969392898686

8285879899969788858283

8285899495939291868488

8892939495898786838192

9296979896939584828184

8585818280808185909395

8486819899989796958487

8081858283848790958688

8082818485868382818082

8788899899979698949287第48頁/共107頁第四十八頁，共108頁。山區(qū)地形地貌圖程序原始地貌圖程序：x=0:.5:5;y=0:.5:6;[xx,yy]=meshgrid(x,y);z=[8990878592919693908782929698999591898684828496989592908885848381858081828995969392898686828587989996978885828382858994959392918684888892939495898786838192929697989693958482818485858182808081859093958486819899989796958487808185828384879095868880828184858683828180828788899899979698949287];mesh(xx,yy,z)加密后的地貌圖x=0:.5:5;y=0:.5:6;z=[8990878592919693908782929698999591898684828496989592908885848381858081828995969392898686828587989996978885828382858994959392918684888892939495898786838192929697989693958482818485858182808081859093958486819899989796958487808185828384879095868880828184858683828180828788899899979698949287];xi=linspace(0,5,50);%加密橫坐標(biāo)數(shù)據(jù)到50個(gè)yi=linspace(0,6,80);%加密縱坐標(biāo)數(shù)據(jù)到60個(gè)[xii,yii]=meshgrid(xi,yi);%生成網(wǎng)格數(shù)據(jù)zii=interp2(x,y,z,xii,yii,'cubic');%插值mesh(xii,yii,zii)%加密后的地貌圖第49頁/共107頁第四十九頁，共108頁。山區(qū)地形地貌圖結(jié)果第50頁/共107頁第五十頁，共108頁。

通過此例對最近鄰點(diǎn)插值、雙線性插值方法和雙三次插值方法的插值效果進(jìn)行比較。ToMATLAB(moutain)返回第51頁/共107頁第五十一頁，共108頁。第52頁/共107頁第五十二頁，共108頁。第53頁/共107頁第五十三頁，共108頁。第54頁/共107頁第五十四頁，共108頁。第55頁/共107頁第五十五頁，共108頁。第56頁/共107頁第五十六頁，共108頁。

插值函數(shù)griddata格式為:

cz

=griddata（x，y，z，cx，cy，‘method’）用MATLAB作散點(diǎn)數(shù)據(jù)的插值計(jì)算

要求cx取行向量，cy取為列向量。被插值點(diǎn)插值方法插值節(jié)點(diǎn)被插值點(diǎn)的函數(shù)值‘nearest’

最鄰近插值‘linear’

雙線性插值‘cubic’

雙三次插值'v4'-Matlab提供的插值方法缺省時(shí),雙線性插值第57頁/共107頁第五十七頁，共108頁。>>x=-3+6*rand(200,1);y=-2+4*rand(200,1);>>z=(x.^2-2*x).*exp(-x.^2-y.^2-x.*y);>>[x1,y1]=meshgrid(-3:.2:3,-2:.2:2);>>z1=griddata(x,y,z,x1,y1,'cubic');>>surf(x1,y1,z1),axis([-3,3,-2,2,-0.7,1.5])>>z2=griddata(x,y,z,x1,y1,'v4');>>figure;surf(x1,y1,z2),axis([-3,3,-2,2,-0.7,1.5])例：在x為［3，3］，y為［－2，2］矩形區(qū)域隨機(jī)選擇一組坐標(biāo)，用’v4’與’cubic’插值法進(jìn)行處理，并對誤差進(jìn)行比較。第58頁/共107頁第五十八頁，共108頁。第59頁/共107頁第五十九頁，共108頁。誤差分析>>z0=(x1.^2-2*x1).*exp(-x1.^2-y1.^2-x1.*y1);>>surf(x1,y1,abs(z0-z1)),axis([-3,3,-2,2,0,0.15])>>figure;surf(x1,y1,abs(z0-z2)),axis([-3,3,-2,2,0,0.15])第60頁/共107頁第六十頁，共108頁。第61頁/共107頁第六十一頁，共108頁。例：

在x為［3，3］，y為［－2，2］矩形區(qū)域隨機(jī)選擇一組坐標(biāo)中，對分布不均勻數(shù)據(jù)，進(jìn)行插值分析。>>x=-3+6*rand(200,1);y=-2+4*rand(200,1);>>z=(x.^2-2*x).*exp(-x.^2-y.^2-x.*y);%生成已知數(shù)據(jù)>>plot(x,y,'x')%樣本點(diǎn)的二維分布>>figure,plot3(x,y,z,'x'),axis([-3,3,-2,2,-0.7,1.5])，gridon第62頁/共107頁第六十二頁，共108頁。去除在(-1,-1/2)點(diǎn)為圓心，以0.5為半徑的圓內(nèi)的點(diǎn)。>>x=-3+6*rand(200,1);y=-2+4*rand(200,1);%重新生成樣本點(diǎn)>>z=(x.^2-2*x).*exp(-x.^2-y.^2-x.*y);>>ii=find((x+1).^2+(y+0.5).^2>0.5^2);%找出不在圓內(nèi)的點(diǎn)坐標(biāo)>>x=x(ii);y=y(ii);z=z(ii);plot(x,y,'x')>>t=[0:.1:2*pi,2*pi];x0=-1+0.5*cos(t);y0=-0.5+0.5*sin(t);>>line(x0,y0)第63頁/共107頁第六十三頁，共108頁。新樣本點(diǎn)擬合曲面>>[x1,y1]=meshgrid(-3:.2:3,-2:.2:2);>>z1=griddata(x,y,z,x1,y1,'v4');>>surf(x1,y1,z1),axis([-3,3,-2,2,-0.7,1.5])第64頁/共107頁第六十四頁，共108頁。誤差分析>>z0=(x1.^2-2*x1).*exp(-x1.^2-y1.^2-x1.*y1);>>surf(x1,y1,abs(z0-z1)),axis([-3,3,-2,2,0,0.1])>>contour(x1,y1,abs(z0-z1),30);holdon,plot(x,y,‘x’);line(x0,y0)％誤差的二維等高線圖第65頁/共107頁第六十五頁，共108頁。命令3

interp3

三維網(wǎng)格生成用meshgrid()函數(shù),調(diào)用格式:[x,y,z]=meshgrid(x1,y1,z1)

其中x1,y1,z1為這三維所需要的分割形式，應(yīng)以向量形式給出，返回x,y,z為網(wǎng)格的數(shù)據(jù)生成，均為三維數(shù)組。

griddata3()三維非網(wǎng)格形式的插值擬合命令4

interpnn維網(wǎng)格生成用ndgrid()函數(shù),調(diào)用格式:[x1,x2,…,xn]=ndgrid[v1,v2,…,vn]griddatan()n維非網(wǎng)格形式的插值擬合interp3()、interpn()調(diào)用格式同interp2()函數(shù)一致；griddata3()、griddatan()調(diào)用格式同griddata()函數(shù)一致。第66頁/共107頁第六十六頁，共108頁。例：通過函數(shù)生成一些網(wǎng)格型樣本點(diǎn)，試根據(jù)樣本點(diǎn)進(jìn)行擬合，并給出擬合誤差。>>[x,y,z]=meshgrid(-1:0.2:1);[x0,y0,z0]=meshgrid(-1:0.05:1);>>V=exp(x.^2.*z+y.^2.*x+z.^2.*y…).*cos(x.^2.*y.*z+z.^2.*y.*x);>>V0=exp(x0.^2.*z0+y0.^2.*x0…+z0.^2.*y0).*cos(x0.^2.*y0.*z0+z0.^2.*y0.*x0);>>V1=interp3(x,y,z,V,x0,y0,z0,'spline');err=V1-V0;max(err(:))ans=0.0419第67頁/共107頁第六十七頁，共108頁。slice(x0,y0,z0,V1,[-0.5,0.3,0.9],[0.6,-0.1],[-1,-0.5,0.5,1])title('SlivesforFourDimFigures')第68頁/共107頁第六十八頁，共108頁。[x,y,z]=meshgrid(-3:.1:3,-4:.2:4,-4:.4:4);v=x.*exp(-x.^2-y.^2-z.^2);xslice=[-2.2,.5,3];yslice=[3,0.9];zslice=[-4,1];slice(x,y,z,v,xslice,yslice,zslice);第69頁/共107頁第六十九頁，共108頁。

例在某海域測得一些點(diǎn)(x,y)處的水深z由下表給出，船的吃水深度為5英尺，在矩形區(qū)域（75，200）*（-50，150）里的哪些地方船要避免進(jìn)入。第70頁/共107頁第七十頁，共108頁。ToMATLABhd1返回4.作出水深小于5的海域范圍,即z=5的等高線.第71頁/共107頁第七十一頁，共108頁。第72頁/共107頁第七十二頁，共108頁。第73頁/共107頁第七十三頁，共108頁。%程序一：插值并作海底曲面圖

x=[129.0140.0103.588.0185.5195.0105.5157.5107.577.081.0162.0162.0117.5];y=[7.5141.523.0147.022.5137.585.5-6.5-813.056.5-66.584.0-33.5];z=[48686889988949];x1=75:1:200;y1=-50:1:150;[x1,y1]=meshgrid(x1,y1);z1=griddata(x,y,z,x1,y1,'v4');meshc(x1,y1,z1)第74頁/共107頁第七十四頁，共108頁。海底曲面圖第75頁/共107頁第七十五頁，共108頁。%程序二：插值并作出水深小于5的海域范圍。x1=75:1:200;y1=-50:1:150;[x1,y1]=meshgrid(x1,y1);z1=griddata(x,y,z,x1,y1,'v4');

%插值z1(z1>=5)=nan;

%將水深大于5的置為nan，這樣繪圖就不會顯示出來meshc(x1,y1,z1)

第76頁/共107頁第七十六頁，共108頁。水深小于5的海域范圍第77頁/共107頁第七十七頁，共108頁。拉格朗日插值法拉格朗日插值法是基于基函數(shù)的插值方法，插值多項(xiàng)式可表示為其中稱為i次基函數(shù)：第78頁/共107頁第七十八頁，共108頁。在MATLAB中編程實(shí)現(xiàn)拉格朗日插值法函數(shù)為：Language。功能：求已知數(shù)據(jù)點(diǎn)的拉格朗日多項(xiàng)式；調(diào)用格式：f=

Language(x,y)或f=

Language(x,y,x0)。其中，x為已知數(shù)據(jù)點(diǎn)的x坐標(biāo)向量；

y為已知數(shù)據(jù)點(diǎn)的y坐標(biāo)向量；

x0為插值點(diǎn)的x坐標(biāo)；

f為求得的拉格朗日多項(xiàng)式或x0處的插值。第79頁/共107頁第七十九頁，共108頁。

functionf=Language(x,y,x0)%求已知數(shù)據(jù)點(diǎn)的拉格朗日多項(xiàng)式%已知數(shù)據(jù)點(diǎn)的x坐標(biāo)向量:x%已知數(shù)據(jù)點(diǎn)的y坐標(biāo)向量:y%為插值點(diǎn)的x坐標(biāo):x0%求得的拉格朗日多項(xiàng)式或x0處的插值:fsymst;if(length(x)==length(y))n=length(x);elsedisp('x和y的維數(shù)不相等！');return;end%檢錯(cuò)f=0.0;for(i=1:n)l=y(i);for(j=1:i-1)l=l*(t-x(j))/(x(i)-x(j));end;for(j=i+1:n)l=l*(t-x(j))/(x(i)-x(j));%計(jì)算拉格朗日基函數(shù)

end;f=f+l;%計(jì)算拉格朗日插值函數(shù)

simplify(f);%化簡

if(i==n)if(nargin==3)f=subs(f,'t',x0);%計(jì)算插值點(diǎn)的函數(shù)值

elsef=collect(f);%將插值多項(xiàng)式展開

f=vpa(f,6);%將插值多項(xiàng)式的系數(shù)化成6位精度的小數(shù)

endendend第80頁/共107頁第八十頁，共108頁。例4-3根據(jù)下表的數(shù)據(jù)點(diǎn)求出其拉格朗日插值多項(xiàng)式，并計(jì)算當(dāng)x=1.6時(shí)y的值。解：>>

x=[11.21.82.54];>>y=[0.84150.93200.97380.5985-0.7568];>>f=language(x,y)f=1.05427*t-.145485e-1*t^2-.204917*t^3+.328112e-1*t^4-.261189e-1>>f=language(x,y,1.6)f=0.9992x11.21.82.54y0.84150.93200.97380.5985-0.7568第81頁/共107頁第八十一頁，共108頁。利用均差的牛頓插值法函數(shù)f的零階均差定義為，一階定義均差為一般地，函數(shù)f的k階均差定義為：利用均差的牛頓插值法多項(xiàng)式為

第82頁/共107頁第八十二頁，共108頁。系數(shù)的計(jì)算過程如表所示表

均差計(jì)算表格一階均差二階均差三階均差……n階均差……………第83頁/共107頁第八十三頁，共108頁。在MATLAB中編程實(shí)現(xiàn)利用均差牛頓插值法函數(shù)為：Newton。功能：求已知數(shù)據(jù)點(diǎn)的均差形式的牛頓插值多項(xiàng)式；調(diào)用格式：f=

Newton(x,y)或f=

Newton(x,y,x0)。其中，x為已知數(shù)據(jù)點(diǎn)的x坐標(biāo)向量；

y為已知數(shù)據(jù)點(diǎn)的y坐標(biāo)向量；

x0為插值點(diǎn)的x坐標(biāo)；

f為求得的牛頓插值法多項(xiàng)式或x0處的插值。第84頁/共107頁第八十四頁，共108頁。functionf=Newton(x,y,x0)%求已知數(shù)據(jù)點(diǎn)的均差形式牛頓插值多項(xiàng)式%已知數(shù)據(jù)點(diǎn)的x坐標(biāo)向量:x%已知數(shù)據(jù)點(diǎn)的y坐標(biāo)向量:y%為插值點(diǎn)的x坐標(biāo):x0%求得的均差形式牛頓插值多項(xiàng)式或x0處的插值:fsymst;if(length(x)==length(y))n=length(x);c(1:n)=0.0;elsedisp('x和y的維數(shù)不相等！');return;end第85頁/共107頁第八十五頁，共108頁。f=y(1);y1=0;l=1;

for(i=1:n-1)for(j=i+1:n)y1(j)=(y(j)-y(i))/(x(j)-x(i));endc(i)=y1(i+1);l=l*(t-x(i));f=f+c(i)*l;simplify(f);y=y1;

if(i==n-1)if(nargin==3)f=subs(f,'t',x0);elsef=collect(f);%將插值多項(xiàng)式展開

f=vpa(f,6);endendend第86頁/共107頁第八十六頁，共108頁。例根據(jù)下表的數(shù)據(jù)點(diǎn)求出其均差形式牛頓插值多項(xiàng)式，并計(jì)算當(dāng)x=2.0時(shí)y的值。解：>>x=[11.21.82.54];>>y=[11.443.246.2516];>>f=Newton(x,y)f=.182711e-14-.482154e-14*t+1.00000*t^2-.169177e-14*t^3+.211471e-15*t^4>>f=Newton(x,y,2.0)f=4x11.21.82.54y11.443.246.2516第87頁/共107頁第八十七頁，共108頁。利用差分的牛頓插值法差分分為向前差分、后向差分和中心差分三種，它們的記法及定義如下為：

其中：代表向前差分；代表向后差分；代表向后差分。第88頁/共107頁第八十八頁，共108頁。

假設(shè)。為了方便，可構(gòu)造如表所示的差分表（）。表

差分計(jì)算表格………………第89頁/共107頁第八十九頁，共108頁。向前牛頓插值向前牛頓插值多項(xiàng)式可表示如下：其中叫做步長，，且的取值范圍為。第90頁/共107頁第九十頁，共108頁。在MATLAB中編程實(shí)現(xiàn)向前牛頓插值法函數(shù)為：Newtonforward。功能：求已知數(shù)據(jù)點(diǎn)的向前牛頓插值法多項(xiàng)式；調(diào)用格式：f=

Newtonforward(x,y)或

f=

Newtonforward

(x,y,x0)。其中，x為已知數(shù)據(jù)點(diǎn)的x坐標(biāo)向量；

y為已知數(shù)據(jù)點(diǎn)的y坐標(biāo)向量；

x0為插值點(diǎn)的x坐標(biāo)；

f為求得的向前牛頓插值法多項(xiàng)式或x0處的插值。第91頁/共107頁第九十一頁，共108頁。functionf=Newtonforward(x,y,x0)%求已知數(shù)據(jù)點(diǎn)的向前差分牛頓插值多項(xiàng)式%已知數(shù)據(jù)點(diǎn)的x坐標(biāo)向量:x%已知數(shù)據(jù)點(diǎn)的y坐標(biāo)向量:y%為插值點(diǎn)的x坐標(biāo):x0%求得的向前差分牛頓插值多項(xiàng)式或x0處的插值:fsymst;if(length(x)==length(y))n=length(x);c(1:n)=0.0;elsedisp('x和y的維數(shù)不相等！');return;end第92頁/共107頁第九十二頁，共108頁。f=y(1);y1=0;xx=linspace(x(1),x(n),(x(2)-x(1)));if(xx~=x)disp('節(jié)點(diǎn)之間不是等距的！');return;endfor(i=1:n-1)for(j=1:n-i)y1(j)=y(j+1)-y(j);endc(i)=y1(1);l=t;for(k=1:i-1)l=l*(t-k);end;f=f+c(i)*l/factorial(i);simplify(f);y=y1;if(i==n-1)if(nargin==3)f=subs(f,'t',(x0-x(1))/(x(2)-x(1)));elsef=collect(f);f=vpa(f,6);endendend第93頁/共107頁第九十三頁，共108頁。向前牛頓插值向后牛頓插值多項(xiàng)式可表示如下：其中叫做步長，，且的取值范圍為。第94頁/共107頁第九十四頁，共108頁。在MATLAB中編程實(shí)現(xiàn)向后牛頓插值法函數(shù)為：Newtonback。功能：求已知數(shù)據(jù)點(diǎn)的向前牛頓插值法多項(xiàng)式；調(diào)用格式：f=

Newtonback

(x,y)或

f=

Newtonback(x,y,x0)。其中，x為已知數(shù)據(jù)點(diǎn)的x坐標(biāo)向量；

y為已知數(shù)據(jù)點(diǎn)的y坐標(biāo)向量；

x0為插值點(diǎn)的x坐標(biāo)；

f為求得的向前牛頓插值法多項(xiàng)式或x0處的插值。第95頁/共107頁第九十五頁，共108頁。functionf=Newtonback(x,y,x0)%求已知數(shù)據(jù)點(diǎn)的向后差分牛頓插值多項(xiàng)式%已知數(shù)據(jù)點(diǎn)的x坐標(biāo)向量:x%已知數(shù)據(jù)點(diǎn)的y坐標(biāo)向量:y%為插值點(diǎn)的x坐標(biāo):x0%求得的向前差分牛頓插值多項(xiàng)式或x0處的插值:fsymst;if(length(x)==length(y))n=length(x);c(1:n)=0.0;elsedisp('x和y的維數(shù)不相等！');return;end第96頁/共107頁第九十六頁，共108頁。f=y(n);y1=0;xx=linspace(x(1),x(n),(x(2)-x(1)));if(xx~=x)disp('節(jié)點(diǎn)之間不是等距的！');return;endfor(i=1:n-1)for(j=i+1:n)y1(j)=y(j)-y(j-1);endc(i)=y1(n);l=t;for(k=1:i-1)l=l*(t+k);end;

f=f+c(i)*l/factorial(i);simplify(f);y=y1;if(i==n-1)if(nargin==3)f=subs(f,'t',(x0-x(n))/(x(2)-x(1)));elsef=collect(f);f=vpa(f,6);endendend第97頁/共107頁第九十七頁，共108頁。例5根據(jù)下表的數(shù)據(jù)點(diǎn)求出其差分形式的牛頓插值多項(xiàng)式，并計(jì)算當(dāng)x=1.55時(shí)y的值。解：>>x=[11.21.41.61.8];>>y=[0.84150.9320 0.98540.99960.9738];>>f=Newtonforward(x,y)f=.841500+.108025*t-.169042e-1*t^2-.675000e-3*t^3+.541667e-4*t^4>>f=Newtonforward(x,y,1.55)f=0.9998f=Newtonback(x,y)f=.973800-.457417e-1*t-.198042e-1*t^2+.191667e-3*t^3+.541667e-4*t^4>>f=Newtonback(x,y,1.55)f=0.9998x11.21.41.61.8y0.84150.93200.98540.99960.9738第98頁/共107頁第九十八頁，共108頁。Hermite插值Hermite插值滿足在節(jié)點(diǎn)上等于給定函數(shù)值，而且在節(jié)點(diǎn)上的導(dǎo)數(shù)值也等于給定的導(dǎo)數(shù)值。對于高階導(dǎo)數(shù)的情況，Hermite插值多項(xiàng)式比較復(fù)雜，在實(shí)際中，常常遇到的是函數(shù)值與一階導(dǎo)數(shù)給定的情況。在此情況下，n個(gè)節(jié)點(diǎn)x1，x2，…，xn的Hermite插值多項(xiàng)式的表達(dá)式如下：其中，，，第99頁/共107頁第九十九頁，共108頁。在MATLAB中編程實(shí)現(xiàn)Hermite插值法函數(shù)為：Hermite。功能：求已知數(shù)據(jù)點(diǎn)的Hermite插值法多項(xiàng)式；調(diào)用格式：f=

Hermite

(x,y,y_1)或

f=

Hermite

(x,y,y_1,x0)。其中，x為已知數(shù)據(jù)點(diǎn)的x坐標(biāo)向量；

y為已知數(shù)據(jù)點(diǎn)的y坐標(biāo)向量；

y_1為已知數(shù)據(jù)點(diǎn)導(dǎo)數(shù)向量；

x0為插值點(diǎn)的x坐標(biāo)；

f為求得的Hermite插值法多項(xiàng)式或x0處的插值。第100頁/共107頁第一百頁，共108頁。functionf=Hermite(x,y,y_1,x0)%求已知數(shù)據(jù)點(diǎn)的向后差分牛頓插值多項(xiàng)式%已知數(shù)據(jù)點(diǎn)的x坐標(biāo)向量:x%已知數(shù)據(jù)點(diǎn)的y坐標(biāo)向量:y%已知數(shù)據(jù)點(diǎn)的導(dǎo)數(shù)向量:y_1%求得的Hermite插值多項(xiàng)式或x0處的插值:fsymst;f=0.0;if(length(x)==length(y))if(length(y)==length(y_1))n=length(x);elsedisp('y和y的導(dǎo)數(shù)的維數(shù)不相等！');return;endelsedisp('x和y的維數(shù)不相等！');return;end第101頁/共107頁第一百零一頁，共108頁。fori=1:nh=1.0;a=0.0;forj=1:nif(j~=i)h=h*(t-x(j))^2/((x(i)-x(j))^2);a=a+1/(x(i)-x(j));endend

f=f+h*((x(i)-t)*(2*a*y(i)-y_1(i))+y(i));

if(i==n)if(nargin==4)f=subs(f,'t',x0);elsef=vpa(f,6);endendend第102頁/共107頁第一百零二頁，共108頁。例6根據(jù)