插值與數據擬合建模

插值與數據擬合就是通過一些已知數據去確定某類函數的參數或尋找某個近似函數，使所得的函數與已知數據具有較高的精度，并且能夠使用數學分析的工具分析數據所反映的對象的性質．幾種常用的方法：

1、一般插值法

2、樣條插值法

3、最小二乘曲線

4、曲面的擬合數據擬合與插值建模

上大學二年級的小華正在做作業(yè)，“爸爸，計算這道題要用到sin，可是我的計算器壞了，怎么辦?！碑敼こ處煹睦蠌垙暮窈竦囊晦f書底下抽出一本數學用表來，“給你，這是我念大學時用的，那時候啊，計算器聽都沒聽說過。”小華拿著表翻了一會兒，無奈的說：“表上每才有一個函數值，這里只sin和sin““表中沒有的，都可以用插值方法計算”“插值！我們的數學實驗課就要學了，不過，今天我要先自己想個辦法，用這個算出sin”這本四位數學用表給出sin＝0.576，sin=0.5783。小華認為在sin到sin這樣小的范圍內，正弦可以近似為線性函數，于是很容易地得到Sin=0.576+(0.5783-0.5760)×0.6=0.5774

聰明的小華用的這個辦法是一種插值方法——分段線性插值。實際上，插值可以理解為，要根據一個用表格表示的函數，計算表中沒有的函數值。表中有的，如（sin，0.5760）(sin,0.5783)稱為節(jié)點；要計算的，如sin，稱為插值點，結果（0.5774）即為插值。小華作的線性函數為插值函數，插值函數所表示的直線當然要通過節(jié)點。

插值最初來源于天體計算——由若干觀測值（即節(jié)點）計算任意時刻星球的位置（即插值點和插值）——的需要?，F在，雖然人們已很少需要用它從函數表計算函數值了，但是插值仍然在諸如機械加工等工程技術和數據處理等科學研究中有著許多直接的應用，另一方面，插值又是數值微分、數值積分、常微分方程數值等數值計算的基礎。

幾天后，小華在物理實驗里又碰到一個看起來非常類似的問題：有一只對溫度敏感的電阻，已經測得了一組溫度T和電阻R數據。

現在想知道600C時的電阻多大。溫度t(0C)20.532.751.073.095.7電阻R()7658268739421032

小華征求老師的意見，老師給了他兩點提示：一是在直角坐標系中把5個點（T，R）畫一下，看看電阻R和溫度T之間大致有什么關系；二是測量數據總有相當大的誤差，這與用函數表作插值計算應該有不同之處吧（雖然函數表也存在舍入誤差，但很小，可以認為表中數值是精確的）通過圖形小華看到，R與T大致呈直線關系，于是用手畫了一條靠近5個點的直線，又想起中學物理學過，金屬材料的電阻率與溫度成正比，從而確定R與T的關系應該是

R=at+b

其中a，b為待定常數。

正是由于測量誤差的存在，由R=at+b表示的直線不可能通過全部5個點，所以，與插值曲線要通過全部節(jié)點不同，小華打算作一條盡量靠近所有的點的直線，求出a，b待定常數，由此計算t=600C

的R就十分簡單了。

根據一組（二組）數據，即平面上的若干點，確定一個一元函數，即曲線，使這些節(jié)點與曲線總體來說盡量接近，這就是曲線擬合。函數值與曲線擬合都是要根據一組數據構造一個函數作為近似，由于近似的要求不同，二者的數學方法是完全不同的。

數據擬合1.擬合的基本原理；2.最小二乘擬合；3.用Matlab作最小二乘擬合；4.如何用擬合解決實際問題。

t(h)0.250.511.523468c(g/ml)19.2118.1515.3614.1012.899.327.455.243.01

對某人用快速靜脈注射方式一次性注射某種藥物300mg后，經過時間t采集血樣，測得血藥濃度c如下表：求血藥濃度隨時間的變化規(guī)律c(t).半對數坐標系(semilogy)下的圖形Log10c(t)=at+b引例1：血藥濃度的變化規(guī)律曲線擬合問題的提法

已知一組（二維）數據，即平面上n個點（xi,yi)i=1,…n,

尋求一個函數（曲線）y=f(x),

使f(x)

在某種準則下與所有數據點最為接近，即曲線擬合得最好。

+++++++++xyy=f(x)(xi,yi)

i

i為點（xi,yi)與曲線y=f(x)的距離最小二乘擬合

第一步：先選定一類函數f(x，a1，a2，

…，am)

其準則為（最小二乘準則）：使n個點（xi,yi)與曲線y=f(x，a1，a2，

…，am)的距離

i的平方和最小

。其中

a1,a2,…am

為待定常數。f可以為一些簡單的“基函數”（如冪函數，三角函數等等）的線性組合：第二步：確定參數a1,a2,…am，記問題歸結為，求

a1,a2,…am

使

J(a1,a2,…am)

最小。這樣的擬合稱為最小二乘擬合。

除了最小二乘準則（即各點誤差的平方和最?。?，你認為還可以用怎樣的擬合準則？比較起來，最小二乘準則有什么優(yōu)點？思考最小二乘擬合函數f(x，a1，…am)的選取++++++++++++++++++++f=a1+a2xf=a1+a2x+a3x2f=a1+a2x+a3x2f=a1exp(a2x)+++++f=a1exp(a2x)1.通過機理分析建立數學模型來確定f；2.將數據(xi,yi)i=1,…n作圖，通過直觀判斷確定f：2.作一般的最小二乘曲線擬合，可利用已有程序curvefit,其調用格式為：

a=curvefit(‘f’,a0,x,y)

1.作多項式f(x)=a1xm+…+amx+am+1函數擬合,可利用已有程序polyfit,其調用格式為:a=polyfit(x,y,m)用MATLAB作最小二乘擬合數據點擬合多項式次數系數注：f為擬合函數y=f(a,x)的函數M—文件，f(a,x)為擬合函數。數據點待定常數a的初值函數M文件用MATLAB作多項式最小二乘擬合2.用命令polyfit(x,y,m)得到a1=3.3940,a2=702.49181.選取函數R=

a1t+a2溫度t(0C)20.532.751.073.095.7電阻R()7658268739421032例.由數據擬合R=f(t)用MATLAB作最小二乘曲線擬合例：用函數f(x)=a1*exp(-a2*x)+a3*exp(-a4*x)擬合下列數據點：xdata=[0:.1:2]ydata=[5.89553.56392.51731.97901.89901.39381.13591.00961.03430.84350.68560.61000.53920.39460.39030.54740.34590.13700.22110.17040.2636]用命令curvefit(‘f’,a0,x,y)

擬合的應用——參數辨識

數學建模的方法：機理分析和測試分析。機理分析是根據對客觀事物特性的認識，找出反映內部機理的數量規(guī)律，建立的模型常有明確的物理意義。測試分析將研究的對象看作一個“黑箱”，通過對實驗數據的統計分析，找出與數據擬合得最好的模型。機理分析——>模型結構實驗數據——>未知參數范例：薄膜滲透率的測定

一、問題：某種醫(yī)用薄膜，具有從高濃度的溶液向低濃度的溶液擴散的功能，在試制時需測定薄膜被物質分子穿透的能力。測定方法：用面積為S的薄膜將容器分成體積分別為的兩部份，在兩部分中分別注滿該物質的兩種不同濃度的溶液。此時該物質分子就會從高濃度溶液穿過薄膜向低濃度溶液中擴散。平均每單位時間通過單位面積薄膜的物質分子量與膜兩側溶液的濃度差成正比，比例系數K表征了薄膜被該物質分子穿透的能力，稱為滲透率。定時測量容器中薄膜某一側的溶液濃度，以此確定K。VAVBS二、問題分析考察時段[t，t+Δt]薄膜兩側容器中該物質質量的變化。設，對容器的B部分溶液濃度的測試結果如下表：（濃度單位）

1）在容器的一側，物質質量的增加是由于另一側的物質向該側滲透的結果，因此物質質量的增量應等于另一側的該物質向這側的滲透量。

以容器A側為例，在時段[t，t+Δt]物質質量的增量為：分別表示在時刻t膜兩側溶液設的濃度，濃度單位:

由于平均每單位時間通過單位面積薄膜的物質分子量與膜兩側溶液的濃度差成正比，比例系數為K。

因此，在時段[t，t+Δt]，從B側滲透至A側的該物質的質量為：于是有：兩邊除以Δt，并令Δt→0取極限再稍加整理即得：分別表示在初始時刻兩側溶液的濃度其中（1）2)注意到整個容器的溶液中含有該物質的質量不變,與初始時刻該物質的含量相同，因此從而：加上初值條件：代入式（1）得：便可得出CB(t)的變化規(guī)律，從而根據實驗數據進行擬合，估計出參數K,。三、數學模型假設：1）薄膜兩側的溶液始終是均勻的；2）平均每單位時間通過單位面積薄膜的物質分子量與膜兩側溶液的濃度差成正比。3）薄膜是雙向同性的即物質從膜的任何一側向另一側滲透的性能是相同的?；诩僭O和前面的分析，B側的濃度CB(t)應滿足如下微分方程和初始條件：四、求解方法：1.函數擬合法前面得到的模型是一個帶初值的一階線性微分方程，解之得：問題歸結為利用CB在時刻tj的測量數據Cj(j=1,2,...,N)來辨識K和。引入從而

用函數CB(t)來擬合所給的實驗數據，從而估計出其中的參數a,b,K。將代入上式有：用MATLAB軟件進行計算.1）編寫函數M-文件nongdu.mfunctionf=nongdu(x,tdata)f=x(1)+x(2)*exp(-0.02*x(3)*tdata);其中x(1)=a;x(2)=b;x(3)=k;2)在工作空間中執(zhí)行以下命令(test1.m)tdata=linspace(100,1000,10);cdata=[4.544.995.355.655.906.10...6.266.396.506.59];x0=[0.2,0.05,0.05];x=curvefit(‘nongdu’,x0,tdata,cdata)3)輸出結果:x=0.007-0.0030.1012

即k=0.1012,a=0.007,b=-0.003,進一步求得：2.非線性規(guī)劃法

利用CB在時刻tj的測量數據Cj(j=1,2,...,N)來辨識K和。問題可轉化為求函數即求函數的最小值點（K，a，b）。3.導函數擬合法前面得到的微分方程為：令上式變?yōu)椋哼@可以看作隨CB的變化規(guī)律(j=1,2,...,N)若知道一組數據則可用最小二乘擬合的方法來求出函數中的未知參數K和h。即為求參數K,a使下列誤差函數達到最小：

該問題等價于用函數f(K,a,CB)=K(0.01a-0.02CB)來擬合數據(j=1,2,...,N)用MATLAB軟件進行計算.％求數據點(j=1,2,...,N)tdata=linspace(100,1000,10);cdata=1e-05.*[454499535565590...

610626639650659];[d,ifail]=e01bef(tdata,cdata);[cj,dcj]=e01bgf(tdata,cdata,d,tdata);1）編寫函數M-文件baomof.mfunctionf=baomof(x,cdata)f=x(1)*(0.01*x(2)-0.02*cdata)其中x(1)=K;x(2)=h2)編寫命令M文件(baomo21.m)3)輸出結果:x=0.10090.014

即k=0.1009,h=0.014％作函數擬合x0=[0.2,0.1];x=curvefit('baomof',x0,cdata,dcj')4.線性化迭代法前面帶初始條件的一階線性微分方程的解為其中：

如果得到了參數K的一個較好的近似值K*，則將CB(t)關于K在K*處展開，略去K的二次及以上的項得CB(t)的一個近似式通過極小化確定a,b,d,再由

K=d/0.02b得到K*的修正值

K。K*K*-K,得到K的一個新的近似值，用同樣的方法再求新的修正值

K。這個過程可以不斷重復，直到修正值足夠小為止。1）當K的初值取為k=0.3時，出現奇異情況，迭代不收斂；2）當K的初值取為k=0.2時，經四次迭代，已經收斂到一個很好的解。迭代結果如下表。五、結果及誤差分析

幾種方法得出的結果及相應的誤差總結于下表，誤差為計算數據與實驗數據之差的平方和。注：導函數擬合法得出的參數值精度有限，線性化迭代法要求參數的初值比較接近精確值。因此可將導函數擬合法和線性化迭代法結合起來使用，把前者得到的參數K的值作為迭代法中K的初值，這樣可使迭代法收斂或收斂更快。3）取K的初值為k=0.1009,只一次迭代就得到2）中的最后結果。函數擬合法的擬合效果求解參數辨識模型的方法：函數擬合；非線性規(guī)劃；導函數擬合；線性化迭代；其它方法。用Logistic模擬水稻葉伸長生長時間11.82.63.44.14.85.46.16.87.48.1重量0.30.50.91.42.53.24.37.610.114.418.5時間8.89.410.110.811.712.413.114.415.115.7

重量23.025.230.433.738.841.743.744.845.545.3

生長觀測記錄數據模型表達式：程序！關于polyfit命令命令：p=polyfit(x,y,n)（1）x與y為模擬數據（2）n為擬合多項式的次數（3）當n=1時為用最小二乘法進行直線擬合（4）得到的向量p為長度n+1向量，對應p的分量依次是次數從高到底各多項式系數用Richard模擬

水稻葉伸長生長關于inline函數例如：y=inline(‘sin(x)-cos(x)’,’x’)輸入y(0)，可得：-1作圖：x=0:0.1:2*pi;plot(x,y(x))1、插值問題：不知道某一函數f(x)在待定范圍[a,b]上的具體表達式，而只能通過實驗測量得到該函數在一系列點a≤x1,x2,...,xn≤b上的值

y0,y1,y2,...,yn，需要找一個簡單的函數P(x)

來近似地代替f(x)，要求滿足：

P(xi)=yi(i=1,2,...,n)2、插值函數：P(x)3、插值法：求插值函數P(x)的方法

插值二、常用插值函數1、多項式函數2、樣條函數1、多項式插值方法（1）n次代數插值（2）拉格朗日插值幾點說明：（1）拉格朗日插值基函數僅與節(jié)點有關，而與被插值函數f(x)無關。（2）拉格朗日插值多項式僅由數對(xi,yi)(i＝

1,2,…,n)確定，而與數對排列次序無關。（3）多項式插值除了上述插值法外還有其它插值法，如newton插值法、hermite插值法等。2、樣條插值方法（1）樣條函數——m次半截冪函數（2）k次B樣條或k次基本樣條函數的定義（一）廣泛使用的樣條函數（1）廣泛采用：二次樣條、三次樣條及B樣條（2）力學意義：

A:二次樣條在力學上解釋為集中力偶作用下的彈性細梁撓度曲線。

B:彈性細梁受集中載荷作用形成的撓度曲線，在小撓度的情況下，恰好表示為三次樣條函數，集中載荷的作用點，恰好就是三次樣條函數的節(jié)點。（1）二次樣條的定義

設[a,b]的一個劃分：a=x0<x1,x2,...,xn=b，函數f(x)各節(jié)點的值分別為：

f(xi)=yi(i=1,2,...,n)

如果二次樣條函數：滿足：S(xi)=yi(i=1,2,...,n)（2）三次樣條函數的定義

設[a,b]的一個劃分：a=x0<x1,x2,...,xn=b，函數f(x)各節(jié)點的值分別為：

f(xi)=yi(i=1,2,...,n)

如果三次樣條函數：3滿足：S(xi)=yi(i=1,2,...,n)例：某實驗對一根長10米的鋼軌進行熱源的溫度傳播測試。用x表示測量點0:2.5:10(米)，用h表示測量時間0:30:60(秒)，用T表示測試所得各點的溫度(℃)。試用線性插值求出在一分鐘內每隔20秒、鋼軌每隔1米處的溫度TI。

命令如下：

x=0:2.5:10;

h=[0:30:60]';

T=[95,14,0,0,0;88,48,32,12,6;67,64,54,48,41];

xi=[0:10];

hi=[0:20:60]';

TI=interp2(x,h,T,xi,hi)例：設z=x2+y2，對z函數在[0,1]×[0,2]區(qū)域內進行插值。為了說明這個維數的插值，再考慮一個問題。設人們對平板上的溫度分布估計感興趣，給定的溫度值取自平板表面均勻分布的格柵。采集了下列的數據：

?width=1:5;%indexforwidthofplate(i.e.,thex-dimension)

?depth=1:3;%indexfordepthofplate(i,e,,they-dimension)

?temps=[8281808284;7963616581;8484828586]%temperaturedata temps= 8281808284 7963616581 8484828586

如同在標引點上測量一樣，矩陣temps表示整個平板的溫度分布。temps的列與下標depth或y-維相聯系，行與下標width或x-維相聯系(見圖6)。為了估計

在中間點的溫度，我們必須對它們進行辨識。

?wi=1:0.2:5;%estimateacrosswidthofplate

?d=2;%atadepthof2

?zlinear=interp2(width,depth,temps,wi,d);%linearinterpolation

?zcubic=interp2(width,depth,temps,wi,d,'cubic');%cubicinterpolation

?plot(wi,zlinear,'-',wi,zcubic)%plotresults

?xlabel('WidthofPlate'),ylabel('DegreesCelsius')

?title(['TemperatureatDepth='num2str(d)])圖6在深度d=2處的平板溫度

另一種方法，我們可以在兩個方向插值。先在三維坐標畫出原始數據，看一下該數據的粗糙程度(見圖7)。

?mesh(width,depth,temps)%usemeshplot

?xlabel(‘WidthofPlate’),ylabel(‘DepthofPlate’)

?zlabel(‘DegreesCelsius’),axis(‘ij’),grid

[xi,yi]=meshgrid(width,depth);

zi=interp2(width,depth,temps,xi,yi,‘cubic’);

mesh(xi,yi,zi)

圖7平板溫度然后在兩個方向上插值，以平滑數據(見圖8)。

?di=1:0.2:3;%choosehigherresolutionfordepth

?wi=1:0.2:5;%choosehigherresolutionforwidth

?zcubic=interp2(width,depth,temps,wi,di,'cubic');%cubic

?mesh(wi,di,zcubic)

?xlabel('WidthofPlate'),ylabel('DepthofPlate')

?zlabel('DegreesCelsius'),axis('ij'),grid圖8二維插值后的平板溫度

上面的例子清楚地證明了，二維插值更為復雜，只是因為有更多的量要保持跟蹤。interp2的基本形式是interp2(x,y,z,xi,yi,method)。這里x和y是兩個獨立變量，z是一個應變量矩陣。x和y對z的關系是

z(i,:)=f(x,y(i))和z(:,j)=f(x(j),y).

也就是，當x變化時，z的第i行與y的第i個元素y(i)相關，當y變化時，z的第j列與x的第j個元素x(j)相關。xi是沿x-軸插值的一個數值數組；yi是沿y-軸插值的一個數值數組。

可選的參數method可以是‘linear’，‘cubic’或‘nearest’。在這種情況下，cubic不意味著3次樣條，而是使用3次多項式的另一種算法。linear方法是線性插值，僅用作連接圖上數據點。nearest方法只選擇最接近各估計點的粗略數據點。

插值的優(yōu)點：利用已知點確定未知點粗糙——精確集合大的——簡化的例：某觀測站測得某日6:00時至18:00時之間每隔2小時的室內外溫度(℃)，用3次樣條插值分別求得該日室內外6:30至17:30時之間每隔2小時各點的近似溫度(℃)。

設時間變量h為一行向量，溫度變量t為一個兩列矩陣，其中第一列存放室內溫度，第二列儲存室外溫度。命令如下：

h=6:2:18;

t=[18,20,22,25,30,28,24;15,19,24,28,34,32,30]';

XI=6.5:2:17.5

YI=interp1(h,t,XI,‘spline’)

%用3次樣條插值計算

在討論二維插值之前，強調interp1所強制的二個強約束是很重要的。首先，人們不能要求有獨立變量范圍以外的結果，例如，interp1(hours,temps,13.5)導致一個錯誤，因為hours在1到12之間變化。其次，獨立變量必須是單調的。即獨立變量在值上必須總是增加的或總是減小的。

在我們的例子里，hours是單調的。然而，如果我們已經定義獨立變量為一天的實際時間，

?time_of_day=[7:121:6]%startat

7AM,endat6PM

time_of_day=

789101112123456

則獨立變量將不是單調的，因為time_of_day增加到12，然后跌到1，再然后增加。如果用time_of_day代替interp1中的hours，將會返回一個錯誤。同樣的理由，人們不能對temps插值來找出產生某溫度的時間(小時)，因為temps不是單調的。案例3估計水箱的水流量模型長度單位：E(＝30．24cm)容積單位：G(=3.785L(升))