1、數(shù)值計(jì)算3插值和曲線擬合插值法是實(shí)用的數(shù)值方法，是函數(shù)逼近的重要方法。在生產(chǎn)和科學(xué)實(shí)驗(yàn)中，自變量x與因變量y的函數(shù)y = f(x)的關(guān)系式有時(shí)不能直接寫出表達(dá)式，而只能得到函數(shù)在若干個(gè)點(diǎn)的函數(shù)值或?qū)?shù)值。當(dāng)要求知道觀測(cè)點(diǎn)之外的函數(shù)值時(shí)，需要估計(jì)函數(shù)值在該點(diǎn)的值。如何根據(jù)觀測(cè)點(diǎn)的值，構(gòu)造一個(gè)比較簡(jiǎn)單的函數(shù)y=(x)，使函數(shù)在觀測(cè)點(diǎn)的值等于已知的數(shù)值或?qū)?shù)值。用簡(jiǎn)單函數(shù)y=(x)在點(diǎn)x處的值來(lái)估計(jì)未知函數(shù)y=f(x)在x點(diǎn)的值。尋找這樣的函數(shù)(x)，辦法是很多的。(x)可以是一個(gè)代數(shù)多項(xiàng)式，或是三角多項(xiàng)式，也可以是有理分式；(x)可以是任意光滑（任意階導(dǎo)數(shù)連續(xù)）的函數(shù)或是分段函數(shù)。函數(shù)類的不同，自

2、然地有不同的逼近效果。在許多應(yīng)用中，通常要用一個(gè)解析函數(shù)（一、二元函數(shù)）來(lái)描述觀測(cè)數(shù)據(jù)。根據(jù)測(cè)量數(shù)據(jù)的類型：1測(cè)量值是準(zhǔn)確的，沒(méi)有誤差。2測(cè)量值與真實(shí)值有誤差。這時(shí)對(duì)應(yīng)地有兩種處理觀測(cè)數(shù)據(jù)方法：1插值或曲線擬合。2回歸分析（假定數(shù)據(jù)測(cè)量是精確時(shí)，一般用插值法，否則用曲線擬合）。MATLAB中提供了眾多的數(shù)據(jù)處理命令。有插值命令，有擬合命令，有查表命令。一維插值 插值定義為對(duì)數(shù)據(jù)點(diǎn)之間函數(shù)的估值方法，這些數(shù)據(jù)點(diǎn)是由某些集合給定。當(dāng)人們不能很快地求出所需中間點(diǎn)的函數(shù)值時(shí)，插值是一個(gè)有價(jià)值的工具。例如，當(dāng)數(shù)據(jù)點(diǎn)是某些實(shí)驗(yàn)測(cè)量的結(jié)果或是過(guò)長(zhǎng)的計(jì)算過(guò)程時(shí)，就有這種情況。 interp1(x,y,xi,m

3、ethod) x和y為既有數(shù)據(jù)的向量，其長(zhǎng)度必須相同。 xi為要插值的數(shù)據(jù)點(diǎn)向量。 method插值方法,nearest/linear/cubic/spline之一，分別為最近點(diǎn)插值/線性插值/分段三次Hermite插值/三次樣條插值。 例x=1.0 2.0 3.0 4.0 5.0; %輸入變量數(shù)據(jù)x y=11.2 16.5 20.4 26.3 30.5; %輸入變量數(shù)據(jù)y x1=2.55; %輸入待插值點(diǎn)x y11=interp1(x,y,x1,nearest) %最近點(diǎn)插值方法的插值結(jié)果y12=interp1(x,y,x1,linear) %線性插值方法的插值結(jié)果y13=interp1(x

4、,y,x1,cubic) %三次Hermite插值方法的插值結(jié)果y14=interp1(x,y,x1,spline) %樣條插值方法的插值結(jié)果 y11 = 20.4000y12 = 18.6450y13 = 18.6028y14 = 18.4874 plot(x,y) 或許最簡(jiǎn)單插值的例子是MATLAB的作圖。按缺省，MATLAB用直線連接所用的數(shù)據(jù)點(diǎn)以作圖。這個(gè)線性插值猜測(cè)中間值落在數(shù)據(jù)點(diǎn)之間的直線上。當(dāng)然，當(dāng)數(shù)據(jù)點(diǎn)個(gè)數(shù)的增加和它們之間距離的減小時(shí)，線性插值就更精確。例如： x1=linspace(0, 2*pi, 60);x2=linspace(0, 2*pi, 6);plot(x1, s

5、in(x1), x2, sin(x2), - )xlabel( x ), ylabel( sin(x) ), title( Linear Interpolation ) 上圖是sin函數(shù)的兩個(gè)圖，一個(gè)在數(shù)據(jù)點(diǎn)之間用60個(gè)點(diǎn)，它比另一個(gè)只用6個(gè)點(diǎn)更光滑和更精確。 根據(jù)所作的假設(shè)，有多種插值。而且，可以在一維以上空間中進(jìn)行插值。即如果有反映兩個(gè)變量函數(shù)的插值，z=f(x, y)，那么就可在x之間和在y之間，找出z的中間值進(jìn)行插值。MATLAB在一維函數(shù)interp1和在二維函數(shù)interp2中，提供了許多的插值選擇。其中的每個(gè)函數(shù)將在下面闡述。 為了說(shuō)明一維插值，考慮下列問(wèn)題，12小時(shí)內(nèi)，一小時(shí)測(cè)

6、量一次室外溫度。數(shù)據(jù)存儲(chǔ)在兩個(gè)MATLAB變量中。hours=1:12; % index for hour data was recordedtemps=5 8 9 15 25 29 31 30 22 25 27 24; % recorded temperaturesplot(hours, temps, hours, temps, + ) % view temperaturestitle( Temperature )xlabel( Hour ), ylabel( Degrees Celsius ) 在線性插值下室外溫度曲線正如上圖看到的，MATLAB畫出了數(shù)據(jù)點(diǎn)線性插值的直線。為了計(jì)算在任意給

7、定時(shí)間的溫度，人們可試著對(duì)可視的圖作解釋。另外一種方法，可用函數(shù)interp1。t=interp1(hours, temps, 9.3) % estimate temperature at hour=9.3 t = 22.9000 t=interp1(hours, temps, 4.7) % estimate temperature at hour=4.7 t = 22 t=interp1(hours, temps, 3.2 6.5 7.1 11.7) % find temp at many points! t = 10.2000 30.0000 30.9000 24.9000 interp1

8、的缺省用法是由interp1(x, y, xo)來(lái)描述，這里x是獨(dú)立變量(橫坐標(biāo))，y是應(yīng)變量(縱坐標(biāo))，xo是進(jìn)行插值的一個(gè)數(shù)值數(shù)組。另外，該缺省的使用假定為線性插值。若不采用直線連接數(shù)據(jù)點(diǎn)，我們可采用某些更光滑的曲線來(lái)擬合數(shù)據(jù)點(diǎn)。最常用的方法是用一個(gè)3階多項(xiàng)式，即3次多項(xiàng)式，來(lái)對(duì)相繼數(shù)據(jù)點(diǎn)之間的各段建模，每個(gè)3次多項(xiàng)式的頭兩個(gè)導(dǎo)數(shù)與該數(shù)據(jù)點(diǎn)相一致。這種類型的插值被稱為3次樣條或簡(jiǎn)稱為樣條。函數(shù)interp1也能執(zhí)行3次樣條插值。 hours=1:12; temps=5 8 9 15 25 29 31 30 22 25 27 24; plot(hours, temps, hours, tem

9、ps, + ) title( Temperature )xlabel(時(shí)間軸), ylabel(溫度變化軸) x1=9.3 2.5;t1=interp1(hours,temps,x1,spline) t1 = 21.8577 8.2364 x2=3.2 6.5 7.1 11.7;t3=interp1(hours,temps,x2,spline) t3 = 9.6734 30.0427 31.1755 25.3820 t3 = 9.6734 30.0427 31.1755 25.3820 注意，樣條插值得到的結(jié)果，與上面所示的線性插值的結(jié)果不同。因?yàn)椴逯凳且粋€(gè)估計(jì)或猜測(cè)的過(guò)程，其意義在于，應(yīng)用不

10、同的估計(jì)規(guī)則導(dǎo)致不同的結(jié)果。 一個(gè)最常用的樣條插值是對(duì)數(shù)據(jù)平滑。也就是，給定一組數(shù)據(jù)，使用樣條插值在更細(xì)的間隔求值。例如， hours=1:12; temps=5 8 9 15 25 29 31 30 22 25 27 24;h=1:0.1:12; % estimate temperature every 1/10 hourt=interp1(hours,temps,h,spline) ;plot(hours, temps, - , hours, temps, + , h, t) % plot comparative resultstitle( Springfield Temperature

11、)xlabel( Hour ), ylabel( Degrees Celsius ) 在不同插值下室外溫度曲線 在上圖中，虛線是線性插值，實(shí)線是平滑的樣條插值，標(biāo)有 + 的是原始數(shù)據(jù)。如要求在時(shí)間軸上有更細(xì)的分辨率，并使用樣條插值，我們有一個(gè)更平滑、但不一定更精確地對(duì)溫度的估計(jì)。尤其應(yīng)注意，在數(shù)據(jù)點(diǎn)，樣條解的斜率不突然改變。作為這個(gè)平滑插值的回報(bào)，3次樣條插值要求更大量的計(jì)算，因?yàn)楸仨氄业?次多項(xiàng)式以描述給定數(shù)據(jù)之間的特征。 在討論二維插值之前，了解interp1所強(qiáng)制的二個(gè)強(qiáng)約束是很重要的。首先，人們不能要求有獨(dú)立變量范圍以外的結(jié)果，例如，interp1(hours, temps, 13.5

12、)導(dǎo)致一個(gè)錯(cuò)誤，因?yàn)閔ours在1到12之間變化。其次，獨(dú)立變量必須是單調(diào)的。即獨(dú)立變量在值上必須總是增加的或總是減小的。在我們的例子里，hours是單調(diào)的。然而，如果我們已經(jīng)定義獨(dú)立變量為一天的實(shí)際時(shí)間， time_of_day=7:12 1:6 % start at 7AM,end at 6PM time_of_day = 7 8 9 10 11 12 1 2 3 4 5 6則獨(dú)立變量將不是單調(diào)的，因?yàn)閠ime_of_day增加到12，然后跌到1，再然后增加。如果用time_of_day代替interp1中的hours，將會(huì)返回一個(gè)錯(cuò)誤。同樣的理由，人們不能對(duì)temps插值來(lái)找出產(chǎn)生某溫度

13、的時(shí)間(小時(shí))，因?yàn)閠emps不是單調(diào)的。二維插值interp2(x,y,z,xi,yi,method) x,y,z為既有數(shù)據(jù)的向量，其長(zhǎng)度必須相同。 xi,yi為要插值的數(shù)據(jù)點(diǎn)。 method為插值方法，可為nearest、linear、cubic之一， 分別為最近點(diǎn)插值、雙線性插值、雙三次插值%例：x=0.0 0.5 1.0 1.5 2.0 2.5 3.0; %輸入變量數(shù)據(jù)x y=0.0 0.5 1.0 1.5 2.0 2.5 3.0; %輸入變量數(shù)據(jù)y %輸入變量數(shù)據(jù)z z=100 99 100 99 100 99 100 100 99 99 99 100 99 100 99 99 98

14、 98 100 99 100 100 98 97 97 99 100 103 101 100 98 98 100 102 104 102 103 101 100 102 106 106 99 102 100 100 103 108 103; x1=2.3;y1=2.8; %輸入待插值點(diǎn)x1,y1 z11=interp2(x,y,z,x1,y1,nearest) %最近點(diǎn)插值方法的插值結(jié)果z12=interp2(x,y,z,x1,y1,bilinear) %雙線性插值插值方法的插值結(jié)果z13=interp2(x,y,z,x1,y1,bicubic) %雙三次插值插值方法的插值結(jié)果 z11 = 1

15、08z12 = 105.3600z13 = 105.9744 二維插值是基于與一維插值同樣的基本思想。然而，正如名字所隱含的，二維插值是對(duì)兩變量的函數(shù)z=f(x, y) 進(jìn)行插值。為了說(shuō)明這個(gè)附加的維數(shù)，考慮一個(gè)問(wèn)題。設(shè)人們對(duì)平板上的溫度分布估計(jì)感興趣，給定的溫度值取自平板表面均勻分布的格柵。 采集了下列的數(shù)據(jù)：width=1:5; % index for width of plate (i.e.,the x-dimension)depth=1:3; % index for depth of plate (i,e,the y-dimension)temps=82 81 80 82 84; 79

16、 63 61 65 81; 84 84 82 85 86 % temperature data temps = 82 81 80 82 84 79 63 61 65 81 84 84 82 85 86 如同在標(biāo)引點(diǎn)上測(cè)量一樣，矩陣temps表示整個(gè)平板的溫度分布。temps的列與下標(biāo)depth或y-維相聯(lián)系，行與下標(biāo)width或x-維相聯(lián)系(見(jiàn)圖)。為了估計(jì)在中間點(diǎn)的溫度，我們必須對(duì)它們進(jìn)行辨識(shí)。wi=1:0.2:5; % estimate across width of plated=2; % at a depth of 2zlinear=interp2(width, depth, temp

17、s, wi, d) ; %linear interpolationzcubic=interp2(width, depth, temps, wi,d, cubic) ; % cubic interpolationplot(wi, zlinear, - , wi, zcubic,r) % plot resultsxlabel( Width of Plate ), ylabel( Degrees Celsius )title( Temperature at Depth = num2str(d) ) 在深度d=2處的平板溫度 另一種方法，我們可以在兩個(gè)方向插值。先在三維坐標(biāo)畫出原始數(shù)據(jù)，看一下該數(shù)據(jù)的

18、粗糙程度(見(jiàn)圖11.7)。mesh(width, depth, temps) % use mesh plotxlabel( Width of Plate ), ylabel( Depth of Plate )zlabel( Degrees Celsius ), axis(ij), grid 平板溫度 然后在兩個(gè)方向上插值，以平滑數(shù)據(jù)。width=1:5; depth=1:3; temps=82 81 80 82 84; 79 63 61 65 81; 84 84 82 85 86;di=1:0.2:3; % choose higher resolution for depthwi=1:0.2:

19、5; % choose higher resolution for widthzcubic=interp2(width,depth,temps,wi,di,cubic) ; % cubicmesh(wi, di, zcubic)xlabel( Width of Plate ), ylabel( Depth of Plate )zlabel( Degrees Celsius ),axis(ij),grid 上面的例子清楚地證明了，二維插值更為復(fù)雜，只是因?yàn)橛懈嗟牧恳３指?。interp2的基本形式是interp2(x, y, z, xi, yi, method)。這里x和y是兩個(gè)獨(dú)立變量，z

20、是一個(gè)應(yīng)變量矩陣。x和y對(duì)z的關(guān)系是z(i, :) = f(x, y(i) 和 z(:, j) = f(x(j), y). 也就是，當(dāng)x變化時(shí)，z的第i行與y的第i個(gè)元素y(i)相關(guān)，當(dāng)y變化時(shí)，z的第j列與x的第j個(gè)元素x(j)相關(guān)，。xi是沿x-軸插值的一個(gè)數(shù)值數(shù)組；yi是沿y-軸插值的一個(gè)數(shù)值數(shù)組。圖11.8 二維插值后的平板溫度 可選的參數(shù)method可以是 linear，cubic或nearest。在這種情況下，cubic不意味著3次樣條，而是使用3次多項(xiàng)式的另一種算法。linear方法是線性插值，僅用作連接圖上數(shù)據(jù)點(diǎn)。nearest方法只選擇最接近各估計(jì)點(diǎn)的粗略數(shù)據(jù)點(diǎn)。在所有的情況

21、下，假定獨(dú)立變量x和y是線性間隔和單調(diào)的。關(guān)于這些方法的更多的信息，可請(qǐng)求在線幫助。二曲線擬合使用多項(xiàng)式曲線擬合 p=polyfit(x,y,n) p為最小二乘意義上擬合多項(xiàng)式的相關(guān)系數(shù)。 x,y為數(shù)據(jù)點(diǎn)向量。n為多項(xiàng)式階次。 y=polyval(p,x) y為自變量為x、系數(shù)為p的多項(xiàng)式值. x為矩陣或向量時(shí)對(duì)其每一點(diǎn)進(jìn)行操作. p為一向量，各元素為多項(xiàng)式系數(shù)(降冪排列). 曲線擬合涉及回答兩個(gè)基本問(wèn)題：最佳擬合意味著什么？應(yīng)該用什么樣的曲線？可用許多不同的方法定義最佳擬合，并存在無(wú)窮數(shù)目的曲線。所以，從這里開始，我們走向何方？正如它證實(shí)的那樣，當(dāng)最佳擬合被解釋為在數(shù)據(jù)點(diǎn)的最小誤差平方和，且

22、所用的曲線限定為多項(xiàng)式時(shí)，那么曲線擬合是相當(dāng)簡(jiǎn)捷的。數(shù)學(xué)上，稱為多項(xiàng)式的最小二乘曲線擬合。如果這種描述使你混淆，再研究圖。虛線和標(biāo)志的數(shù)據(jù)點(diǎn)之間的垂直距離是在該點(diǎn)的誤差。對(duì)各數(shù)據(jù)點(diǎn)距離求平方，并把平方距離全加起來(lái)，就是誤差平方和。這條虛線是使誤差平方和盡可能小的曲線，即是最佳擬合。最小二乘這個(gè)術(shù)語(yǔ)僅僅是使誤差平方和最小的省略說(shuō)法。在MATLAB中，函數(shù)polyfit求解最小二乘曲線擬合問(wèn)題。為了闡述這個(gè)函數(shù)的用法，讓我們以上圖中的數(shù)據(jù)開始。x=0 .1 .2 .3 .4 .5 .6 .7 .8 .9 1; y=-.447 1.978 3.28 6.16 7.08 7.34 7.66 9.56

23、9.48 9.30 11.2; 為了用polyfit，我們必須給函數(shù)賦予上面的數(shù)據(jù)和我們希望最佳擬合數(shù)據(jù)的多項(xiàng)式的階次或度。如果我們選擇n=1作為階次，得到最簡(jiǎn)單的線性近似。通常稱為線性回歸。相反，如果我們選擇n=2作為階次，得到一個(gè)2階多項(xiàng)式。現(xiàn)在，我們選擇一個(gè)2階多項(xiàng)式。n=2; % polynomial orderp=polyfit(x, y, n) p = -9.8108 20.1293 -0.0317 polyfit 的輸出是一個(gè)多項(xiàng)式系數(shù)的行向量。其解是y = 9.8108x2 20.1293x0.0317。為了將曲線擬合解與數(shù)據(jù)點(diǎn)比較，讓我們把二者都繪成圖。xi=linspace

24、(0, 1, 100); % x-axis data for plotting z=polyval(p, xi); 為了計(jì)算在xi數(shù)據(jù)點(diǎn)的多項(xiàng)式值，調(diào)用MATLAB的函數(shù)polyval。 plot(x, y, o , x, y, xi, z, : ) 畫出了原始數(shù)據(jù)x和y，用o標(biāo)出該數(shù)據(jù)點(diǎn)，在數(shù)據(jù)點(diǎn)之間，再用直線重畫原始數(shù)據(jù)，并用點(diǎn) : 線，畫出多項(xiàng)式數(shù)據(jù)xi和z。 xlabel( x ), ylabel( y=f(x) ), title( Second Order Curve Fitting ) 將圖作標(biāo)志。這些步驟的結(jié)果表示于前面的圖中。例設(shè)0.00.30.81.11.62.3，y=0.5

25、00.821.141.251.351.40，試求二次多項(xiàng)式擬合系數(shù)，并據(jù)此計(jì)算x1=0.9 1.2時(shí)對(duì)應(yīng)的y1。x=0.0 0.3 0.8 1.1 1.6 2.3; %輸入變量數(shù)據(jù)xy=0.50 0.82 1.14 1.25 1.35 1.40; %輸入變量數(shù)據(jù)yp=polyfit(x,y,2) %對(duì)x,y用二次多項(xiàng)式擬合,得到系數(shù)px1=0.9 1.2; %輸入點(diǎn)x1y1=polyval(p,x1) %估計(jì)x1處對(duì)應(yīng)的y1 p = -0.2387 0.9191 0.5318y1 = 1.1656 1.2909 使用指定函數(shù)進(jìn)行曲線擬合例2使用試根據(jù)土的壓縮實(shí)驗(yàn)數(shù)據(jù)（數(shù)據(jù)），用雙曲線模型確定模型參數(shù)，（壓力