插值法與數(shù)據(jù)擬合法

第頁第七講插值方法及數(shù)據(jù)擬合§7.1引言在工程和科學(xué)實(shí)驗(yàn)中，常常需要從一組實(shí)驗(yàn)觀測數(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插值方法引例1已經(jīng)測得在北緯32.3海洋不同深度處的溫度如下表：表7.1.1深度x(m)46671495014221634水溫y(C)7.044.283.402.542.13根據(jù)這些數(shù)據(jù)，我們希望能合理地估計(jì)出其它深度（如500米、600米、1000米…）處的水溫。解決這個問題，可以通過構(gòu)造一個及給定數(shù)據(jù)相適應(yīng)的函數(shù)來解決，這是一個被稱為插值的問題。2．插值問題的基本提法對于給定的函數(shù)表xx0x1…xny=f(x)y0y1…yn其中f(x)在區(qū)間[a,b]上連續(xù)，x0，x1，…，xn為[a,b]上n1個互不相同的點(diǎn)，要求在一個性質(zhì)優(yōu)良、便于計(jì)算的函數(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ū)間，稱點(diǎn)x0，x1，…，xn為插值節(jié)點(diǎn)，稱函數(shù)類{P(x)}為插值函數(shù)類，稱式(7.1.1)為插值條件，稱函數(shù)P(x)為插值函數(shù)，稱f(x)為被插函數(shù)。求插值函數(shù)P(x)的方法稱為插值法。§7.1.2數(shù)據(jù)擬合引例2在某化學(xué)反應(yīng)中，已知生成物的濃度及時間有關(guān)。今測得一組數(shù)據(jù)如下：表7.1.2時間t（分）12345678濃度y1034.006.408.008.809.229.509.709.86時間t（分）9116濃度y10310.0010.2010.3210.3210.5010.5510.5810.60根據(jù)這些數(shù)據(jù)，我們希望尋找一個y=f(t)的近似表達(dá)式（如建立濃度y及時間t之間的經(jīng)驗(yàn)公式等）。從幾何上看，就是希望根據(jù)給定的一組點(diǎn)（1,4.00），…,（16,10.60），求函數(shù)y=f(t)的圖象的一條擬合曲線。數(shù)據(jù)擬合問題的基本提法對于給定的函數(shù)表xx0x1…xny=f(x)y0y1…yn其中f(x)在區(qū)間[a,b]上連續(xù)，x0，x1，…，xn為[a,b]上n1個互不相同的點(diǎn)，要求找一個簡單合理的函數(shù)近似表達(dá)式(x)，使(x)及f(x)在某種準(zhǔn)則下最為接近，這就是最基本的數(shù)據(jù)擬合問題（見圖7.1.2）。通常，我們稱(x)為給定數(shù)據(jù)點(diǎn)的擬合函數(shù)。圖7.1.1插值問題示意圖圖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，存在多項(xiàng)式P(x)，使得即：有界區(qū)間上的連續(xù)函數(shù)被多項(xiàng)式一致逼近?！?.1.4實(shí)際應(yīng)用中兩種方法的選擇在實(shí)際應(yīng)用中，究竟選擇哪種方法比較恰當(dāng)？總的原則是根據(jù)實(shí)際問題的特點(diǎn)來決定采用哪一種方法。具體說來，可從以下兩方面來考慮：1．如果給定的數(shù)據(jù)是少量的且被認(rèn)為是嚴(yán)格精確的，那么宜選擇插值方法。采用插值方法可以保證插值函數(shù)及被插函數(shù)在插值節(jié)點(diǎn)處完全相等。2．如果給定的數(shù)據(jù)是大量的測試或統(tǒng)計(jì)的結(jié)果，并不是必須嚴(yán)格遵守的，而是起定性地控制作用的，那么宜選用數(shù)據(jù)擬合的方法。這是因?yàn)椋环矫鏈y試或統(tǒng)計(jì)數(shù)據(jù)本身往往帶有測量誤差，如果要求所得的函數(shù)及所給數(shù)據(jù)完全吻合，就會使所求函數(shù)保留著原有的測量誤差；另一方面，測試或統(tǒng)計(jì)數(shù)據(jù)通常很多，如果采用插值方法，不僅計(jì)算麻煩，而且逼近效果往往較差?！?.2一維數(shù)據(jù)的基本插值方法簡介插值函數(shù)類的取法很多，可以是代數(shù)多項(xiàng)式，也可以是三角多項(xiàng)式或有理函數(shù)；可以是[a,b]上任意光滑函數(shù)，也可以是分段光滑函數(shù)。在此介紹最基本、最常用的兩種插值方法：分段多項(xiàng)式插值及三次樣條插值，及其Matlab實(shí)現(xiàn)。§7.2.1一維數(shù)據(jù)的分段多項(xiàng)式插值對于給定的一維數(shù)據(jù)xx0x1…xny=f(x)y0y1…yn分段多項(xiàng)式插值就是求一個分段（共n段）多項(xiàng)式P(x)，使其滿足P(xi)=yi（i=0,1,…,n）或更高的要求。一般地，分段多項(xiàng)式插值中的多項(xiàng)式都是低次多項(xiàng)式（不超過三次）。分段線性插值y分段線性插值函數(shù)P1(x)是一個分段一次多項(xiàng)式（分段線f(x)性函數(shù)）。在幾何上就是用折線代替曲線，如圖7.2.1，故分段線性插值亦稱為折線插值。其插值公式為P(x)，x[xi,xi+1](7.2.1)0x0x1xn-1xnx2．分段二次插值圖7.2.1分段線性插值示意圖分段二次插值函數(shù)P2(x)是一個分段二次多項(xiàng)式。在幾何上就是分段拋物線代替曲線y=f(x)，故分段二次插值又稱為分段拋物插值。其插值公式為，x[xi-1,xi+1](7.2.2)3．三次Hermite插值三次Hermite插值問題的基本提法一：已知一維數(shù)據(jù)xx0x1y=f(x)y0y1y=f(x)m0m1求一個三次多項(xiàng)式P3(x)，使之滿足P3(xi)=yi，P3(xi)=mi，i=0,1(7.2.3)構(gòu)造三次插值基函數(shù)0(x)，1(x)，0(x)，1(x)，使之滿足(7.2.4)利用這四個插值基函數(shù)，取三次多項(xiàng)式P3(x)為P3(x)=0(x)y01(x)y10(x)m01(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求一個三次多項(xiàng)式P3(x)，使之滿足P3(xi)=yi，i=0,1,2，P3(x1)=mi(7.2.7)構(gòu)造三次插值基函數(shù)0(x)，1(x)，2(x)，1(x)，使之滿足(7.2.8)利用這四個插值基函數(shù)，取三次多項(xiàng)式P3(x)為P3(x)=0(x)y01(x)y12(x)y21(x)m1(7.2.9)將插值條件(7.2.7)式代入，可推得：(7.2.10)(7.2.9)、(7.2.10)兩式構(gòu)成了三次Hermite插值基本提法二的插值公式?！?.2.2一維數(shù)據(jù)的三次樣條插值上述介紹的分段多項(xiàng)式插值，其優(yōu)點(diǎn)為計(jì)算簡單、穩(wěn)定性好、收斂性有保證，且易于在計(jì)算機(jī)上實(shí)現(xiàn)。但它也明顯存在著缺陷。它只能保證在每個小區(qū)間段[xi,xi+1]內(nèi)光滑，在各小區(qū)間連接點(diǎn)xi處連續(xù)，卻不能保證整條曲線的光滑、光順性，難以滿足某些工程的要求。對于象高速飛機(jī)的機(jī)翼形線，船體放樣等型值線往往要求有二階光滑度，即有二階連續(xù)導(dǎo)數(shù)。而由60年代開始，首先起源及航空、造船業(yè)等工程設(shè)計(jì)的實(shí)際需要而發(fā)展起來的樣條插值，既保留了分段多項(xiàng)式插值的各種優(yōu)點(diǎn)，又提高了插值函數(shù)的光滑度。在此，僅介紹應(yīng)用最廣且具有二階連續(xù)導(dǎo)數(shù)的三次樣條插值方法。三次樣條插值問題的基本提法對于給定的一維數(shù)據(jù)xx0x1…Xny=f(x)y0y1…Yn求一個三次多項(xiàng)式S(x)滿足條件（1）S(xi)=yi，i=0,1,…,n；（2）S(x)具有二階連續(xù)導(dǎo)數(shù)，特別在節(jié)點(diǎn)xi上應(yīng)滿足連續(xù)性要求，即對i=0,1,…,n有三次樣條插值函數(shù)給定區(qū)間[a,b]的一個劃分：a=x0<x1<…<xn=b，設(shè)函數(shù)y=f(x)在節(jié)點(diǎn)xi上的值為yi=f(xi)，i=0,1,…,n。如果S(x)于[a,b]有二階連續(xù)導(dǎo)數(shù)，且在每個小區(qū)間[xi,xi+1]上是三次多項(xiàng)式，則稱S(x)是節(jié)點(diǎn)x0，x1，…，xn上的三次樣條函數(shù)。如果S(x)在節(jié)點(diǎn)xi上還滿足插值條件S(xi)=yi，i=0,1,…,n，(7.2.11)則稱S(x)為三次樣條插值函數(shù)。對應(yīng)于劃分的三次樣條插值函數(shù)的表達(dá)式為(7.2.12)其中，。邊界條件在式(7.2.12)給出的三次多項(xiàng)式中，共含有n3個待定系數(shù)。而由插值條件(7.2.11)式，可列出n1個方程，方程組中未知數(shù)的個數(shù)比方程個數(shù)多2，還需附加2個條件才能進(jìn)行求解。通?？稍趨^(qū)間端點(diǎn)x0=a和xn=b處各附加一個條件（稱為邊界條件或邊值條件）去確定S(x)。邊界條件類型很多，較基本而又常見的有三類：第一邊值條件，即給出邊界點(diǎn)的一階導(dǎo)數(shù)值S(x0)=y0，S(xn)=yn(7.2.13)第二邊值條件，即給出邊界點(diǎn)的二階導(dǎo)數(shù)值S(x0)=y0，S(xn)=yn(7.2.14)特別地，當(dāng)S(x0)=S(xn)=0時，稱為自然邊界條件。滿足自然邊界條件的三次樣條插值函數(shù)稱為自然樣條插值函數(shù)。第三邊值條件（混合邊值條件）(7.2.15)其中1、2、1、2、1、2為定數(shù)。當(dāng)1、2為零時，則為第一邊值條件，當(dāng)1、2為零時，則為第二邊值條件?！?.2.3一維數(shù)據(jù)插值的Matlab實(shí)現(xiàn)1．一維數(shù)據(jù)插值的Matlab實(shí)現(xiàn)Matlab中一維數(shù)據(jù)的插值函數(shù)為：interp1()。其調(diào)用格式為yi=interp1(x,y,xi,'methos'),其中x，y——為插值節(jié)點(diǎn)，均為向量；xi——任取的被插值點(diǎn)可以是一個數(shù)值，也可以是一個向量；yi——為被插值點(diǎn)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)，它們是等價的。2．例：§7.1.1引例1的求解解法一直接借助于Matlab命令，分別采用最臨近插值方法、線性插值方法、三次插值方法和三次樣條插值方法進(jìn)行插值。Matlab程序如下：x=[466,714,950,1422,1634];y=[7.04,4.28,3.40,2.54,2.13];x1=linspace(400,1660,100);y1=interp1(x,y,x1,'nearest');y2=interp1(x,y,x1,'linear');y3=interp1(x,y,x1,'cubic');y4=interp1(x,y,x1,'spline');subplot(2,2,1),plot(x,y,'o',x1,y1),title('最鄰近插值曲線')subplot(2,2,2),plot(x,y,'o',x1,y2),title('線性插值曲線')subplot(2,2,3),plot(x,y,'o',x1,y3),title('三次插值曲線')subplot(2,2,4),plot(x,y,'o',x1,y4),title('三次樣條插值曲線')x2=500:100:1600w1=interp1(x,y,x2,'nearest')w2=interp1(x,y,x2,'linear')w3=interp1(x,y,x2,'cubic')w4=interp1(x,y,x2,'spline')要想求得任一深度值xi處的水溫，只需在上述程序運(yùn)行完畢之后，繼續(xù)運(yùn)行如下命令即可：xi=500:100:1600yi1=interp1(x,y,xi,'nearest')yi2=interp1(x,y,xi,'linear')yi3=interp1(x,y,xi,'cubic')yi4=interp1(x,y,xi,'spline')下面是運(yùn)行上述程序后得到的四種不同插值方法在500米到1600米不同深度處的水溫值（見表7.2.1），以及四種不同的插值函數(shù)的曲線圖（見圖7.2.2）。從圖中可以看出，三次插值和三次樣條插值得到的插值曲線的光滑性更要好些。表7.2.1深度(m)59001000最鄰近插值的水溫(C)7.04004.28004.28004.28003.40003.4000線性插值的水溫(C)6.66165.54874.43583.95933.58643.3089三次插值的水溫(C)6.48875.21444.37033.84333.52033.2905三次樣條插值的水溫(C)6.48875.21444.37033.84333.52033.2905深度(m)110012001600最鄰近插值的水溫(C)3.40002.54002.54002.54002.54002.1300線性插值的水溫(C)3.12672.94452.76232.58012.38922.1958三次插值的水溫(C)3.09252.91422.74602.57792.40032.2032三次樣條插值的水溫(C)3.09252.91422.74602.57792.40032.2032圖7.2.2四種不同插值函數(shù)的曲線圖如果要想求得任一深度值xi處的水溫，只需在上述程序運(yùn)行完畢之后，繼續(xù)運(yùn)行如下命令即可：xi=500:100:1600yi1=interp1(x,y,xi,'nearest')yi2=interp1(x,y,xi,'linear')yi3=interp1(x,y,xi,'cubic')yi4=interp1(x,y,xi,'spline')解法二求三次樣條插值函數(shù)選取第一邊值條件，并用差商和分別代替y0和yn，則所求三次樣條插值函數(shù)為(7.2.16)且滿足（見表7.1.1）：S3(xi)=yi，(i=0,1,…,4)，S3(x0)=0.0111，S3(x4)=0.0019(7.2.17)其中，。將已知數(shù)據(jù)代入(7.2.17)式，得到一個7階線性方程組。解之得（用alpha=linsolve(a,b)命令）：0=3.8854，1=0.0363，2=1.7650104，3=3.2126107，1=5.4789107，2=2.2710107，3=4.5398109，即0=3.8854，1=0.0363，2=0.0002，3=1=2=3=0（用alpha=a\b命令也可得到這個結(jié)果）。并計(jì)算得在500米到1600米不同深度處的水溫值如下：表7.2.2深度(m)59001000水溫(C)6.64925.44114.39523.77463.48823.3157深度(m)110012001600水溫(C)3.14372.96572.77992.58452.37912.1763Matlab程序如下：x=[466,714,950,1422,1634];y=[7.04,4.28,3.40,2.54,2.13];a=[1,x(1),(x(1)^2)/2,(x(1)^3)/6,0,0,0;...1,x(2),(x(2)^2)/2,(x(2)^3)/6,0,0,0;...1,x(3),(x(3)^2)/2,(x(3)^3)/6,((x(3)-x(2))^3)/6,0,0;...1,x(4),(x(4)^2)/2,(x(4)^3)/6,((x(4)-x(2))^3)/6,((x(4)-x(3))^3)/6,0;...1,x(5),(x(5)^2)/2,(x(5)^3)/6,((x(5)-x(2))^3)/6,((x(5)-x(3))^3)/6,((x(5)-x(4))^3)/6;...0,1,x(1),(x(1)^2)/2,0,0,0;...0,1,x(5),(x(5)^2)/2,((x(5)-x(2))^2)/2,((x(5)-x(3))^2)/2,((x(5)-x(4))^2)/2];b=[7.04;4.28;3.40;2.54;2.13;-0.0111;-0.0019];alpha=a\bx1=500:100:1600;fori=1:length(x1)ifx1(i)>=466&x1(i)<=714a1=[1,x1(i),(x1(i)^2)/2,(x1(i)^3)/6,0,0,0];y1(i)=a1*alpha;elseifx1(i)>714&x1(i)<=950a3=[1,x1(i),(x1(i)^2)/2,(x1(i)^3)/6,((x1(i)-x(2))^3)/6,0,0];y1(i)=a3*alpha;elseifx1(i)>950&x1(i)<=1422a4=[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;elsea5=[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ī)律分布的還是散亂的、隨機(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ī)則分布點(diǎn)上的數(shù)據(jù)，轉(zhuǎn)化為數(shù)據(jù)分布有規(guī)律的情形來處理。二維數(shù)據(jù)插值的方法也有很多。在此，針對給定數(shù)據(jù)有規(guī)律分布和散亂分布兩種情形，簡單介紹雙三次樣條插值方法和改進(jìn)的Shepard方法（反距離平方法）的基本概念和基本思想，及其Matlab實(shí)現(xiàn)。§7.3.1雙三次樣條插值雙三次樣條插值方法，是用來解決規(guī)則區(qū)域上給定數(shù)據(jù)有規(guī)律分布的插值問題的常用方法。設(shè)R：[a,b][c,d]是xy平面上的一個矩形區(qū)域。在x軸和y軸上分別取定分割x：a=x0<x1<…<xn<b，y：c=y0<y1<…<ym<d。由此導(dǎo)出R上的一個矩形網(wǎng)格分割：xy，將R分成mn個子矩形Rij：[xi1,xi][yi1,yi]。它的兩條鄰邊長分別是hi=xixi1(i=1,2,…,n)，gj=yjyj1(j=1,2,…,m)。直線x=xi(i=1,2,…,n)和y=yj(j=1,2,…,m)稱為分割的兩族網(wǎng)格線。網(wǎng)格線的交點(diǎn)(xi,yj)(i=1,2,…,n，j=1,2,…,m)稱為分割的節(jié)點(diǎn)，總共有(m1)(n1)個節(jié)點(diǎn)，見圖7.3.1。ydRyiRijyi1c0axi1xib圖7.3.1矩形分割設(shè)函數(shù)S(x,y)定義在矩形域R上且滿足下列條件：（1）在每個子矩形Rij(i=1,2,…,n，j=1,2,…,m)上，S(x,y)是一個關(guān)于x和y的都是三次的多項(xiàng)式函數(shù)，即；(7.3.1)（2）在整個R上，函數(shù)S(x,y)的偏導(dǎo)數(shù)都是連續(xù)的。則稱S(x,y)為雙三次樣條函數(shù)。如果給定一數(shù)組{sij}(i=1,2,…,n，j=1,2,…,m)，雙三次樣條函數(shù)S(x,y)還滿足插值條件S(xi,yi)=sij(i=1,2,…,n，j=1,2,…,m)，就稱S(x,y)為插值雙三次樣條函數(shù)。實(shí)際上，雙三次樣條函數(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è)計(jì)出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個點(diǎ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ù)的點(diǎn)處，定義其函數(shù)值由已知數(shù)據(jù)點(diǎn)及該點(diǎn)距離的近或遠(yuǎn)作加權(quán)平均決定。記，則二元函數(shù)（曲面）定義為，其中(7.3.2)如此定義的曲面是全局相關(guān)的，對曲面的任一點(diǎn)作數(shù)據(jù)計(jì)算都要涉及到全體數(shù)據(jù)，這在大量實(shí)測數(shù)據(jù)插值和擬合中是很慢的。此外，F(xiàn)(x，y)在每個插值點(diǎn)(xk，yk)附近產(chǎn)生一個小的“平臺”，使曲面不具有光順性。當(dāng)fk=c（常數(shù)）時，F(xiàn)(x，y)=c。當(dāng)(xk，yk)皆位于非常數(shù)的同一平面上時，F(xiàn)(x，y)并不能保證再生這個平面。但因?yàn)檫@種作法思想簡單，人們對它作了種種改進(jìn)，例如，取適當(dāng)?shù)某?shù)R>0，令由于()是可微函數(shù)，使得如下定義的F(x，y)在性能上有所改善，其中。(7.3.3)按照上述的思想，可從給定的數(shù)據(jù)恢復(fù)出規(guī)則分布點(diǎn)上的數(shù)據(jù)，接下來就可應(yīng)用雙三次樣條插值或其它的二維數(shù)據(jù)插值方法來處理。§7.3.3二維數(shù)據(jù)插值的Matlab實(shí)現(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é)點(diǎn)，均為向量；zi——為被插值點(diǎn)(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é)點(diǎn)，均為向量；px(j)，py(i)——為被插值節(jié)點(diǎn)；pf(i,j)——為被插值點(diǎn)(px(j),py(i))處的插值結(jié)果；其它輸出參數(shù)涉及插值算法，可以不用了解。e01sef的輸出fnodes和rnw為確定插值的參數(shù)，它們是e01sff需要的輸入?yún)?shù)，因此兩函數(shù)需配合使用。3．例例1氣旋變化情況的可視化表7.3.1是氣象學(xué)家測量得到的氣象資料，它們分別表示在南半球地區(qū)按不同緯度、不同月份的平均氣旋數(shù)字。根據(jù)這些數(shù)據(jù)，繪制出氣旋分布曲面圖形。表7.3.10—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.3.2四種插值方法的可視化圖形最鄰近插值的Matlab程序?yàn)椋簒=1:12;y=5:10:85;z=[2.4,18.7,20.8,22.1,37.3,48.2,25.6,5.3,0.31.6,21.4,18.5,20.1,28.8,36.6,24.2,5.3,02.4,16.2,18.2,20.5,27.8,35.5,25.5,5.4,03.2,9.2,16.6,25.1,37.2,40,24.6,4.9,0.31.0,2.8,12.9,29.2,40.3,37.6,21.1,4.9,00.5,1.7,10.1,32.6,41.7,35.4,22.2,7.1,00.4,1.4,8.3,33.0,46.2,35,20.2,5.3,0.10.2,2.4,11.2,31.0,39.9,34.7,21.2,7.3,0.20.5,5.8,12.5,28.6,25.9,35.7,22.6,7,0.30.8,9.2,21.1,32.0,40.3,39.5,28.5,8.6,02.4,10.3,23.9,28.1,38.2,40,25.3,6.3,0.13.6,16,25.5,25.6,43.4,41.9,24.3,6.6,0.3];[xi,yi]=meshgrid(1:12,5:1:85);zi=interp2(x,y,z,xi,yi,'nearest');figuresurf(xi,yi,zi)xlabel('月份'),ylabel('緯度'),zlabel('氣旋'),axis([012090050])title('南半球氣旋可視化圖形')雙線性插值、雙三次插值、雙三次樣條插值的Matlab程序?yàn)椋悍謩e將最鄰近線性插值程序中的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=interp2(x,y,z,xi,yi,'spline')例2水道測量數(shù)據(jù)（AMCM86A—題）在某海域測得一些點(diǎn)（x,y）處的水深z（單位：英尺）由表7.3.2給出，水深數(shù)據(jù)是在低潮時測得的。船的吃水深度為5英尺，問在矩形區(qū)域(75，200)(50，150)里的哪些地方船要避免進(jìn)入。表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.5y6.581.03.056.566.584.038.5z9988949解（1）假設(shè)由題目給出的信息是很少的，除了14個位置的水深之外一無所知。顯然，題目要求我們找出水深不到5英尺的區(qū)域。為了討論方便，下面三個假設(shè)是合理的：=1\*GB3①所給數(shù)據(jù)是精確的；=2\*GB3②討論區(qū)域的海底曲面是光滑的，更確切地說，可以認(rèn)為曲面的一階、二階導(dǎo)數(shù)是連續(xù)的。因?yàn)槲覀兛梢哉J(rèn)為討論區(qū)域?yàn)闇\水海域，由于長期的海水水流作用，形成的是以礫石或沙為主要組成部分的海底，不存在珊瑚礁、水底峽谷、山脊等不可意料的突變地形。=3\*GB3③水深是一個按區(qū)域來劃分的變量，在某個位置的水深及其周圍區(qū)域的水深是相互依賴的，但這種依賴作用隨距離的增大而減小。就我們討論的問題來說，每一個給定數(shù)據(jù)點(diǎn)影響周圍的每一個未知點(diǎn)，一個給定數(shù)據(jù)點(diǎn)離未知點(diǎn)越近，作用就越大。（2）問題分析根據(jù)假設(shè)，海底曲面是連續(xù)光滑的，不存在珊瑚礁、水底峽谷、山脊等不可意料的突變地形，因而很自然的想法就是用某種光滑的擬合曲面去逼近已知的14個數(shù)據(jù)點(diǎn)或以14個已知的數(shù)據(jù)點(diǎn)為基礎(chǔ)，利用二維插值補(bǔ)充一些點(diǎn)的水深，以求得水深不超過5米的區(qū)域。在此，我們采用二維插值方法，應(yīng)用Matlab程序，作出矩形區(qū)域(75，200)(50，150)范圍內(nèi)的海底地形圖、水深不超過5米的危險區(qū)域的平面圖以及水深不超過5米的危險區(qū)域的海底地貌圖，并求出水深不超過5米的危險海域范圍。（3）問題求解采用改進(jìn)的Shepard方法，利用Matlab軟件作出作出矩形區(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,'+'),title('測量點(diǎn)分布圖');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);fori=1:nyforj=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)范圍內(nèi)的海底地形圖');[a,b]=find(pf>=-5);amin=min(a);amax=max(a);bmin=min(b);bmax=max(b);xmin=75+((200-75)/100)*bminxmax=75+((200-75)/100)*bmaxymin=-50+((150+50)/200)*aminymax=-50+((150+50)/200)*amaxfork=1:length(b)i0(k)=75+((200-75)/100)*b(k);endfork=1:length(a)j0(k)=-50+((150+50)/200)*a(k);endsubplot(2,2,3),plot(i0,j0,'+'),title('水深不超過5米的危險區(qū)域的平面圖');[i1,j1]=find(pf<-5);fork=1:length(i1)pf(i1(k),j1(k))=-5;endsubplot(2,2,4),meshc(px,py,pf),title('水深不超過5米的危險區(qū)域的海底地貌圖');rotate3d§7.4一維數(shù)據(jù)擬合的最小二乘法簡介數(shù)據(jù)擬合問題的形式多種多樣，解決的方法也有許多。在此，我們只簡單介紹以最小二乘法為準(zhǔn)則的一維數(shù)據(jù)擬合方法?！?.4.1數(shù)據(jù)擬合的最小二乘法簡介對于給定的一組測量數(shù)據(jù)xx1x2…xny=f(x)y1y2…yn設(shè)y=(x)為一擬合函數(shù)，記i=(xi)yi(i=1,2,…,n)，(7.4.1)則稱i為擬合函數(shù)(x)在xi點(diǎn)處的偏差或殘量。為使(x)在整體上盡可能及給定數(shù)據(jù)的函數(shù)f(x)近似，我們常采用偏差的平方和達(dá)到最小，即(7.4.2)來保證每個偏差的絕對值|i|都很小。這一原則稱為最小二乘原則，根據(jù)最小二乘原則確定擬合函數(shù)(x)的方法稱為最小二乘法。線性最小二乘擬合我們知道，函數(shù)系{xk|k=0,1,…,m}的線性組合(x)=a0a1xa2x2…amxm為m次多項(xiàng)式。一般地，若函數(shù)系{k(x)|k=0,1,…,m}是線性無關(guān)的，則其線性組合(7.4.3)稱為函數(shù)系{k(x)|k=0,1,…,m}的廣義多項(xiàng)式。如三角多項(xiàng)式就是函數(shù)系{1,cosx,sinx,cos2x,sin2x,…,cosmx,sinmx}的廣義多項(xiàng)式。設(shè){k(x)|k=0,1,…,m}為一線性無關(guān)的函數(shù)系，取擬合函數(shù)為(7.4.3)式給出的廣義多項(xiàng)式，使得(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）即亦即，（k=0,1,…,m）(7.4.4)式(7.4.4)是關(guān)于a0,a1,a2,…,am的線性方程組，稱為正規(guī)方程組。從正規(guī)方程組(7.4.4)中解出a0,a1,a2,…,am，于是就得到了最小二乘擬合函數(shù)(x)。非線性最小二乘擬合如果擬合函數(shù)(x)=(x,a0,a1,a2,…,am)的待定參數(shù)a0,a1,a2,…,am不能全部以線性形式出現(xiàn)，如指數(shù)擬合函數(shù)等，這便是非線性最小二乘擬合問題。一般地，非線性最小二乘擬合問題是一個非線性函數(shù)的極小化問題，可用非線性優(yōu)化方法求解。最小二乘擬合函數(shù)的選擇最小二乘法中，擬合函數(shù)的選擇是很重要的。可以通過對給定數(shù)據(jù)的分析來選擇，也可以直接由實(shí)際問題給定。最常用的是多項(xiàng)式和樣條函數(shù)，尤其是當(dāng)不知道該選擇什么樣的擬合函數(shù)時，通常可以考慮選擇樣條函數(shù)。另外，對同一問題，也可選擇不同的函數(shù)進(jìn)行最小二乘擬合，比較各自誤差的大小，從中選出誤差較小的作為擬合函數(shù)?！?.4.2一維數(shù)據(jù)擬合的Matlab實(shí)現(xiàn)多項(xiàng)式函數(shù)擬合擬合函數(shù)的命令為：polyfit()，其調(diào)用格式為a=polyfit(xdata,ydata,m),其中m——表示多項(xiàng)式的最高階數(shù)；xdata，ydata——為將要擬合的數(shù)據(jù)，它們都是以數(shù)組方式輸入；a——輸出參數(shù)，為擬合多項(xiàng)式的系數(shù)a=[a0,a1,a2,…,am]。多項(xiàng)式在x處的值y可用如下命令格式計(jì)算：y=polyval(a,x)。一般的曲線擬合擬合函數(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ù)的初值。若要求點(diǎn)x處的函數(shù)值y，可用程序f=Fun(p,x)計(jì)算。例：§7.1.2引例2的求解解：（1）選擇擬合曲線作出給定數(shù)據(jù)的散點(diǎn)圖如下：圖7.4.1給定數(shù)據(jù)的散點(diǎn)圖通過對散點(diǎn)圖的分析可以看出，數(shù)據(jù)點(diǎn)的分布為一條單調(diào)上升的曲線。具有這種特性的曲線很多，我們選取如下三種函數(shù)來擬合：=1\*GB3①多項(xiàng)式(x)=a0a1x…amxm，m為適當(dāng)選取的正整數(shù)；=2\*GB3②；=3\*GB3③（a>0,b<0）。（2）擬合運(yùn)算首先，分別用二、三、六次多項(xiàng)式擬合，計(jì)算得輸出參數(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.0711x4.3252x2(2)(x)=0.00600.1963x2.1346x22.5952x3(3)(x)=0.0004x0.0103x20.1449x31.1395x44.9304x50.0498x6；其次，再用有理分式擬合，計(jì)算得輸出參數(shù)分別為p=[0.0841，0.1392]擬合函數(shù)為最后，用指數(shù)函數(shù)擬合，計(jì)算得輸出參數(shù)分別為p=[11.3578，1.0873]擬合函數(shù)為三種方式五個種函數(shù)的擬合曲線見圖7.4.2—7.4.4。圖7.4.2多項(xiàng)式函數(shù)擬合曲線圖圖7.4.3有理分式函數(shù)擬合曲線圖圖7.4.4指數(shù)函數(shù)擬合曲線圖（3）誤差分析和給定的16組數(shù)據(jù)比較，三種方式五個函數(shù)擬合的誤差見下表：表7.4.1五個函數(shù)擬合的誤差表偏差平方和平均偏差平方和最大偏差二次多項(xiàng)式擬合4.44160.27761.3518三次多項(xiàng)式擬合1.22920.07680.6067六次多項(xiàng)式擬合0.11690.00730.2684有理分式擬合0.57320.03580.4772指數(shù)函數(shù)擬合0.17770.01110.2544其中：偏差平方和為，平均偏差平方和為。從上表中求得的誤差情況來看，好似六次多項(xiàng)式函數(shù)擬合的最好，指數(shù)函數(shù)擬合次之，然后分別是有理分式函數(shù)擬合、三次多項(xiàng)式函數(shù)擬合和二次多項(xiàng)式函數(shù)擬合。但是，就這個實(shí)際問題的本質(zhì)來說，化學(xué)反應(yīng)中生成物的濃度到一定時間后應(yīng)基本穩(wěn)定，即當(dāng)t時，f(t)常數(shù)。而我們有并且三個多項(xiàng)式函數(shù)擬合曲線的趨勢也可從圖7.4.5看出。因而，我們有理由認(rèn)為指數(shù)函數(shù)擬合的效果最好，有理分式函數(shù)擬合的效果次之，本問題不宜使用多項(xiàng)式擬合。（4）結(jié)論通過以上的計(jì)算和分析，我們得出如下結(jié)論：本問題可用指數(shù)函數(shù)或有理分式函數(shù)來擬合，其擬合函數(shù)分別為和擬合誤差為表7.4.2偏差平方和平均偏差平方和最大偏差指數(shù)函數(shù)擬合0.17770.01110.4772有理分式擬合0.57320.03580.2544圖7.4.5（5）Matlab程序=1\*GB3①多項(xiàng)式擬合的Matlab程序：t=1:16;y=[4,6.4,8,8.4,9.28,9.5,9.7,9.86,10,10.2,10.32,10.42,10.5,10.55,10.58,10.6];figure(1)%作出給定數(shù)據(jù)的散點(diǎn)圖plot(t,y,'o')p1=polyfit(t,y,2)p2=polyfit(t,y,3)p3=polyfit(t,y,6)x1=linspace(0,16,160);y1=polyval(p1,x1);y2=polyval(p2,x1);y3=polyval(p3,x1);figure(2)%作出三條擬合曲線圖plot(t,y,'.',x1,y1,'-',x1,y2,':',x1,y3,'-.')u1=polyval(p1,t);%以下計(jì)算誤差u2=polyval(p2,t);u3=polyval(p3,t);d1=(u1-y).^2;d2=(u2-y).^2;d3=(u3-y).^2;delta11=sum(d1)delta21=sum(d2)delta31=sum(d3)delta12=sum(d1)/16delta22=sum(d2)/16delta32=sum(d3)/16delta13=sqrt(abs(max(d1)))delta23=sqrt(abs(max(d2)))delta33=sqrt(abs(max(d3)))x2=linspace(0,25,60);%畫三條擬合曲線的趨勢圖w1=polyval(p1,x2);w2=polyval(p2,x2);w3=polyval(p3,x2);figure(3)plot(t,y,'.',x2,w1,'-',x2,w2,':',x2,w3,'-.')title('三條擬合曲線的趨勢圖')=2\*GB3②有理分式擬合的Matlab程序：tdata=1:16;ydata=[4,6.4,8,8.4,9.28,9.5,9.7,9.86,10,10.2,10.32,10.42,10.5,10.55,10.58,10.6];figure(1)plot(tdata,ydata,'o')p0=[0.1,0.1];cde:\matlab\model\%input'p=curvefit('e111fun',p0,tdata,ydata)fori=1:16phi(i)=p(1)*exp(p(2)/tdata(i));endfigure(2)plot(tdata,phi)figure(3)plot(tdata,ydata,'.',tdata,phi)d=(phi-ydata).^2;delta1=sum(d)delta2=sum(d)/16delta3=abs(max(d))=3\*GB3③指數(shù)函數(shù)擬合的Matlab程序：tdata=1:16;ydata=[4,6.4,8,8.4,9.28,9.5,9.7,9.86,10,10.2,10.32,10.42,10.5,10.55,10.58,10.6];figure(1)plot(tdata,ydata,'o')p0=[0.1,0.1];cde:\matlab\model\%input'p=curvefit('e110fun',p0,tdata,ydata)fori=1:16phi(i)=tdata(i)/(p(1)*tdata(i)+p(2));endfigure(2)plot(tdata,phi)figure(3)plot(tdata,ydata,'.',tdata,phi)title('有理分式擬合')d=(phi-ydata).^2;delta1=sum(d)delta2=sum(d)/16delta3=abs(max(d))=4\*GB3④調(diào)用的M文件程序：e110fun.m文件functionphi=e90fun(p,tdata)fori=1:16phi(i)=tdata(i)/(p(1)*tdata(i)+p(2));endphi;e111fun.m文件functionphi=e91fun(p,tdata)fori=1:16phi(i)=p(1)*exp(p(2)/tdata(i));endphi;參考文獻(xiàn)：易大義等，數(shù)值方法，浙江科學(xué)技術(shù)出版社，1984。徐濤，數(shù)值計(jì)算方法，吉林科學(xué)技術(shù)出版社，2019。壽紀(jì)麟，數(shù)學(xué)建?！椒胺独?，西安交通大學(xué)出版社，1993。葉其孝主編，大學(xué)生數(shù)學(xué)建模競賽輔導(dǎo)教材，湖南教育出版社，1993。李尚志主編，數(shù)學(xué)建模競賽教程，江蘇教育出版社，1996。傅鵬等編著，數(shù)學(xué)實(shí)驗(yàn)，科學(xué)出版社，2000。謝云蓀主編，數(shù)學(xué)實(shí)驗(yàn)，科學(xué)出版社，2000。應(yīng)用Matlab求解水道測量數(shù)據(jù)問題摘要：水道測量數(shù)據(jù)問題是一個給定數(shù)據(jù)散亂、隨機(jī)分布的二維離散數(shù)據(jù)的插值問題。本文以水道測量數(shù)據(jù)問題為例，應(yīng)用Matlab軟件提供的求解三維網(wǎng)格點(diǎn)數(shù)據(jù)的函數(shù)，對求解決給定數(shù)據(jù)散亂、隨機(jī)分布的二維離散數(shù)據(jù)插值問題，給出了一個簡便易行的方法。問題重述水道測量數(shù)據(jù)問題是1986年美國大學(xué)生數(shù)學(xué)建模競賽的A題，由加州海軍研究生院數(shù)學(xué)系的RichardFranke提供。問題如下：在某海域測得一些點(diǎn)（x,y）處的水深z（單位：英尺）由表1給出，水深數(shù)據(jù)是在低潮時測得的。船的吃水深度為5英尺，問在矩形區(qū)域(75，200)(50，150)里的哪些地方船要避免進(jìn)入。表1水道水深測量數(shù)據(jù)（單位：英尺）x129.0140.0108.588.0185.5195.0105.5y7.5141.528.0147.022.5137.585.5z4868688x157.5107.577.0