數學建模教程擬合與插值在大量的應用領域中，人們經常面臨這樣的問題：給定一批數據點，需確定滿足特定要求的曲線或曲面。對這個問題有兩種方法。一種是插值法，數據假定是正確的，要求以某種方法描述數據點之間所發生的情況。另一種方法是曲線擬合或回歸。人們設法找出某條光滑曲線，它最佳地擬合數據，但不必要經過任何數據點。本專題的目的之一是：了解插值和擬合的基本內容及方法；

函數插值與曲線擬合都是要根據一組數據構造一個函數作為近似，由于近似的要求不同，二者的數學方法上是完全不同的。

假設已獲得某函數關系的成批離散實驗數據或觀測數據，擬合問題就是為這樣的大量離散數據建立對應的、近似的連續模型的一種應用基礎問題。所建立的模型的基本形式是一條曲線（一元曲線），稱為擬合曲線或經驗公式。

通常采用“誤差的平方和最小”的原則，即最小二乘擬合問題。

它不要求目標模型（即擬合曲線）精確地過已知的各離散點，只要求目標模型符合已知離散點分布的總體輪廓，并與已知離散點的誤差按某種意義盡量地小。一、擬合問題擬合問題引例1溫度t(0C)20.532.751.073.095.7電阻R()7658268739421032已知熱敏電阻數據：求600C時的電阻R。

設

R=at+ba,b為待定系數擬合問題引例2

t(h)0.250.511.523468c(g/ml)19.2118.1515.3614.1012.899.327.455.243.01已知一室模型快速靜脈注射下的血藥濃度數據(t=0注射300mg)求血藥濃度隨時間的變化規律c(t).作半對數坐標系(semilogy)下的圖形已知一組觀測數據：要求在某特定函數類中尋找一個函數作為的近似函數，使得二者在節點產生的殘差按某種度量標準為最小。常用原則：殘差平方和最小曲線擬合問題的提法線性最小二乘擬合函數的選取

++++++++++++++++++++φ=a1+a2x+++++φ=a1+a2x+a3x2+++++φ=a1+a2x+a3x2φ=a1+a2/xφ=aebxφ=ae-bx1.通過機理分析建立數學模型來確定；2.將數據(xi,yi)i=1,…n作圖，通過直觀判斷確定

：曲線擬合問題最常用的解法——線性最小二乘法的基本思路

第二步:確定a1,a2,…an

的準則（最小二乘準則）：使n個點（xi,yi)與曲線y=φ(x)的距離i的平方和最小

。記

問題歸結為，求

a1,a2,…an

使

J(a1,a2,…an)

最小。第一步:先選定一組函數

使其中

a1,a2,…an

為待定系數。方程組沒有通常意義下的解，這類方程組稱為超定方當線性方程組的方程個數多于未知數的個數時，設線性方程組為程組或矛盾方程組。最小二乘法的求解：預備知識若能找到一組向量令最小，其中使得則稱為該超定方程組的最小二乘解。由多元函數取極值的必要條件有求其最小值。即故得即稱為正則方程組。該方程組的解即為超定方程組的最小二乘解。(2)求解正則方程組得最小二乘解。用最小二乘法解超定方程組的步驟：(1)計算和，得正則方程組。例解得最小二乘解為得方程組解：已知試驗數據，用最小二乘法求擬合直線0.00.20.40.60.80.91.92.83.34.2故得擬合直線可線性化模型的最小二乘擬合很多實際問題中，變量間并非線性關系，但擬合曲線可視為的形式，指數函數如雙曲線即將非線性化問題轉化為線性問題。令則有例給定如下觀測數據，試用指數曲線進行擬合。1.01.251.51.752.05.15.796.537.458.46解：令則有1.01.251.51.752.01.6291.7561.8762.0082.135故解此超定方程組得則擬合曲線為多變量數據擬合有時變量間關系為多元函數關系，有如下觀測數據觀測次數…12m假定變量y與n個變量xi間為線性關系，可設擬合方程為第i組觀測數據對應的殘差為下面考慮用最小二乘原理確定擬合方程的系數ai

。按照最小二乘原理，應使

最小。令求解該超定方程組的最小二乘解即可。用MATLAB解擬合問題1、線性最小二乘擬合2、非線性最小二乘擬合用MATLAB作線性最小二乘擬合1.作多項式f(x)=a1xm+…+amx+am+1擬合,可利用已有命令:a=polyfit(x,y,m)3.對超定方程組可得最小二乘意義下的解。，用2.多項式在x處的值y的計算命令：y=polyval（a，x）輸出擬合多項式系數a=[a1,…,am,am+1]’

（數組）輸入同長度數組X，Y擬合多項式

次數即要求出二次多項式:中的使得:例對下面一組數據作二次多項式擬合1）輸入命令：x=0:0.1:1;y=[-0.447,1.978,3.28,6.16,7.08,7.34,7.66,9.56,9.48,9.30,11.2];R=[(x.^2)',x',ones(11,1)]；

A=R\y'解法1．解超定方程的方法2）計算結果：Ａ=[-9.8108,20.1293,-0.0317]1）輸入命令：x=0:0.1:1;y=[-0.447,1.978,3.28,6.16,7.08,7.34,7.66,9.56,9.48,9.30,11.2];A=polyfit(x,y,2)z=polyval(A,x);plot(x,y,'k+',x,z,'r')%作出數據點和擬合曲線的圖形2）計算結果：Ａ=[-9.8108,20.1293,-0.0317]解法2．用多項式擬合的命令1.lsqcurvefit已知數據點：xdata=（xdata1，xdata2，…，xdatan）

ydata=（ydata1，ydata2，…，ydatan）用MATLAB作非線性最小二乘擬合兩個求非線性最小二乘擬合的函數：lsqcurvefit、lsqnonlin。相同點和不同點：兩個命令都要先建立M-文件fun.m，定義函數f(x)，但定義f(x)的方式不同。

lsqcurvefit用以求含參量x（向量）的向量值函數F(x,xdata)=（F（x，xdata1），…，F（x，xdatan））T使得輸入格式:(1)x=lsqcurvefit(‘fun’,x0,xdata,ydata);(2)x=lsqcurvefit(‘fun’,x0,xdata,ydata,lb,ub);

(3)x=lsqcurvefit(‘fun’,x0,xdata,ydata,lb,ub,options);(4)[x,resnorm]=lsqcurvefit(‘fun’,x0,xdata,ydata,…);(5)[x,resnorm,residual]=lsqcurvefit(‘fun’,x0,xdata,ydata,…);

fun是一個事先建立的定義函數F(x,xdata)

的M-文件,自變量為x和xdata說明：x=lsqcurvefit(‘fun’,x0,xdata,ydata,options);迭代初值已知數據點選項見無約束優化

lsqnonlin用以求含參量x（向量）的向量值函數

f(x)=(f1(x),f2(x),…,fn(x))T

，使得

最小。其中fi（x）=f（x，xdatai，ydatai）

=F(x,xdatai)-ydatai2.lsqnonlin已知數據點：xdata=（xdata1，xdata2，…，xdatan）

ydata=（ydata1，ydata2，…，ydatan）輸入格式：

1）x=lsqnonlin（‘fun’，x0）；

2）x=lsqnonlin

（‘fun’，x0，lb，ub）；

3）x=lsqnonlin

（‘fun’，x0，，lb，ub，options）；

4）[x，resnorm]=lsqnonlin

（‘fun’，x0，…）；

5）[x，resnorm

，residual]=lsqnonlin

（‘fun’，x0，…）；說明：x=lsqnonlin

（‘fun’，x0，options）；fun是一個事先建立的定義函數f(x)的M-文件，自變量為x迭代初值選項見無約束優化例2用下面一組數據擬合中的參數a，b，k該問題即解的最優化問題：

1）編寫M-文件curvefun1.m

functionf=curvefun1(x,tdata)f=x(1)+x(2)*exp(-0.02*x(3)*tdata)%其中x(1)=a;x(2)=b；x(3)=k;2）輸入命令tdata=100:100:1000cdata=1e-03*[4.54,4.99,5.35,5.65,5.90,6.10,6.26,6.39,6.50,6.59];x0=[0.2,0.05,0.05];

x=lsqcurvefit('curvefun1',x0,tdata,cdata)f=curvefun1(x,tdata)

F(x，tdata)=，x=(a，b，k)解法1.用命令lsqcurvefit3）運算結果：f=0.00430.00510.00560.00590.00610.00620.00620.00630.00630.0063x=0.0063-0.00340.25424）結論：a=0.0063,b=-0.0034,k=0.25421）編寫M-文件curvefun2.m

functionf=curvefun2(x)

tdata=100:100:1000;

cdata=1e-03*[4.54,4.99,5.35,5.65,5.90,6.10,6.26,6.39,6.50,6.59];f=x(1)+x(2)*exp(-0.02*x(3)*tdata)-cdata2）輸入命令:

x0=[0.2,0.05,0.05];x=lsqnonlin('curvefun2',x0)f=curvefun2(x)函數curvefun2的自變量是x，cdata和tdata是已知參數，故應將cdata

tdata的值寫在curvefun2.m中解法2用命令lsqnonlin

x=（a，b，k）3）運算結果為

f=1.0e-003*（0.2322-0.1243-0.2495-0.2413-0.1668-0.07240.02410.11590.20300.2792）

x=0.0063-0.00340.2542可以看出，兩個命令的計算結果是相同的。4）結論：即擬合得a=0.0063b=-0.0034k=0.2542設有一個年產240噸的食品加工廠,需要統計其生產費用,但由于該廠各項資料不全,無法統計。在這種情況下,統計部門收集了設備、生產能力和該廠大致相同的五個食品加工廠的產量與生產費用資料如下表：如何根據已知五個食品加工廠的資料較準確地估計該廠的生產費用呢？引例插值問題（2）有的函數甚至沒有表達式，只是一種表格函數，（1）如果函數表達式本身比較復雜，且需要多次實際問題中經常要涉及到函數值的計算問題：重復計算時，計算量會很大；

方便且表達簡單的函數來近似代替，這就是插值問題。實際對于這兩種情況，我們都需要尋找一個計算而我們需要的函數值可能不在該表格中。問題：插值與擬合的區別？一維插值一、插值的定義二、插值的方法三、用Matlab解插值問題拉格朗日插值牛頓插值反插值三次樣條插值分段插值法二維插值一、二維插值定義二、網格節點插值法三、用Matlab解插值問題最鄰近插值分片線性插值雙線性插值網格節點數據的插值散點數據的插值一維插值的定義已知n+1個節點其中互不相同，不妨設求任一插值點處的插值構造一個(相對簡單的)函數通過全部節點,即再用計算插值，即代數插值的唯一性是惟一的。定理個不同節點，滿足插值條件的n次插值多項式當選擇代數多項式作為插值函數時，稱為代數多項式插值。

代數插值一、線性插值（n=1，一次插值）求解L1(x)=a1

x+a0已知使得L1(xi)

=

yi

.(i=0,1)如果令則稱l0(x),l1(x)為x0,x1上的線性插值基函數。這時根據點斜式得到并稱其為一次Lagrange插值多項式。f(x)

≈L1(x)=y0l0(x)+y1l1(x)二、拋物線插值（n=2，二次插值）求解L2(x)=a2x2+a1

x+a0使得L2(xi)=yi

,i=0,1,2.已知拋物插值仿照線性函數的構造方法，構造且其中要求均為二次多項式。設即由求故同理例題求線性插值函數。已知的函數表解：多項式可使用類似方法。分析由個不同節點構造次其中且上述多項式稱為拉格朗日插值多項式。拉格朗日插值多項式插值余項和誤差估計且依賴于。其中設f(x)在[a,b]有n+1階導數，則n次插值多定理項式對任意，有插值余項拉格朗日多項式插值的這種振蕩現象叫Runge現象采用拉格朗日多項式插值：選取不同插值節點個數n+1，其中n為插值多項式的次數，當n分別取2,4,6,8,10時，繪出插值結果圖形.例分段線性插值xjxj-1xj+1x0xnxoy例用分段線性插值法求插值,并觀察插值誤差.在[-6,6]中平均選取41個點作插值,結果如圖示比分段線性插值更光滑。xyxi-1xiab

在數學上，光滑程度的定量描述是：函數(曲線)的k階導數存在且連續，則稱該曲線具有k階光滑性。光滑性的階次越高，則越光滑。是否存在較低次的分段多項式達到較高階光滑性的方法？三次樣條插值就是一個很好的例子。三次樣條插值所謂樣條函數，從數學上說，就是按照一定光滑性它在每個內節點上具有直到m-1階連續導數。要求“裝配”起來的分段多項式。它在每個子區間內是m次多項式，三次樣條插值三次樣條函數滿足：在每個子區間上是一個三次多項式。(1)(2)在每個內節點上具有直到二階的連續導數。(3)上可設故在每一個小區間求解使用三彎矩法例用三次樣條插值選取11個基點計算插值用MATLAB作插值計算一維插值函數：yi=interp1(x，y，xi，'method')插值方法被插值點插值節點xi處的插值結果‘nearest’

：最鄰近插值‘linear’

：線性插值；‘spline’

：三次樣條插值；‘cubic’

：立方插值。缺省時：分段線性插值。注意：所有的插值方法都要求x是單調的，并且xi不能夠超過x的范圍。例：在1-12的11小時內，每隔1小時測量一次溫度，測得的溫度依次為：5，8，9，15，25，29，31，30，22，25，27，24。試估計每隔1/10小時的溫度值。hours=1:12;temps=[589152529313022252724];h=1:0.1:12;t=interp1(hours,temps,h,'spline');(直接輸出數據將是很多的)plot(hours,temps,'+',h,t,hours,temps,'r:')%作圖xlabel('Hour'),ylabel('DegreesCelsius’)xy機翼下輪廓線例已知飛機下輪廓線上數據如下，求x每改變0.1時的y值。二維插值的定義xyO第一種（網格節點）：

已知mn個節點其中互不相同，不妨設構造一個二元函數通過全部已知節點,即再用計算插值，即第二種（散亂節點）：yx0已知n個節點其中互不相同，構造一個二元函數通過全部已知節點,即再用計算插值，即注意：最鄰近插值一般不連續。具有連續性的最簡單的插值是分片線性插值。最鄰近插值xy(x1,y1)(x1,y2)(x2,y1)(x2,y2)O二維或高維情形的最鄰近插值，與被插值點最鄰近的節點的函數值即為所求。將四個插值點（矩形的四個頂點）處的函數值依次簡記為：

分片線性插值xy(xi,yj)(xi,yj+1)(xi+1,yj)(xi+1,yj+1)Of(xi,yj)=f1，f(xi+1,yj)=f2，f(xi+1,yj+1)=f3，f(xi,yj+1)=f4插值函數為：第二片(上三角形區域)：(x,y)滿足插值函數為：注意：(x,y)當然應該是在插值節點所形成的矩形區域內。顯然，分片線性插值函數是連續的；分兩片的函數表達式如下：第一片(下三角形區域)：(x,y)滿足雙線性插值是一片一片的空間二次曲面構成。雙線性插值函數的形式如下：其中有四個待定系數，利用該函數在矩形的四個頂點（插值節點）的函數值，得到四個代數方程，正好確定四個系數。雙線性插值xy(x1,y1)(x1,y2)(x2,y1)(x2,y2)O

要求x0,y0單調；x，y可取為矩陣，或x取行向量，y取為列向量，x,y的值分別不能超出x0,y0的范圍。z=interp2(x0,y0,z0,x,y,’method’)被插值點插值方法用MATLAB作網格節點數據的插值插值節點被插值點的函數值‘nearest’最鄰近插值‘linear’雙線性插值‘cubic’雙三次插值缺省時,雙線性插值例：測得平板表面3*5網格點處的溫度分別為：828180828479636165818484828586試作出平板表面的溫度分布曲面z=f(x,y)的圖形。輸入以下命令：x=1:5;y=1:3;temps=[8281808284;7963616581;8484828586];mesh(x,y,temps)1.先在三維坐標畫出原始數據，畫出粗糙的溫度分布曲圖.再輸入以下命令:xi=1:0.2:5;yi=1:0.2:3;zi=interp2(x,y,temps,xi',yi,'cubic');mesh(xi,yi,zi)畫出插值后的溫度分布曲面圖.2．以平滑數據,在x、y方向上每隔0.2個單位的地方進行插值.

通過此例對最近鄰點插值、雙線性插值方法和雙三次插值方法的插值效果進行比較。原始數據圖最鄰近插值雙線性插值雙三次插值

插值函數griddata格式為:

cz

=griddata（x，y，z，cx，cy，‘method’）用MATLAB作散點數據的插值計算

要求cx取行向量，cy取為列向量。被插值點插值方法插值節點被插值點的函數值‘nearest’最鄰近插值‘linear’雙線性插值‘cubic’雙三次插值'v4'-Matlab提供的插值方法缺省時,雙線性插值一室模型：將整個機體看作一個房室，稱中心室，室內血藥濃度是均勻的。快速靜脈注射后，濃度立即上升；然后迅速下降。當濃度太低時，達不到預期的治療效果；當濃度太高，又可能導致藥物中毒或副作用太強。臨床上，每種藥物有一個最小有效濃度c1和一個最大有效濃度c2。