數(shù)學(xué)建模課件-數(shù)值計(jì)算方法

數(shù)學(xué)建模教程擬合與插值在大量的應(yīng)用領(lǐng)域中，人們經(jīng)常面臨這樣的問題：給定一批數(shù)據(jù)點(diǎn)，需確定滿足特定要求的曲線或曲面。對(duì)這個(gè)問題有兩種方法。一種是插值法，數(shù)據(jù)假定是正確的，要求以某種方法描述數(shù)據(jù)點(diǎn)之間所發(fā)生的情況。另一種方法是曲線擬合或回歸。人們設(shè)法找出某條光滑曲線，它最佳地?cái)M合數(shù)據(jù)，但不必要經(jīng)過任何數(shù)據(jù)點(diǎn)。本專題的目的之一是：了解插值和擬合的基本內(nèi)容及方法；

函數(shù)插值與曲線擬合都是要根據(jù)一組數(shù)據(jù)構(gòu)造一個(gè)函數(shù)作為近似，由于近似的要求不同，二者的數(shù)學(xué)方法上是完全不同的。

假設(shè)已獲得某函數(shù)關(guān)系的成批離散實(shí)驗(yàn)數(shù)據(jù)或觀測數(shù)據(jù)，擬合問題就是為這樣的大量離散數(shù)據(jù)建立對(duì)應(yīng)的、近似的連續(xù)模型的一種應(yīng)用基礎(chǔ)問題。所建立的模型的基本形式是一條曲線（一元曲線），稱為擬合曲線或經(jīng)驗(yàn)公式。

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

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

設(shè)

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

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

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

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

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

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

。記

問題歸結(jié)為，求

a1,a2,…an

使

J(a1,a2,…an)

最小。第一步:先選定一組函數(shù)

使其中

a1,a2,…an

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

。按照最小二乘原理，應(yīng)使

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

（數(shù)組）輸入同長度數(shù)組X，Y擬合多項(xiàng)式

次數(shù)即要求出二次多項(xiàng)式:中的使得:例對(duì)下面一組數(shù)據(jù)作二次多項(xiàng)式擬合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）計(jì)算結(jié)果：Ａ=[-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')%作出數(shù)據(jù)點(diǎn)和擬合曲線的圖形2）計(jì)算結(jié)果：Ａ=[-9.8108,20.1293,-0.0317]解法2．用多項(xiàng)式擬合的命令1.lsqcurvefit已知數(shù)據(jù)點(diǎn)：xdata=（xdata1，xdata2，…，xdatan）

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

lsqcurvefit用以求含參量x（向量）的向量值函數(shù)F(x,xdata)=（F（x，xdata1），…，F(xiàn)（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是一個(gè)事先建立的定義函數(shù)F(x,xdata)

的M-文件,自變量為x和xdata說明：x=lsqcurvefit(‘fun’,x0,xdata,ydata,options);迭代初值已知數(shù)據(jù)點(diǎn)選項(xiàng)見無約束優(yōu)化

lsqnonlin用以求含參量x（向量）的向量值函數(shù)

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

，使得

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

=F(x,xdatai)-ydatai2.lsqnonlin已知數(shù)據(jù)點(diǎn)：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是一個(gè)事先建立的定義函數(shù)f(x)的M-文件，自變量為x迭代初值選項(xiàng)見無約束優(yōu)化例2用下面一組數(shù)據(jù)擬合中的參數(shù)a，b，k該問題即解的最優(yōu)化問題：

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）運(yùn)算結(jié)果：f=0.00430.00510.00560.00590.00610.00620.00620.00630.00630.0063x=0.0063-0.00340.25424）結(jié)論：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)函數(shù)curvefun2的自變量是x，cdata和tdata是已知參數(shù)，故應(yīng)將cdata

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

x=（a，b，k）3）運(yùn)算結(jié)果為

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可以看出，兩個(gè)命令的計(jì)算結(jié)果是相同的。4）結(jié)論：即擬合得a=0.0063b=-0.0034k=0.2542設(shè)有一個(gè)年產(chǎn)240噸的食品加工廠,需要統(tǒng)計(jì)其生產(chǎn)費(fèi)用,但由于該廠各項(xiàng)資料不全,無法統(tǒng)計(jì)。在這種情況下,統(tǒng)計(jì)部門收集了設(shè)備、生產(chǎn)能力和該廠大致相同的五個(gè)食品加工廠的產(chǎn)量與生產(chǎn)費(fèi)用資料如下表：如何根據(jù)已知五個(gè)食品加工廠的資料較準(zhǔn)確地估計(jì)該廠的生產(chǎn)費(fèi)用呢？引例插值問題（2）有的函數(shù)甚至沒有表達(dá)式，只是一種表格函數(shù)，（1）如果函數(shù)表達(dá)式本身比較復(fù)雜，且需要多次實(shí)際問題中經(jīng)常要涉及到函數(shù)值的計(jì)算問題：重復(fù)計(jì)算時(shí)，計(jì)算量會(huì)很大；

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

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

x+a0已知使得L1(xi)

=

yi

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

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

x+a0使得L2(xi)=yi

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

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

：最鄰近插值‘linear’

：線性插值；‘spline’

：三次樣條插值；‘cubic’

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

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

分片線性插值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插值函數(shù)為：第二片(上三角形區(qū)域)：(x,y)滿足插值函數(shù)為：注意：(x,y)當(dāng)然應(yīng)該是在插值節(jié)點(diǎn)所形成的矩形區(qū)域內(nèi)。顯然，分片線性插值函數(shù)是連續(xù)的；分兩片的函數(shù)表達(dá)式如下：第一片(下三角形區(qū)域)：(x,y)滿足雙線性插值是一片一片的空間二次曲面構(gòu)成。雙線性插值函數(shù)的形式如下：其中有四個(gè)待定系數(shù)，利用該函數(shù)在矩形的四個(gè)頂點(diǎn)（插值節(jié)點(diǎn)）的函數(shù)值，得到四個(gè)代數(shù)方程，正好確定四個(gè)系數(shù)。雙線性插值xy(x1,y1)(x1,y2)(x2,y1)(x2,y2)O

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

通過此例對(duì)最近鄰點(diǎn)插值、雙線性插值方法和雙三次插值方法的插值效果進(jìn)行比較。原始數(shù)據(jù)圖最鄰近插值雙線性插值雙三次插值

插值函數(shù)griddata格式為:

cz

=griddata（x，y，z，cx，cy，‘method’）用MATLAB作散點(diǎn)數(shù)據(jù)的插值計(jì)算

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