插值與擬合§7.1引言

在工程和科學(xué)實驗中，常常需要從一組實驗觀測數(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ù)與要求的不同而異，通常可采用兩種方法：插值與數(shù)據(jù)擬合?！?.1.1插值引例7.1.1已經(jīng)測得在北緯32.3

海洋不同深度處的溫度如下表：表7.1.1根據(jù)這些數(shù)據(jù)，我們希望能合理地估計出其它深度（如500米、600米、1000米…）處的水溫。解決這個問題，可以通過構(gòu)造一個與給定數(shù)據(jù)相適應(yīng)的函數(shù)來解決，這是一個被稱為插值的問題。深度x(m)46671495014221634水溫y(C)7.044.283.402.542.13解決這個問題，可以通過構(gòu)造一個與給定數(shù)據(jù)相適應(yīng)的函數(shù)來解決，這是一個被稱為插值的問題。

插值問題的基本提法：對于給定的函數(shù)表其中f(x)在區(qū)間[a,b]上連續(xù)，x0，x1，…，xn為[a,b]上n1個互不相同的點，要求在一個性質(zhì)優(yōu)良、便于計算的函數(shù)類{P(x)}中，選出一個使P(xi

)=yi，i=0,1,…,n(7.1.1)成立的函數(shù)P(x)作為f(x)的近似，這就是最基本的插值問題（見圖7.1.1）。xx0x1…xny=f(x)y0y1…yn為便于敘述，通常稱區(qū)間[a,b]為插值區(qū)間，稱點x0，x1，…，xn為插值節(jié)點，稱函數(shù)類{P(x)}為插值函數(shù)類，稱式(7.1.1)為插值條件，稱函數(shù)P(x)為插值函數(shù)，稱f(x)為被插函數(shù)。求插值函數(shù)P(x)的方法稱為插值法。圖7.1.1插值問題示意圖

引例7.1.2在某化學(xué)反應(yīng)中，已知生成物的濃度與時間有關(guān)。今測得一組數(shù)據(jù)如下：

表7.1.2根據(jù)這些數(shù)據(jù)，我們希望尋找一個y=f(t)的近似表達(dá)式（如建立濃度y與時間t之間的經(jīng)驗公式等）。從幾何上看，就是希望根據(jù)給定的一組點（1,4.00），…,（16,10.60），求函數(shù)y=f(t)的圖象的一條擬合曲線。時間t（分）12345678濃度y1034.006.408.008.809.229.509.709.86時間t（分）910111213141516濃度y10310.0010.2010.3210.3210.5010.5510.5810.60

數(shù)據(jù)擬合問題的基本提法：對于給定的函數(shù)表其中f(x)在區(qū)間[a,b]上連續(xù)，x0，x1，…，xn為[a,b]上n1個互不相同的點，要求找一個簡單合理的函數(shù)近似表達(dá)式

(x)，使

(x)與f(x)在某種準(zhǔn)則下最為接近，這就是最基本的數(shù)據(jù)擬合問題（見圖7.1.2）。通常，我們稱

(x)為給定數(shù)據(jù)點的擬合函數(shù)。xx0x1…xny=f(x)y0y1…yn圖7.1.2數(shù)據(jù)擬合問題示意圖

§7.1.3插值與數(shù)據(jù)擬合的基本理論依據(jù)插值方法與數(shù)據(jù)擬合的基本理論依據(jù)，就是數(shù)學(xué)分析中的Weierstrass定理：設(shè)函數(shù)f(x)在區(qū)間[a,b]上連續(xù)，則對

>0，存在多項式P(x)，使得即：有界區(qū)間上的連續(xù)函數(shù)被多項式一致逼近?！?.1.4實際應(yīng)用中兩種方法的選擇在實際應(yīng)用中，究竟選擇哪種方法比較恰當(dāng)？總的原則是根據(jù)實際問題的特點來決定采用哪一種方法。具體說來，可從以下兩方面來考慮：1．如果給定的數(shù)據(jù)是少量的且被認(rèn)為是嚴(yán)格精確的，那么宜選擇插值方法。采用插值方法可以保證插值函數(shù)與被插函數(shù)在插值節(jié)點處完全相等。2．如果給定的數(shù)據(jù)是大量的測試或統(tǒng)計的結(jié)果，并不是必須嚴(yán)格遵守的，而是起定性地控制作用的，那么宜選用數(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ù)分段多項式插值就是求一個分段（共n段）多項式P(x)，使其滿足P(xi)=yi（i=0,1,…,n）或更高的要求。一般地，分段多項式插值中的多項式都是低次多項式（不超過三次）。xx0x1…xny=f(x)y0y1…yn1．分段線性插值

分段線性插值函數(shù)P1(x)是一個分段一次多項式（分段線性函數(shù)）。在幾何上就是用折線代替曲線，如圖7.2.1，故分段線性插值亦稱為折線插值。其插值公式為

其中x[xi,xi+1]

圖7.2.1分段線性插值示意圖

2．分段二次插值分段二次插值函數(shù)P2(x)是一個分段二次多項式。在幾何上就是分段拋物線代替曲線y=f(x)，故分段二次插值又稱為分段拋物插值。其插值公式其中x[xi-1

,xi+1]

3．三次Hermite插值三次Hermite插值問題的基本提法一：已知一維數(shù)據(jù)求一個三次多項式P3(x)，使之滿足P3(xi)=yi，P3(xi)=mi，i=0,1(7.2.3)xx0x1y=f(x)y0y1y

=f(x)m0m1下面的(7.2.5)、(7.2.6)兩式構(gòu)成里三次Hermite插值基本提法一的插值公式P3(x)=0(x)y0

1(x)y1

0(x)m0

1(x)m1(7.2.5)三次Hermite插值問題的基本提法二：已知一維數(shù)據(jù)求一個三次多項式P3(x)，使之滿足P3(xi)=yi，i=0,1,2，P3(x1)=mi(7.2.3)xx0x1x2y=f(x)y0y1y2y

=f(x)m1下面的(7.2.9)、(7.2.10)兩式構(gòu)成里三次Hermite插值基本提法二的插值公式P3(x)=0(x)y0

1(x)y1

0(x)m0

1(x)m1(7.2.9)§7.2.2一維數(shù)據(jù)的三次樣條插值上述介紹的分段多項式插值，其優(yōu)點為計算簡單、穩(wěn)定性好、收斂性有保證，且易于在計算機(jī)上實現(xiàn)。但它也明顯存在著缺陷。它只能保證在每個小區(qū)間段[xi,xi+1]內(nèi)光滑，在各小區(qū)間連接點xi

處連續(xù)，卻不能保證整條曲線的光滑、光順性，難以滿足某些工程的要求。對于象高速飛機(jī)的機(jī)翼形線，船體放樣等型值線往往要求有二階光滑度，即有二階連續(xù)導(dǎo)數(shù)。而由60年代開始，首先起源與航空、造船業(yè)等工程設(shè)計的實際需要而發(fā)展起來的樣條插值，既保留了分段多項式插值的各種優(yōu)點，又提高了插值函數(shù)的光滑度。在此，僅介紹應(yīng)用最廣且具有二階連續(xù)導(dǎo)數(shù)的三次樣條插值方法。

1．三次樣條插值問題的基本提法對于給定的一維數(shù)據(jù)求一個三次多項式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有xx0x1…xny=f(x)y0y1…yn

2．三次樣條插值函數(shù)給定區(qū)間[a,b]的一個劃分：a=x0<x1<…<xn=b，設(shè)函數(shù)y=f(x)在節(jié)點xi

上的值為yi=f(xi)，i=0,1,…,n。如果S(x)于[a,b]有二階連續(xù)導(dǎo)數(shù)，且在每個小區(qū)間[xi

,xi+1]上是三次多項式，則稱S(x)是節(jié)點x0，x1，…，xn上的三次樣條函數(shù)。如果S(x)在節(jié)點xi上還滿足插值條件S(xi)=yi，i=0,1,…,n，(7.2.11)則稱S(x)為三次樣條插值函數(shù)。對應(yīng)于劃分的三次樣條插值函數(shù)的表達(dá)式為其中

3．邊界條件在式(7.2.12)給出的三次多項式中，共含有n3個待定系數(shù)。而由插值條件(7.2.11)式，可列出n1個方程，方程組中未知數(shù)的個數(shù)比方程個數(shù)多2，還需附加2個條件才能進(jìn)行求解。通?？稍趨^(qū)間端點x0=a和xn=b處各附加一個條件（稱為邊界條件或邊值條件）去確定S(x)。邊界條件類型很多，較基本而又常見的有三類：(1)第一邊值條件，即給出邊界點的一階導(dǎo)數(shù)值S(x0)=y0，

S(xn)=yn(7.2.13)

(2)第二邊值條件，即給出邊界點的二階導(dǎo)數(shù)值S(x0)=y0，

S(xn)=yn(7.2.14)特別地，當(dāng)S(x0)=S(xn)=0時，稱為自然邊界條件。滿足自然邊界條件的三次樣條插值函數(shù)稱為自然樣條插值函數(shù)。(3)第三邊值條件（混合邊值條件）

其中1、2、1、2、1、2

為定數(shù)。當(dāng)1、2為零時，則為第一邊值條件，當(dāng)1、2

為零時，則為第二邊值條件。§7.2.3一維數(shù)據(jù)插值的Matlab實現(xiàn)Matlab中一維數(shù)據(jù)的插值函數(shù)為：interp1()。其調(diào)用格式為yi=interp1(x,y,xi,'methos'),其中，x，y——為插值節(jié)點，均為向量；

xi——任取的被插值點可以是一個數(shù)值，也可以是一個向量；

yi——為被插值點xi處的插值結(jié)果；

‘methos’——為采用的插值方法：'nearest'：表示最臨近插值，

'linear'：表示線性插值，

'cubic'：表示三次插值，

'spline'：表示三次樣條插值。注意：（1）上述‘methos’中所有的插值方法都要求x是單調(diào)的，并且xi不能超過x的取值范圍；（2）三次樣條插值函數(shù)的調(diào)用格式有兩種：yi=interp1(x,y,xi,‘spline’)yi=spline(x,y,xi)它們是等價的。注意：（1）上述‘methos’中所有的插值方法都要求x是單調(diào)的，并且xi不能超過x的取值范圍；（2）三次樣條插值函數(shù)的調(diào)用格式有兩種：yi=interp1(x,y,xi,‘spline’)yi=spline(x,y,xi)它們是等價的?！?.3二維數(shù)據(jù)的基本插值方法簡介

對于二維數(shù)據(jù)的插值，首先要考慮兩個問題：一是二維區(qū)域是任意區(qū)域還是規(guī)則區(qū)域，二是給定的數(shù)據(jù)是有規(guī)律分布的還是散亂的、隨機(jī)分布的。第一個問題比較容易處理。目前的插值方法基本上是基于規(guī)則區(qū)域的，對于不規(guī)則區(qū)域，只需將其，劃分為規(guī)則區(qū)域或擴(kuò)充為規(guī)則區(qū)域來討論即可。對于第二個問題，當(dāng)給定的數(shù)據(jù)是有規(guī)律分布時，方法較多也較成熟；而給定的數(shù)據(jù)是散亂的、隨機(jī)分布時，沒有固定的方法，但一般的處理思想是：從給定的數(shù)據(jù)出發(fā)，依據(jù)一定的規(guī)律恢復(fù)出規(guī)則分布點上的數(shù)據(jù)，轉(zhuǎn)化為數(shù)據(jù)分布有規(guī)律的情形來處理。

二維數(shù)據(jù)插值的方法也有很多。在此，針對給定數(shù)據(jù)有規(guī)律分布和散亂分布兩種情形，簡單介紹雙三次樣條插值方法和改進(jìn)的Shepard方法（反距離平方法）的基本概念和基本思想，及其Matlab實現(xiàn)。§7.3.1雙三次樣條插值雙三次樣條插值方法，是用來解決規(guī)則區(qū)域上給定數(shù)據(jù)有規(guī)律分布的插值問題的常用方法。實際上，雙三次樣條函數(shù)是由兩個一維三次樣條函數(shù)作直積產(chǎn)生的。對任意固定的y0[c,d]，S(x,y0)是關(guān)于x的三次樣條函數(shù)，同理，對任意固定的x0[a,b]，S(x0,y)是關(guān)于y的三次樣條函數(shù)。從而，根據(jù)一維三次樣條函數(shù)的算法可以設(shè)計出S(x,y)的具體算法?！?.3.2改進(jìn)的Shepard方法改進(jìn)的Shepard方法，也稱反距離加權(quán)平均法，這是解決規(guī)則區(qū)域上給定數(shù)據(jù)散亂、隨機(jī)分布的插值問題的一個常用的方法。問題：設(shè)T=[a,b][c,d]上散亂分布N個點V1，V2，…，VN

，其中Vk

=(xk，yk)處給出數(shù)據(jù)fk

，k=1,2,…,N。要尋求T上的二元函數(shù)F(x，y)，使F(xk，yk)=fk

，k=1,2,…,N。一個典型的容易想到的方法是“反距離加權(quán)平均”方法，又稱之為Shepard方法。這方法的基本思想是，在非給定數(shù)據(jù)的點處，定義其函數(shù)值由已知數(shù)據(jù)點與該點距離的近或遠(yuǎn)作加權(quán)平均決定。按照上述的思想，可從給定的數(shù)據(jù)恢復(fù)出規(guī)則分布點上的數(shù)據(jù)，接下來就可應(yīng)用雙三次樣條插值或其它的二維數(shù)據(jù)插值方法來處理?！?.3.3二維數(shù)據(jù)插值的Matlab實現(xiàn)1．規(guī)則區(qū)域上給定數(shù)據(jù)有規(guī)律分布的二維插值數(shù)據(jù)形式為：y1y2…ynx1x11z12…z1nx2z21z22…z2n……………xmzm1zm2…zmn

插值函數(shù)為：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ù)散亂或隨機(jī)分布的二維插值數(shù)據(jù)形式為：(x1,y1)(x2,y2)…(xn,yn)z1z2…zn

插值函數(shù)為：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．實例例7.3.1氣旋變化情況的可視化表7.3.1是氣象學(xué)家測量得到的氣象資料，它們分別表示在南半球地區(qū)按不同緯度、不同月份的平均氣旋數(shù)字。根據(jù)這些數(shù)據(jù)，繪制出氣旋分布曲面圖形。0—1010—2020—3030—4040—5050—6060—7070—8080—901月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.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.4一維數(shù)據(jù)擬合的最小二乘法簡介

數(shù)據(jù)擬合問題的形式多種多樣，解決的方法也有許多。在此，我們只簡單介紹以最小二乘法為準(zhǔn)則的一維數(shù)據(jù)擬合方法。§7.4.1數(shù)據(jù)擬合的最小二乘法簡介對于給定的一組測量數(shù)據(jù)設(shè)y=

(x)為一擬合函數(shù)，記i

=(xi)

yi(i=1,2,…,n)，(7.4.1)則稱i

為擬合函數(shù)(x)在xi

點處的偏差或殘量。xx1x2…xny=f(x)y1y2…yn為使(x)在整體上盡可能與給定數(shù)據(jù)的函數(shù)f(x)近似，我們常采用偏差的平方和達(dá)到最小，即來保證每個偏差的絕對值|i|都很小。這一原則稱為最小二乘原則，根據(jù)最小二乘原則確定擬合函數(shù)(x)的方法稱為最小二乘法。

1.線性最小二乘擬合我們知道，函數(shù)系{xk|k=0,1,…,m}的線性組合(x)=a0

a1xa2x2

…amxm

為m次多項式。一般地，若函數(shù)系{k(x)|k=0,1,…,m}是線性無關(guān)的，則其線性組合稱為函數(shù)系{k(x)|k=0,1,…,m}的廣義多項式。如三角多項式就是函數(shù)系{1,cosx,sinx,cos2x,sin2x,…,cosmx,sinmx}的廣義多項式。設(shè){k(x)|k=0,1,…,m}為一線性無關(guān)的函數(shù)系，取擬合函數(shù)為(7.4.3)式給出的廣義多項式，使得(7.4.2)成立。由于

(x)的待定系數(shù)a0,a1,a2,…,am

全部以線性形式出現(xiàn)，故我們稱之為線性最小二乘擬合。在式(7.4.2)中，目標(biāo)函數(shù)S是關(guān)于參數(shù)a0,a1,a2,…,am的多元函數(shù)，由多元函數(shù)取得最小值的必要條件知，欲使S達(dá)到極小，須滿足即亦即其中k=0,1,…,m，式(7.4.4)是關(guān)于a0,a1,a2,…,am

的線性方程組，稱為正規(guī)方程組。從正規(guī)方程組(7.4.4)中解出a0,a1,a2,…,am，于是就得到了最小二乘擬合函數(shù)(x)。

2．非線性最小二乘擬合如果擬合函數(shù)(x)=(x,a0,a1,a2,…,am)的待定參數(shù)a0,a1,a2,…,am不能全部以線性形式出現(xiàn)，如指數(shù)擬合函數(shù)等，這便是非線性最小二乘擬合問題。一般地，非線性最小二乘擬合問題是一個非線性函數(shù)的極小化問題，可用非線性優(yōu)化方法求解。

3．最小二乘擬合函數(shù)的選擇最小二乘法中，擬合函數(shù)的選擇是很重要的?？梢酝ㄟ^對給定數(shù)據(jù)的分析來選擇，也可以直接由實際問題給定。最常用的是多項式和樣條函數(shù)，尤其是當(dāng)不知道該選擇什么樣的擬合函數(shù)時，通常可以考慮選擇樣條函數(shù)。另外，對同一問題，也可選擇不同的函數(shù)進(jìn)行最小二乘擬合，比較各自誤差的大小，從中選出誤差較小的作為擬合函數(shù)?！?.4.2一維數(shù)據(jù)擬合的Matlab實現(xiàn)

1.多項式函數(shù)擬合擬合函數(shù)的命令為：polyfit()，其調(diào)用格式為：a=polyfit(xdata,ydata,m)，其中

m——表示多項式的最高階數(shù)；

xdata，ydata——為將要擬合的數(shù)據(jù)，它們都是以數(shù)組方式輸入；

a——輸出參數(shù)，為擬合多項式的系數(shù)a=[a0,a1,a2,…,am]。多項式在x處的值y可用如下命令格式計算：y=polyval(a,x)。

2.一般的曲線擬合擬合函數(shù)的命令為：curvefit()，或lsqcurvefit()，其調(diào)用格式為p=curvefit(‘Fun’,p0,xdata,ydata)，或p=lsqcurvefit('Fun',p0,xdata,ydata)，其中

Fun——表示函數(shù)Fun(p,xdta)的M文件；

P0——為函數(shù)的初值。若要求點x處的函數(shù)值y，可用程序f=Fun(p,x)計算。

3．實例例7.4.1求解§7.1.2引例7.1.2。解：(1)選擇擬合曲線作出給定數(shù)據(jù)的散點圖如下：通過對散點圖的分析可以看出，數(shù)據(jù)點的分布為一條單調(diào)上升的曲線。具有這種特性的曲線很多，我們選取如下三種函數(shù)來擬合：①多項式(x)=a0

a1x

…

amxm，m為適當(dāng)選取的正整數(shù)；②；

③（a>0,b<0）。

(2)擬合運(yùn)算首先，分別用二、三、六次多項式擬合，計算得輸出參數(shù)分別為p1=[0.0445，1.0711，4.3252]p2=[0.0060，0.1963，2.1346，2.5952]p3=[0.0000，0.0004，0.0103，0.1449，

1.1395，4.9604，0.0498]擬合函數(shù)分別為(1)(x)=0.04451.0711x

4.3252x2(2)(x)=0.00600.1963x

2.1346x2

2.5952x3(3)(x)=0.0004x

0.0103x2

0.1449x3

1.1395x4

4.9304