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

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

三、插值與擬合的區(qū)別和聯(lián)系第8頁，課件共107頁，創(chuàng)作于2023年2月

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

構(gòu)造一個(相對簡單的)函數(shù)通過全部節(jié)點,即再用計算插值，即返回第13頁，課件共107頁，創(chuàng)作于2023年2月

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

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

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

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

解決此問題的拉格朗日插值多項式公式如下其中Li(x)為n次多項式：拉格朗日(Lagrange)插值第14頁，課件共107頁，創(chuàng)作于2023年2月拉格朗日(Lagrange)插值特別地:兩點一次(線性)插值多項式:三點二次(拋物)插值多項式:第15頁，課件共107頁，創(chuàng)作于2023年2月

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

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

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

三次樣條插值g(x)為被插值函數(shù)。第22頁，課件共107頁，創(chuàng)作于2023年2月例用三次樣條插值選取11個基點計算插值(ych)返回ToMATLABych第23頁，課件共107頁，創(chuàng)作于2023年2月用MATLAB作插值計算一維插值函數(shù)：yi=interp1(x，y，xi，'method')插值方法被插值點插值節(jié)點xi處的插值結(jié)果‘nearest’

：最鄰近插值‘linear’

：線性插值；‘spline’

：三次樣條插值；‘cubic’

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

(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,':')第26頁，課件共107頁，創(chuàng)作于2023年2月第27頁，課件共107頁，創(chuàng)作于2023年2月

例：在1-12的11小時內(nèi)，每隔1小時測量一次溫度，測得的溫度依次為：5，8，9，15，25，29，31，30，22，25，27，24。試估計每隔1/10小時的溫度值。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’)第28頁，課件共107頁，創(chuàng)作于2023年2月第29頁，課件共107頁，創(chuàng)作于2023年2月xy機翼下輪廓線例已知飛機下輪廓線上數(shù)據(jù)如下，求x每改變0.1時的y值。ToMATLAB(plane)返回第30頁，課件共107頁，創(chuàng)作于2023年2月第31頁，課件共107頁，創(chuàng)作于2023年2月二維插值的定義xyO第一種（網(wǎng)格節(jié)點）：第32頁，課件共107頁，創(chuàng)作于2023年2月

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

構(gòu)造一個二元函數(shù)通過全部已知節(jié)點,即再用計算插值，即第33頁，課件共107頁，創(chuàng)作于2023年2月第二種（散亂節(jié)點）：yx0第34頁，課件共107頁，創(chuàng)作于2023年2月已知n個節(jié)點其中互不相同，

構(gòu)造一個二元函數(shù)通過全部已知節(jié)點,即再用計算插值，即返回第35頁，課件共107頁，創(chuàng)作于2023年2月

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

二維或高維情形的最鄰近插值，與被插值點最鄰近的節(jié)點的函數(shù)值即為所求。返回第36頁，課件共107頁，創(chuàng)作于2023年2月

將四個插值點（矩形的四個頂點）處的函數(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第37頁，課件共107頁，創(chuàng)作于2023年2月插值函數(shù)為：第二片(上三角形區(qū)域)：(x,y)滿足插值函數(shù)為：注意：(x,y)當(dāng)然應(yīng)該是在插值節(jié)點所形成的矩形區(qū)域內(nèi)。顯然，分片線性插值函數(shù)是連續(xù)的；分兩片的函數(shù)表達(dá)式如下：第一片(下三角形區(qū)域)：(x,y)滿足返回第38頁，課件共107頁，創(chuàng)作于2023年2月

雙線性插值是一片一片的空間二次曲面構(gòu)成。雙線性插值函數(shù)的形式如下：其中有四個待定系數(shù)，利用該函數(shù)在矩形的四個頂點（插值節(jié)點）的函數(shù)值，得到四個代數(shù)方程，正好確定四個系數(shù)。雙線性插值xy(x1,y1)(x1,y2)(x2,y1)(x2,y2)O返回第39頁，課件共107頁，創(chuàng)作于2023年2月

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

最鄰近插值‘linear’

雙線性插值‘cubic’

雙三次插值缺省時,雙線性插值第40頁，課件共107頁，創(chuàng)作于2023年2月例：測得平板表面3*5網(wǎng)格點處的溫度分別為：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個單位的地方進(jìn)行插值.第41頁，課件共107頁，創(chuàng)作于2023年2月第42頁，課件共107頁，創(chuàng)作于2023年2月再輸入以下命令:xi=1:0.2:5;yi=1:0.2:3;zi=interp2(x,y,temps,xi',yi,'cubic');mesh(xi,yi,zi)畫出插值后的溫度分布曲面圖.ToMATLAB(wendu)第43頁，課件共107頁，創(chuàng)作于2023年2月第44頁，課件共107頁，創(chuàng)作于2023年2月例>>[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),第45頁，課件共107頁，創(chuàng)作于2023年2月選較密的插值點，用默認(rèn)的線性插值算法進(jìn)行插值>>[x1,y1]=meshgrid(-3:.2:3,-2:.2:2);>>z0=interp2(x,y,z,x1,y1);>>surf(x1,y1,z0)第46頁，課件共107頁，創(chuàng)作于2023年2月立方和樣條插值：>>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])第47頁，課件共107頁，創(chuàng)作于2023年2月算法誤差比較>>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))第48頁，課件共107頁，創(chuàng)作于2023年2月山區(qū)地形地貌圖已知某處山區(qū)地形選點測量坐標(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第49頁，課件共107頁，創(chuàng)作于2023年2月山區(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個yi=linspace(0,6,80);%加密縱坐標(biāo)數(shù)據(jù)到60個[xii,yii]=meshgrid(xi,yi);%生成網(wǎng)格數(shù)據(jù)zii=interp2(x,y,z,xii,yii,'cubic');%插值mesh(xii,yii,zii)%加密后的地貌圖第50頁，課件共107頁，創(chuàng)作于2023年2月山區(qū)地形地貌圖結(jié)果第51頁，課件共107頁，創(chuàng)作于2023年2月

通過此例對最近鄰點插值、雙線性插值方法和雙三次插值方法的插值效果進(jìn)行比較。ToMATLAB(moutain)返回第52頁，課件共107頁，創(chuàng)作于2023年2月第53頁，課件共107頁，創(chuàng)作于2023年2月第54頁，課件共107頁，創(chuàng)作于2023年2月第55頁，課件共107頁，創(chuàng)作于2023年2月第56頁，課件共107頁，創(chuàng)作于2023年2月第57頁，課件共107頁，創(chuàng)作于2023年2月

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

cz

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

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

最鄰近插值‘linear’

雙線性插值‘cubic’

雙三次插值'v4'-Matlab提供的插值方法缺省時,雙線性插值第58頁，課件共107頁，創(chuàng)作于2023年2月>>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ū)域隨機選擇一組坐標(biāo)，用’v4’與’cubic’插值法進(jìn)行處理，并對誤差進(jìn)行比較。第59頁，課件共107頁，創(chuàng)作于2023年2月第60頁，課件共107頁，創(chuàng)作于2023年2月誤差分析>>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])第61頁，課件共107頁，創(chuàng)作于2023年2月第62頁，課件共107頁，創(chuàng)作于2023年2月例：

在x為［3，3］，y為［－2，2］矩形區(qū)域隨機選擇一組坐標(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')%樣本點的二維分布>>figure,plot3(x,y,z,'x'),axis([-3,3,-2,2,-0.7,1.5])，gridon第63頁，課件共107頁，創(chuàng)作于2023年2月去除在(-1,-1/2)點為圓心，以0.5為半徑的圓內(nèi)的點。>>x=-3+6*rand(200,1);y=-2+4*rand(200,1);%重新生成樣本點>>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)的點坐標(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)第64頁，課件共107頁，創(chuàng)作于2023年2月新樣本點擬合曲面>>[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])第65頁，課件共107頁，創(chuàng)作于2023年2月誤差分析>>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)％誤差的二維等高線圖第66頁，課件共107頁，創(chuàng)作于2023年2月命令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ù)一致。第67頁，課件共107頁，創(chuàng)作于2023年2月例：通過函數(shù)生成一些網(wǎng)格型樣本點，試根據(jù)樣本點進(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第68頁，課件共107頁，創(chuàng)作于2023年2月slice(x0,y0,z0,V1,[-0.5,0.3,0.9],[0.6,-0.1],[-1,-0.5,0.5,1])title('SlivesforFourDimFigures')第69頁，課件共107頁，創(chuàng)作于2023年2月[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);第70頁，課件共107頁，創(chuàng)作于2023年2月

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

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)第75頁，課件共107頁，創(chuàng)作于2023年2月海底曲面圖第76頁，課件共107頁，創(chuàng)作于2023年2月%程序二：插值并作出水深小于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)

第77頁，課件共107頁，創(chuàng)作于2023年2月水深小于5的海域范圍第78頁，課件共107頁，創(chuàng)作于2023年2月拉格朗日插值法拉格朗日插值法是基于基函數(shù)的插值方法，插值多項式可表示為其中稱為i次基函數(shù)：第79頁，課件共107頁，創(chuàng)作于2023年2月在MATLAB中編程實現(xiàn)拉格朗日插值法函數(shù)為：Language。功能：求已知數(shù)據(jù)點的拉格朗日多項式；調(diào)用格式：f=

Language(x,y)或f=

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

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

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

f為求得的拉格朗日多項式或x0處的插值。第80頁，課件共107頁，創(chuàng)作于2023年2月

functionf=Language(x,y,x0)%求已知數(shù)據(jù)點的拉格朗日多項式%已知數(shù)據(jù)點的x坐標(biāo)向量:x%已知數(shù)據(jù)點的y坐標(biāo)向量:y%為插值點的x坐標(biāo):x0%求得的拉格朗日多項式或x0處的插值:fsymst;if(length(x)==length(y))n=length(x);elsedisp('x和y的維數(shù)不相等！');return;end%檢錯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));%計算拉格朗日基函數(shù)

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

simplify(f);%化簡

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

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

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

endendend第81頁，課件共107頁，創(chuàng)作于2023年2月例4-3根據(jù)下表的數(shù)據(jù)點求出其拉格朗日插值多項式，并計算當(dāng)x=1.6時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第82頁，課件共107頁，創(chuàng)作于2023年2月利用均差的牛頓插值法函數(shù)f的零階均差定義為，一階定義均差為一般地，函數(shù)f的k階均差定義為：利用均差的牛頓插值法多項式為

第83頁，課件共107頁，創(chuàng)作于2023年2月系數(shù)的計算過程如表所示表

均差計算表格一階均差二階均差三階均差……n階均差……………第84頁，課件共107頁，創(chuàng)作于2023年2月在MATLAB中編程實現(xiàn)利用均差牛頓插值法函數(shù)為：Newton。功能：求已知數(shù)據(jù)點的均差形式的牛頓插值多項式；調(diào)用格式：f=

Newton(x,y)或f=

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

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

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

f為求得的牛頓插值法多項式或x0處的插值。第85頁，課件共107頁，創(chuàng)作于2023年2月functionf=Newton(x,y,x0)%求已知數(shù)據(jù)點的均差形式牛頓插值多項式%已知數(shù)據(jù)點的x坐標(biāo)向量:x%已知數(shù)據(jù)點的y坐標(biāo)向量:y%為插值點的x坐標(biāo):x0%求得的均差形式牛頓插值多項式或x0處的插值:fsymst;if(length(x)==length(y))n=length(x);c(1:n)=0.0;elsedisp('x和y的維數(shù)不相等！');return;end第86頁，課件共107頁，創(chuàng)作于2023年2月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);%將插值多項式展開

f=vpa(f,6);endendend第87頁，課件共107頁，創(chuàng)作于2023年2月例根據(jù)下表的數(shù)據(jù)點求出其均差形式牛頓插值多項式，并計算當(dāng)x=2.0時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第88頁，課件共107頁，創(chuàng)作于2023年2月利用差分的牛頓插值法差分分為向前差分、后向差分和中心差分三種，它們的記法及定義如下為：

其中：代表向前差分；代表向后差分；代表向后差分。第89頁，課件共107頁，創(chuàng)作于2023年2月

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

差分計算表格………………第90頁，課件共107頁，創(chuàng)作于2023年2月向前牛頓插值向前牛頓插值多項式可表示如下：其中叫做步長，，且的取值范圍為。第91頁，課件共107頁，創(chuàng)作于2023年2月在MATLAB中編程實現(xiàn)向前牛頓插值法函數(shù)為：Newtonforward。功能：求已知數(shù)據(jù)點的向前牛頓插值法多項式；調(diào)用格式：f=

Newtonforward(x,y)或

f=

Newtonforward

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

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

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

f為求得的向前牛頓插值法多項式或x0處的插值。第92頁，課件共107頁，創(chuàng)作于2023年2月functionf=Newtonforward(x,y,x0)%求已知數(shù)據(jù)點的向前差分牛頓插值多項式%已知數(shù)據(jù)點的x坐標(biāo)向量:x%已知數(shù)據(jù)點的y坐標(biāo)向量:y%為插值點的x坐標(biāo):x0%求得的向前差分牛頓插值多項式或x0處的插值:fsymst;if(length(x)==length(y))n=length(x);c(1:n)=0.0;elsedisp('x和y的維數(shù)不相等！');return;end第93頁，課件共107頁，創(chuàng)作于2023年2月f=y(1);y1=0;xx=linspace(x(1),x(n),(x(2)-x(1)));if(xx~=x)disp('節(jié)點之間不是等距的！');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第94頁，課件共107頁，創(chuàng)作于2023年2月向前牛頓插值向后牛頓插值多項式可表示如下：其中叫做步長，，且的取值范圍為。第95頁，課件共107頁，創(chuàng)作于2023年2月在MATLAB中編程實現(xiàn)向后牛頓插值法函數(shù)為：Newtonback。功能：求已知數(shù)據(jù)點的向前牛頓插值法多項式；調(diào)用格式：f=

Newtonback

(x,y)或

f=

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

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

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

f為求得的向前牛頓插值法多項式或x0處的插值。第96頁，課件共107頁，創(chuàng)作于2023年2月functionf=Newtonback(x,y,x0)%求已知數(shù)據(jù)點的向后差分牛頓插值多項式%已知數(shù)據(jù)點的x坐標(biāo)向量:x%已知數(shù)據(jù)點的y坐標(biāo)向量:y%為插值點的x坐標(biāo):x0%求得的向前差分牛頓插值多項式或x0處的插值:fsymst;if(length(x)==length(y))n=length(x);c(1:n)=0.0;elsedisp('x和y的維數(shù)不相等！');return;end第97頁，課件共107頁，創(chuàng)作于2023年2月f=y(n);y1=0;xx=linspace(x(1),x(n),(x(2)-x(1)));if(xx~=x)disp('節(jié)點之間不是等距的！');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第98頁，課件共107頁，創(chuàng)作于2023年2月例5根據(jù)下表的數(shù)據(jù)點求出其差分形式的牛頓插值多項式，并計算當(dāng)x=1.55時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第99頁，課件共107頁，創(chuàng)作于2023年2月Hermite插值Hermite插值滿足在節(jié)點上等于給定函數(shù)值，而且在節(jié)點上的導(dǎo)數(shù)值也等于給定的導(dǎo)數(shù)值。對于高階導(dǎo)數(shù)的情況，Hermite插值多項式比較復(fù)雜，在實際中，常常遇到的是函數(shù)值與一階導(dǎo)數(shù)給定的情況。在此情況下，n個節(jié)點x1，x2，…，xn的Hermite插值多項式的表達(dá)式如下：其中，，，第100頁，課件共107頁，創(chuàng)作于2023年2月在MATLAB中編程實現(xiàn)Hermite插值法函數(shù)為：Hermite。功能：求已知數(shù)據(jù)點的Hermite插值法多項式；調(diào)用格式：f=

Hermite

(x,y,y_1)或

f=

Hermite

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

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

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

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

f為求得的Hermite插值法多項式或x0處的插值。第101頁，課件共107頁，創(chuàng)作于2023年2月functionf=Hermite(x,y,y_1,x0)%求已知數(shù)據(jù)點的向后差分牛頓插值多項式%已知數(shù)據(jù)點的x坐標(biāo)向量:x%已知數(shù)據(jù)點的y坐標(biāo)向量:y%已知數(shù)據(jù)點的導(dǎo)數(shù)向量:y_1%求得的Hermite插值多項式或x0處的插值:fsymst;f=0.0;if(length(x)==length(y))if(length(y)==length(y_