插值法與數(shù)據(jù)擬合法[分享借鑒]

1、第七講 插值方法與數(shù)據(jù)擬合 7.1 引言在工程和科學實驗中，常常需要從一組實驗觀測數(shù)據(jù) (xi , yi ) (i = 1, 2, , n) 揭示自變量x與因變量y之間的關(guān)系，一般可以用一個近似的函數(shù)關(guān)系式y(tǒng) = f (x) 來表示。函數(shù)f (x) 的產(chǎn)生辦法因觀測數(shù)據(jù)與要求的不同而異，通?？刹捎脙煞N方法：插值與數(shù)據(jù)擬合。 7.1.1 插值方法1 引例1 已經(jīng)測得在北緯32.3 海洋不同深度處的溫度如下表：表7.1.1深度x (m)46671495014221634水溫y (C)7.044.283.402.542.13根據(jù)這些數(shù)據(jù)，我們希望能合理地估計出其它深度（如500米、600米、1000

2、米）處的水溫。 解決這個問題，可以通過構(gòu)造一個與給定數(shù)據(jù)相適應(yīng)的函數(shù)來解決，這是一個被稱為插值的問題。 2插值問題的基本提法對于給定的函數(shù)表xx0x1xny = f (x)y0y1yn其中f (x) 在區(qū)間 a, b 上連續(xù)，x0，x1，xn為 a, b 上n + 1個互不相同的點，要求在一個性質(zhì)優(yōu)良、便于計算的函數(shù)類 P(x) 中，選出一個使P(xi ) = yi，i = 0, 1, , n (7.1.1)成立的函數(shù)P(x) 作為 f (x) 的近似，這就是最基本的插值問題（見圖7.1.1）。為便于敘述，通常稱區(qū)間 a, b 為插值區(qū)間，稱點x0，x1，xn為插值節(jié)點，稱函數(shù)類 P(x) 為

3、插值函數(shù)類，稱式 (7.1.1) 為插值條件，稱函數(shù) P(x) 為插值函數(shù)，稱f (x) 為被插函數(shù)。求插值函數(shù) P(x) 的方法稱為插值法。 7.1.2 數(shù)據(jù)擬合1 引例2 在某化學反應(yīng)中，已知生成物的濃度與時間有關(guān)。今測得一組數(shù)據(jù)如下： 表7.1.2時間t（分）12345678濃度y10-34.006.408.008.809.229.509.709.86時間t（分）910111213141516濃度y10-310.0010.2010.3210.3210.5010.5510.5810.60根據(jù)這些數(shù)據(jù)，我們希望尋找一個y = f (t) 的近似表達式（如建立濃度y與時間t之間的經(jīng)驗公式等）。

4、從幾何上看，就是希望根據(jù)給定的一組點（1, 4.00），,（16, 10.60），求函數(shù)y = f (t) 的圖象的一條擬合曲線。2 數(shù)據(jù)擬合問題的基本提法對于給定的函數(shù)表xx0x1xny = f (x)y0y1yn其中f (x) 在區(qū)間 a, b 上連續(xù)，x0，x1，xn為 a, b 上n + 1個互不相同的點，要求找一個簡單合理的函數(shù)近似表達式 j (x)，使 j (x) 與f (x) 在某種準則下最為接近，這就是最基本的數(shù)據(jù)擬合問題（見圖7.1.2）。 通常，我們稱 j (x) 為給定數(shù)據(jù)點的擬合函數(shù)。 圖7.1.1 插值問題示意圖圖 7.1.2 數(shù)據(jù)擬合問題示意圖 7.1.3 插值方法

5、與數(shù)據(jù)擬合的基本理論依據(jù)插值方法與數(shù)據(jù)擬合的基本理論依據(jù)，就是數(shù)學分析中的Weierstrass定理：設(shè)函數(shù)f (x) 在區(qū)間 a, b 上連續(xù)，則對 e 0，存在多項式P(x)，使得。即：有界區(qū)間上的連續(xù)函數(shù)被多項式一致逼近。 7.1.4 實際應(yīng)用中兩種方法的選擇在實際應(yīng)用中，究竟選擇哪種方法比較恰當？總的原則是根據(jù)實際問題的特點來決定采用哪一種方法。具體說來，可從以下兩方面來考慮：1如果給定的數(shù)據(jù)是少量的且被認為是嚴格精確的，那么宜選擇插值方法。采用插值方法可以保證插值函數(shù)與被插函數(shù)在插值節(jié)點處完全相等。2如果給定的數(shù)據(jù)是大量的測試或統(tǒng)計的結(jié)果，并不是必須嚴格遵守的，而是起定性地控制作用的

6、，那么宜選用數(shù)據(jù)擬合的方法。這是因為，一方面測試或統(tǒng)計數(shù)據(jù)本身往往帶有測量誤差，如果要求所得的函數(shù)與所給數(shù)據(jù)完全吻合，就會使所求函數(shù)保留著原有的測量誤差；另一方面，測試或統(tǒng)計數(shù)據(jù)通常很多，如果采用插值方法，不僅計算麻煩，而且逼近效果往往較差。 7.2 一維數(shù)據(jù)的基本插值方法簡介插值函數(shù)類的取法很多，可以是代數(shù)多項式，也可以是三角多項式或有理函數(shù)；可以是 a, b 上任意光滑函數(shù)，也可以是分段光滑函數(shù)。在此介紹最基本、最常用的兩種插值方法：分段多項式插值與三次樣條插值，及其Matlab實現(xiàn)。 7.2.1 一維數(shù)據(jù)的分段多項式插值對于給定的一維數(shù)據(jù)xx0x1xny = f (x)y0y1yn分段多

7、項式插值就是求一個分段（共n段）多項式P(x)，使其滿足P(xi ) = yi（i = 0, 1, , n）或更高的要求。一般地，分段多項式插值中的多項式都是低次多項式（不超過三次）。1 分段線性插值 y 分段線性插值函數(shù)P1 (x) 是一個分段一次多項式（分段線 f (x)性函數(shù)）。在幾何上就是用折線代替曲線，如圖7.2.1，故分段線性插值亦稱為折線插值。其插值公式為 P(x) ，xxi , xi +1 (7.2.1) 0 x0 x1 xn-1 xn x 2分段二次插值 圖7.2.1 分段線性插值示意圖分段二次插值函數(shù)P2 (x) 是一個分段二次多項式。在幾何上就是分段拋物線代替曲線y =

8、f (x)，故分段二次插值又稱為分段拋物插值。其插值公式為， xxi -1 , xi +1 (7.2.2) 3三次Hermite插值三次Hermite插值問題的基本提法一：已知一維數(shù)據(jù)xx0x1y = f (x)y0y1y = f (x)m0m1求一個三次多項式P3 (x)，使之滿足P3 (xi ) = yi，P3 (xi ) = mi，i = 0, 1 (7.2.3)構(gòu)造三次插值基函數(shù)a0(x)，a1(x)，b0(x)，b1(x)，使之滿足 (7.2.4)利用這四個插值基函數(shù)，取三次多項式P3 (x) 為P3 (x) = a0(x) y0 + a1(x) y1 + b0(x) m0 + b1

9、(x) m1 (7.2.5)將插值條件 (7.2.3) 式代入，可推得： (7.2.6)(7.2.5)、 (7.2.6) 兩式構(gòu)成了三次Hermite插值基本提法一的插值公式。三次Hermite插值問題的基本提法二：已知一維數(shù)據(jù)xx0x1x2y = f (x)y0y1y2y = f (x)m1求一個三次多項式P3 (x)，使之滿足P3 (xi ) = yi，i = 0, 1, 2，P3 (x1 ) = mi (7.2.7)構(gòu)造三次插值基函數(shù)a0(x)，a1(x)，a2 (x)，b1(x)，使之滿足 (7.2.8)利用這四個插值基函數(shù)，取三次多項式P3 (x) 為P3 (x) = a0(x) y

10、0 + a1(x) y1 + a2(x) y2 + b1(x) m1 (7.2.9)將插值條件 (7.2.7) 式代入，可推得： (7.2.10)(7.2.9)、 (7.2.10) 兩式構(gòu)成了三次Hermite插值基本提法二的插值公式。 7.2.2 一維數(shù)據(jù)的三次樣條插值上述介紹的分段多項式插值，其優(yōu)點為計算簡單、穩(wěn)定性好、收斂性有保證，且易于在計算機上實現(xiàn)。但它也明顯存在著缺陷。它只能保證在每個小區(qū)間段 xi , xi +1 內(nèi)光滑，在各小區(qū)間連接點xi 處連續(xù)，卻不能保證整條曲線的光滑、光順性，難以滿足某些工程的要求。對于象高速飛機的機翼形線，船體放樣等型值線往往要求有二階光滑度，即有二階

11、連續(xù)導(dǎo)數(shù)。而由60年代開始，首先起源與航空、造船業(yè)等工程設(shè)計的實際需要而發(fā)展起來的樣條插值，既保留了分段多項式插值的各種優(yōu)點，又提高了插值函數(shù)的光滑度。在此，僅介紹應(yīng)用最廣且具有二階連續(xù)導(dǎo)數(shù)的三次樣條插值方法。1 三次樣條插值問題的基本提法對于給定的一維數(shù)據(jù)xx0x1Xny = f (x)y0y1Yn求一個三次多項式S(x) 滿足條件 （1）S(xi) = yi，i = 0, 1, , n； （2）S(x) 具有二階連續(xù)導(dǎo)數(shù)，特別在節(jié)點xi上應(yīng)滿足連續(xù)性要求，即對i = 0, 1, , n有2 三次樣條插值函數(shù)給定區(qū)間 a, b 的一個劃分D：a = x0 x1 =466&x1(i)714&x

12、1(i)950&x1(i)=1422 a4=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; else a5=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 二維數(shù)據(jù)的基本插值方法簡介對于二維數(shù)據(jù)的插值，首先要考慮兩個問題：一是二維區(qū)域是任意區(qū)域還是規(guī)則區(qū)域，二是給定的數(shù)據(jù)是有規(guī)律分布的還是散亂的、隨機分布的。第一個問

13、題比較容易處理。目前的插值方法基本上是基于規(guī)則區(qū)域的，對于不規(guī)則區(qū)域，只需將其，劃分為規(guī)則區(qū)域或擴充為規(guī)則區(qū)域來討論即可。對于第二個問題，當給定的數(shù)據(jù)是有規(guī)律分布時，方法較多也較成熟；而給定的數(shù)據(jù)是散亂的、隨機分布時，沒有固定的方法，但一般的處理思想是：從給定的數(shù)據(jù)出發(fā)，依據(jù)一定的規(guī)律恢復(fù)出規(guī)則分布點上的數(shù)據(jù)，轉(zhuǎn)化為數(shù)據(jù)分布有規(guī)律的情形來處理。二維數(shù)據(jù)插值的方法也有很多。在此，針對給定數(shù)據(jù)有規(guī)律分布和散亂分布兩種情形，簡單介紹雙三次樣條插值方法和改進的Shepard方法（反距離平方法）的基本概念和基本思想，及其Matlab實現(xiàn)。 7.3.1 雙三次樣條插值 雙三次樣條插值方法，是用來解決規(guī)則區(qū)

14、域上給定數(shù)據(jù)有規(guī)律分布的插值問題的常用方法。設(shè)R：a, bc, d 是xy平面上的一個矩形區(qū)域。在x軸和y軸上分別取定分割Dx：a = x0 x1 xnb，Dy：c = y0 y1 ym 0，令由于w(g) 是可微函數(shù)，使得如下定義的F(x，y) 在性能上有所改善， 其中。 (7.3.3) 按照上述的思想，可從給定的數(shù)據(jù)恢復(fù)出規(guī)則分布點上的數(shù)據(jù)，接下來就可應(yīng)用雙三次樣條插值或其它的二維數(shù)據(jù)插值方法來處理。 7.3.3 二維數(shù)據(jù)插值的Matlab實現(xiàn)1規(guī)則區(qū)域上給定數(shù)據(jù)有規(guī)律分布的二維插值數(shù)據(jù)形式為：y1y 2y nx1x11z12z1nx 2z21z22z2nx mzm1zm2zmn 插值函數(shù)

15、為：interp2( )。其調(diào)用格式為zi = interp2(x, y, z, xi, yi, methos),其中 x，y，z 為插值節(jié)點，均為向量； zi 為被插值點 (xi, yi) 處的插值結(jié)果； methos 為采用的插值方法：nearest：表示最臨近插值， linear：表示雙線性插值， cubic：表示雙三次插值， spline：表示雙三次樣條插值。注意：上述 methos 中所有的插值方法都要求x和y是單調(diào)的網(wǎng)格，x和y可以是等距的也可以是不等距的。2規(guī)則區(qū)域上給定數(shù)據(jù)散亂或隨機分布的二維插值數(shù)據(jù)形式為：(x1, y1)(x2, y2)(xn, yn)z1z2zn插值函數(shù)為

16、：e01sef和e01sff，。通常兩者配合使用，其調(diào)用格式為 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 為插值節(jié)點，均為向量； px(j)，py(i) 為被插值節(jié)點； pf(i, j) 為被插值點 (px(j), py(i) 處的插值結(jié)果；其它輸出參數(shù)涉及插值算法，可以不用了解。e01sef的輸出fnodes和rnw為確定插值的參數(shù)，它們是e01sff需要的輸入?yún)?shù)，因此兩函數(shù)需配合使用。 3例例1 氣旋變化情況的可視化

17、表7.3.1是氣象學家測量得到的氣象資料，它們分別表示在南半球地區(qū)按不同緯度、不同月份的平均氣旋數(shù)字。根據(jù)這些數(shù)據(jù)，繪制出氣旋分布曲面圖形。 表7.3.1010102020303040405050606070708080901月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.

18、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解 下面分別用最鄰近插值、雙線性插值、雙三次插值和雙三次樣條插值，給出不同月份按緯度變化的氣旋值（插值結(jié)果），并作出可視化圖形如下。圖7.3.2 四種插值方法的可視化圖形最鄰近插值的Ma

19、tlab程序為： 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.3 1.6, 21.4, 18.5, 20.1, 28.8, 36.6, 24.2, 5.3, 0 2.4, 16.2, 18.2, 20.5, 27.8, 35.5, 25.5, 5.4, 0 3.2, 9.2, 16.6, 25.1, 37.2, 40, 24.6, 4.9, 0.3 1.0, 2.8, 12.9, 29.2, 40.3, 37.6, 21.1, 4.9, 0 0.5, 1.7, 10.1, 32.6, 41.7, 35

20、.4, 22.2, 7.1, 0 0.4, 1.4, 8.3, 33.0, 46.2, 35, 20.2, 5.3, 0.1 0.2, 2.4, 11.2, 31.0, 39.9, 34.7, 21.2, 7.3, 0.2 0.5, 5.8, 12.5, 28.6, 25.9, 35.7, 22.6, 7, 0.3 0.8, 9.2, 21.1, 32.0, 40.3, 39.5, 28.5, 8.6, 0 2.4, 10.3, 23.9, 28.1, 38.2, 40, 25.3, 6.3, 0.1 3.6, 16, 25.5, 25.6, 43.4, 41.9, 24.3, 6.6, 0.

21、3; xi,yi=meshgrid(1:12,5:1:85); zi=interp2(x,y,z,xi,yi,nearest); figure surf(xi,yi,zi) xlabel(月份),ylabel(緯度),zlabel(氣旋), axis(0 12 0 90 0 50) 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=i

22、nterp2(x,y,z,xi,yi,spline)例2 水道測量數(shù)據(jù)（AMCM86A題）在某海域測得一些點（x, y）處的水深z（單位：英尺）由表7.3.2給出，水深數(shù)據(jù)是在低潮時測得的。船的吃水深度為5英尺，問在矩形區(qū)域 (75，200) (-50，150) 里的哪些地方船要避免進入。表7.3.2 水道水深測量數(shù)據(jù)（單位：英尺）x129.0140.0108.588.0185.5195.0105.5y7.5141.528.0147.022.5137.585.5z4868688x157.5107.577.081.0162.0162.0117.5y-6.5-81.03.056.5-66.584.

23、0-38.5z9988949解 （1）假設(shè) 由題目給出的信息是很少的，除了14個位置的水深之外一無所知。顯然，題目要求我們找出水深不到5英尺的區(qū)域。為了討論方便，下面三個假設(shè)是合理的： 所給數(shù)據(jù)是精確的； 討論區(qū)域的海底曲面是光滑的，更確切地說，可以認為曲面的一階、二階導(dǎo)數(shù)是連續(xù)的。因為我們可以認為討論區(qū)域為淺水海域，由于長期的海水水流作用，形成的是以礫石或沙為主要組成部分的海底，不存在珊瑚礁、水底峽谷、山脊等不可意料的突變地形。 水深是一個按區(qū)域來劃分的變量，在某個位置的水深與其周圍區(qū)域的水深是相互依賴的，但這種依賴作用隨距離的增大而減小。就我們討論的問題來說，每一個給定數(shù)據(jù)點影響周圍的每一

24、個未知點，一個給定數(shù)據(jù)點離未知點越近，作用就越大。 （2）問題分析 根據(jù)假設(shè)，海底曲面是連續(xù)光滑的，不存在珊瑚礁、水底峽谷、山脊等不可意料的突變地形，因而很自然的想法就是用某種光滑的擬合曲面去逼近已知的14個數(shù)據(jù)點或以14個已知的數(shù)據(jù)點為基礎(chǔ)，利用二維插值補充一些點的水深，以求得水深不超過5米的區(qū)域。在此，我們采用二維插值方法，應(yīng)用Matlab程序，作出矩形區(qū)域 (75，200) (-50，150) 范圍內(nèi)的海底地形圖、水深不超過5米的危險區(qū)域的平面圖以及水深不超過5米的危險區(qū)域的海底地貌圖，并求出水深不超過5米的危險海域范圍。 （3）問題求解采用改進的Shepard方法，利用Matlab軟件

25、作出作出矩形區(qū)域 (75，200) (-50，150) 范圍內(nèi)的海底地形圖、水深不超過5米的危險區(qū)域的平面圖以及水深不超過5米的危險區(qū)域的海底地貌圖（見圖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,+)

26、,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);for i=1:ny for j=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)