第九章 插值與擬合

1、第九章 插值與擬合插值：求過(guò)已知有限個(gè)數(shù)據(jù)點(diǎn)的近似函數(shù)。擬合：已知有限個(gè)數(shù)據(jù)點(diǎn)，求近似函數(shù)，不要求過(guò)已知數(shù)據(jù)點(diǎn)，只要求在某種意義下它在這些點(diǎn)上的總偏差最小。插值和擬合都是要根據(jù)一組數(shù)據(jù)構(gòu)造一個(gè)函數(shù)作為近似，由于近似的要求不同，二者的數(shù)學(xué)方法上是完全不同的。而面對(duì)一個(gè)實(shí)際問(wèn)題，究竟應(yīng)該用插值還是擬合，有時(shí)容易確定，有時(shí)則并不明顯。§1 插值方法下面介紹幾種基本的、常用的插值：拉格朗日多項(xiàng)式插值、牛頓插值、分段線性插值、Hermite插值和三次樣條插值。1.1 拉格朗日多項(xiàng)式插值 1.1.1 插值多項(xiàng)式用多項(xiàng)式作為研究插值的工具，稱為代數(shù)插值。其基本問(wèn)題是：已知函數(shù)在區(qū)間上個(gè)不同點(diǎn)處的函

2、數(shù)值，求一個(gè)至多次多項(xiàng)式 （1）使其在給定點(diǎn)處與同值，即滿足插值條件 （2）稱為插值多項(xiàng)式，稱為插值節(jié)點(diǎn)，簡(jiǎn)稱節(jié)點(diǎn)，稱為插值區(qū)間。從幾何上看，次多項(xiàng)式插值就是過(guò)個(gè)點(diǎn)，作一條多項(xiàng)式曲線近似曲線。次多項(xiàng)式（1）有個(gè)待定系數(shù)，由插值條件（2）恰好給出個(gè)方程 （3）記此方程組的系數(shù)矩陣為，則 是范德蒙特(Vandermonde)行列式。當(dāng)互不相同時(shí)，此行列式值不為零。因此方程組（3）有唯一解。這表明，只要個(gè)節(jié)點(diǎn)互不相同，滿足插值要求（2）的插值多項(xiàng)式（1）是唯一的。插值多項(xiàng)式與被插函數(shù)之間的差稱為截?cái)嗾`差，又稱為插值余項(xiàng)。當(dāng)充分光滑時(shí)，其中。1.1.2 拉格朗日插值多項(xiàng)式實(shí)際上比較方便的作法不是解方程

3、（3）求待定系數(shù)，而是先構(gòu)造一組基函數(shù) 是次多項(xiàng)式，滿足令 （4）上式稱為次Lagrange插值多項(xiàng)式，由方程（3）解的唯一性，個(gè)節(jié)點(diǎn)的次Lagrange插值多項(xiàng)式存在唯一。 1.1.3 用Matlab作Lagrange插值Matlab中沒(méi)有現(xiàn)成的Lagrange插值函數(shù)，必須編寫(xiě)一個(gè)M文件實(shí)現(xiàn)Lagrange插值。設(shè)個(gè)節(jié)點(diǎn)數(shù)據(jù)以數(shù)組輸入（注意Matlat的數(shù)組下標(biāo)從1開(kāi)始），個(gè)插值點(diǎn)以數(shù)組輸入，輸出數(shù)組為個(gè)插值。編寫(xiě)一個(gè)名為lagrange.m的M文件：function y=lagrange(x0,y0,x);n=length(x0);m=length(x);for i=1:m z=x(i)

4、; s=0.0; for k=1:n p=1.0; for j=1:n if j=k p=p*(z-x0(j)/(x0(k)-x0(j); end end s=p*y0(k)+s; end y(i)=s;end 1.2 牛頓（Newton）插值在導(dǎo)出Newton公式前，先介紹公式表示中所需要用到的差商、差分的概念及性質(zhì)。1.2.1 差商定義 設(shè)有函數(shù)為一系列互不相等的點(diǎn)，稱為關(guān)于點(diǎn)一階差商（也稱均差）記為，即 稱一階差商的差商 為關(guān)于點(diǎn)的二階差商，記為。一般地，稱 為關(guān)于點(diǎn)的階差商，記為 容易證明，差商具有下述性質(zhì)： 1.2.2 Newton插值公式線性插值公式可表成 稱為一次Newton插值

5、多項(xiàng)式。一般地，由各階差商的定義，依次可得 將以上各式分別乘以，然后相加并消去兩邊相等的部分，即得 記顯然，是至多次的多項(xiàng)式，且滿足插值條件，因而它是的次插值多項(xiàng)式。這種形式的插值多項(xiàng)式稱為Newton插值多項(xiàng)式。稱為Newton插值余項(xiàng)。Newton插值的優(yōu)點(diǎn)是：每增加一個(gè)節(jié)點(diǎn)，插值多項(xiàng)式只增加一項(xiàng)，即因而便于遞推運(yùn)算。而且Newton插值的計(jì)算量小于Lagrange插值。由插值多項(xiàng)式的唯一性可知，Newton插值余項(xiàng)與Lagrange余項(xiàng)也是相等的，即由此可得差商與導(dǎo)數(shù)的關(guān)系其中。1.2.3 差分當(dāng)節(jié)點(diǎn)等距時(shí)，即相鄰兩個(gè)節(jié)點(diǎn)之差（稱為步長(zhǎng)）為常數(shù)，Newton插值公式的形式會(huì)更簡(jiǎn)單。此時(shí)關(guān)

6、于節(jié)點(diǎn)間函數(shù)的平均變化率（差商）可用函數(shù)值之差（差分）來(lái)表示。定義 設(shè)有等距節(jié)點(diǎn)，步長(zhǎng)為常數(shù)，。稱相鄰兩個(gè)節(jié)點(diǎn)處的函數(shù)值的增量為函數(shù)在點(diǎn)處以為步長(zhǎng)的一階差分，記為，即 類似地，定義差分的差分為高階差分。如二階差分為 一般地，階差分為 ，上面定義的各階差分又稱為向前差分。常用的差分還有兩種： 稱為在處以為步長(zhǎng)的向后差分； 稱為在處以為步長(zhǎng)的中心差分。一般地，階向后差分與階中心差分公式為 差分具有以下性質(zhì)： （i）各階差分均可表成函數(shù)值的線性組合，例如 （ii）各種差分之間可以互化。向后差分與中心差分化成向前差分的公式如下： 1.2.4 等距節(jié)點(diǎn)插值公式如果插值節(jié)點(diǎn)是等距的，則插值公式可用差分表示

7、。設(shè)已知節(jié)點(diǎn)，則有 若令，則上式又可變形為上式稱為Newton向前插值公式。1.3 分段線性插值1.3.1 插值多項(xiàng)式的振蕩用Lagrange插值多項(xiàng)式近似，雖然隨著節(jié)點(diǎn)個(gè)數(shù)的增加，的次數(shù)變大，多數(shù)情況下誤差會(huì)變小。但是增大時(shí)，的光滑性變壞，有時(shí)會(huì)出現(xiàn)很大的振蕩。理論上，當(dāng)，在內(nèi)并不能保證處處收斂于。Runge給出了一個(gè)有名的例子：對(duì)于較大的，隨著的增大，振蕩越來(lái)越大，事實(shí)上可以證明，僅當(dāng)時(shí)，才有，而在此區(qū)間外，是發(fā)散的。高次插值多項(xiàng)式的這些缺陷，促使人們轉(zhuǎn)而尋求簡(jiǎn)單的低次多項(xiàng)式插值。1.3.2 分段線性插值簡(jiǎn)單地說(shuō)，將每?jī)蓚€(gè)相鄰的節(jié)點(diǎn)用直線連起來(lái)，如此形成的一條折線就是分段線性插值函數(shù)，記作

8、，它滿足，且在每個(gè)小區(qū)間上是線性函數(shù)?？梢员硎緸?有良好的收斂性，即對(duì)于有，。用計(jì)算點(diǎn)的插值時(shí)，只用到左右的兩個(gè)節(jié)點(diǎn)，計(jì)算量與節(jié)點(diǎn)個(gè)數(shù)無(wú)關(guān)。但越大，分段越多，插值誤差越小。實(shí)際上用函數(shù)表作插值計(jì)算時(shí)，分段線性插值就足夠了，如數(shù)學(xué)、物理中用的特殊函數(shù)表，數(shù)理統(tǒng)計(jì)中用的概率分布表等。1.3.3 用Matlab實(shí)現(xiàn)分段線性插值用Matlab實(shí)現(xiàn)分段線性插值不需要編制函數(shù)程序，Matlab中有現(xiàn)成的一維插值函數(shù)interp1。y=interp1(x0,y0,x,'method')method指定插值的方法，默認(rèn)為線性插值。其值可為：'nearest' 最近項(xiàng)插值'

9、;linear' 線性插值'spline' 立方樣條插值'cubic' 立方插值。所有的插值方法要求x0是單調(diào)的。當(dāng)x0為等距時(shí)可以用快速插值法，使用快速插值法的格式為'*nearest'、'*linear'、'*spline'、'*cubic'。1.4 埃爾米特(Hermite)插值1.4.1 Hermite插值多項(xiàng)式如果對(duì)插值函數(shù)，不僅要求它在節(jié)點(diǎn)處與函數(shù)同值，而且要求它與函數(shù)有相同的一階、二階甚至更高階的導(dǎo)數(shù)值，這就是Hermite插值問(wèn)題。本節(jié)主要討論在節(jié)點(diǎn)處插值函數(shù)與函數(shù)的值及一階

10、導(dǎo)數(shù)值均相等的Hermite插值。設(shè)已知函數(shù)在個(gè)互異節(jié)點(diǎn)上的函數(shù)值 和導(dǎo)數(shù)值 ，要求一個(gè)至多次的多項(xiàng)式，使得 滿足上述條件的多項(xiàng)式稱為Hermite插值多項(xiàng)式。Hermite插值多項(xiàng)式為其中，。1.4.2 用Matlab實(shí)現(xiàn)Hermite插值Matlab中沒(méi)有現(xiàn)成的Hermite插值函數(shù)，必須編寫(xiě)一個(gè)M文件實(shí)現(xiàn)插值。設(shè)個(gè)節(jié)點(diǎn)的數(shù)據(jù)以數(shù)組（已知點(diǎn)的橫坐標(biāo)），（函數(shù)值），（導(dǎo)數(shù)值）輸入（注意Matlat的數(shù)組下標(biāo)從1開(kāi)始），個(gè)插值點(diǎn)以數(shù)組輸入，輸出數(shù)組為個(gè)插值。編寫(xiě)一個(gè)名為hermite.m的M文件：function y=hermite(x0,y0,y1,x);n=length(x0);m=len

11、gth(x);for k=1:m yy=0.0; for i=1:n h=1.0; a=0.0; for j=1:n if j=i h=h*(x(k)-x0(j)/(x0(i)-x0(j)2; a=1/(x0(i)-x0(j)+a; end end yy=yy+h*(x0(i)-x(k)*(2*a*y0(i)-y1(i)+y0(i); end y(k)=yy;end1.5 樣條插值許多工程技術(shù)中提出的計(jì)算問(wèn)題對(duì)插值函數(shù)的光滑性有較高要求，如飛機(jī)的機(jī)翼外形，內(nèi)燃機(jī)的進(jìn)、排氣門(mén)的凸輪曲線，都要求曲線具有較高的光滑程度，不僅要連續(xù)，而且要有連續(xù)的曲率，這就導(dǎo)致了樣條插值的產(chǎn)生。1.5.1 樣條函數(shù)的

12、概念所謂樣條（Spline）本來(lái)是工程設(shè)計(jì)中使用的一種繪圖工具，它是富有彈性的細(xì)木條或細(xì)金屬條。繪圖員利用它把一些已知點(diǎn)連接成一條光滑曲線（稱為樣條曲線），并使連接點(diǎn)處有連續(xù)的曲率。三次樣條插值就是由此抽象出來(lái)的。數(shù)學(xué)上將具有一定光滑性的分段多項(xiàng)式稱為樣條函數(shù)。具體地說(shuō)，給定區(qū)間的一個(gè)分劃 如果函數(shù)滿足： （i）在每個(gè)小區(qū)間上是次多項(xiàng)式；（ii）在上具有階連續(xù)導(dǎo)數(shù)。則稱為關(guān)于分劃的次樣條函數(shù)，其圖形為次樣條曲線。顯然，折線是一次樣條曲線。1.5.2 三次樣條插值利用樣條函數(shù)進(jìn)行插值，即取插值函數(shù)為樣條函數(shù)，稱為樣條插值。例如分段線性插值是一次樣條插值。我們只介紹三次樣條插值，即已知函數(shù)在區(qū)間上

13、的個(gè)節(jié)點(diǎn) 上的值，求插值函數(shù)，使得（i）； （5）（ii）在每個(gè)小區(qū)間上是三次多項(xiàng)式，記為；（iii）在上二階連續(xù)可微。函數(shù)稱為的三次樣條插值函數(shù)。由條件（ii），不妨將記為其中為待定系數(shù)，共個(gè)。由條件（iii） （6）容易看出，（5）、（6）式共含有個(gè)方程，為確定的個(gè)待定參數(shù)，尚需再給出2個(gè)條件。常用的三次樣條函數(shù)的邊界條件有3種類型：（i）。由這種邊界條件建立的樣條插值函數(shù)稱為的完備三次樣條插值函數(shù)。特別地，時(shí)，樣條曲線在端點(diǎn)處呈水平狀態(tài)。如果不知道，我們可以要求與在端點(diǎn)處近似相等。這時(shí)以為節(jié)點(diǎn)作一個(gè)三次Newton插值多項(xiàng)式，以作一個(gè)三次Newton插值多項(xiàng)式，要求由這種邊界條件建立的三

14、次樣條稱為的Lagrange三次樣條插值函數(shù)。（ii）。特別地時(shí)，稱為自然邊界條件。（iii），此條件稱為周期條件。1.5.3 三次樣條插值在Matlab中的實(shí)現(xiàn)在Matlab中數(shù)據(jù)點(diǎn)稱之為斷點(diǎn)。如果三次樣條插值沒(méi)有邊界條件，最常用的方法，就是采用非扭結(jié)（not-a-knot）條件。這個(gè)條件強(qiáng)迫第1個(gè)和第2個(gè)三次多項(xiàng)式的三階導(dǎo)數(shù)相等。對(duì)最后一個(gè)和倒數(shù)第2個(gè)三次多項(xiàng)式也做同樣地處理。Matlab中三次樣條插值也有現(xiàn)成的函數(shù)：y=interp1(x0,y0,x,'spline')；y=spline(x0,y0,x)；pp=csape(x0,y0,conds),pp=csape(x0

15、,y0,conds,valconds)，y=ppval(pp,x)。其中x0,y0是已知數(shù)據(jù)點(diǎn)，x是插值點(diǎn)，y是插值點(diǎn)的函數(shù)值。對(duì)于三次樣條插值，我們提倡使用函數(shù)csape，csape的返回值是pp形式，要求插值點(diǎn)的函數(shù)值，必須調(diào)用函數(shù)ppval。pp=csape(x0,y0)：使用默認(rèn)的邊界條件，即Lagrange邊界條件。pp=csape(x0,y0,conds,valconds)中的conds指定插值的邊界條件，其值可為：'complete' 邊界為一階導(dǎo)數(shù)，一階導(dǎo)數(shù)的值在valconds參數(shù)中給出，若忽略valconds參數(shù)，則按缺省情況處理。'not-a-kn

16、ot' 非扭結(jié)條件 'periodic' 周期條件'second' 邊界為二階導(dǎo)數(shù)，二階導(dǎo)數(shù)的值在valconds參數(shù)中給出，若忽略valconds參數(shù)，二階導(dǎo)數(shù)的缺省值為0, 0。'variational' 設(shè)置邊界的二階導(dǎo)數(shù)值為0,0。對(duì)于一些特殊的邊界條件，可以通過(guò)conds的一個(gè)矩陣來(lái)表示，conds元素的取值為0，1，2。conds(i)=j的含義是給定端點(diǎn)的階導(dǎo)數(shù)，即conds的第一個(gè)元素表示左邊界的條件，第二個(gè)元素表示右邊界的條件，conds=2,1表示左邊界是二階導(dǎo)數(shù)，右邊界是一階導(dǎo)數(shù)，對(duì)應(yīng)的值由valconds給出。詳細(xì)

17、情況請(qǐng)使用幫助help csape。例1 機(jī)床加工待加工零件的外形根據(jù)工藝要求由一組數(shù)據(jù)給出（在平面情況下），用程控銑床加工時(shí)每一刀只能沿方向和方向走非常小的一步，這就需要從已知數(shù)據(jù)得到加工所要求的步長(zhǎng)很小的坐標(biāo)。表中給出的數(shù)據(jù)位于機(jī)翼斷面的下輪廓線上，假設(shè)需要得到坐標(biāo)每改變0.1時(shí)的坐標(biāo)。試完成加工所需數(shù)據(jù)，畫(huà)出曲線，并求出處的曲線斜率和范圍內(nèi)的最小值。 0 3 5 7 9 11 12 13 14 15 0 1.2 1.7 2.0 2.1 2.0 1.8 1.2 1.0 1.6 要求用Lagrange、分段線性和三次樣條三種插值方法計(jì)算。解 編寫(xiě)以下程序：x0=0 3 5 7 9 11 12

18、 13 14 15;y0=0 1.2 1.7 2.0 2.1 2.0 1.8 1.2 1.0 1.6;x=0:0.1:15;y1=lagrange(x0,y0,x); %前面編寫(xiě)的Lagrange插值函數(shù)y2=interp1(x0,y0,x);y3=interp1(x0,y0,x,'spline');pp1=csape(x0,y0);y4=ppval(pp1,x);pp2=csape(x0,y0,'second');y5=ppval(pp2,x);x',y1',y2',y3',y4',y5'subplot(2,2,

19、1)plot(x0,y0,'+',x,y1)title('Lagrange')subplot(2,2,2)plot(x0,y0,'+',x,y2)title('Piecewise linear')subplot(2,2,3)plot(x0,y0,'+',x,y3)title('Spline1')subplot(2,2,4)plot(x0,y0,'+',x,y4)title('Spline2')dx=diff(x);dy=diff(y3);dy_dx=dy./dx;dy

20、_dx0=dy_dx(1)ytemp=y3(131:151);ymin=min(ytemp);index=find(y3=ymin);xmin=x(index);xmin,ymin計(jì)算結(jié)果略?？梢钥闯?，拉格朗日插值的結(jié)果根本不能應(yīng)用，分段線性插值的光滑性較差（特別是在附近彎曲處），建議選用三次樣條插值的結(jié)果。1.6 二維插值前面講述的都是一維插值，即節(jié)點(diǎn)為一維變量，插值函數(shù)是一元函數(shù)（曲線）。若節(jié)點(diǎn)是二維的，插值函數(shù)就是二元函數(shù)，即曲面。如在某區(qū)域測(cè)量了若干點(diǎn)（節(jié)點(diǎn)）的高程（節(jié)點(diǎn)值），為了畫(huà)出較精確的等高線圖，就要先插入更多的點(diǎn)（插值點(diǎn)），計(jì)算這些點(diǎn)的高程（插值）。1.6.1 插值節(jié)點(diǎn)為網(wǎng)格節(jié)

21、點(diǎn)已知個(gè)節(jié)點(diǎn)：（），且；。求點(diǎn)處的插值。Matlab中有一些計(jì)算二維插值的程序。如 z=interp2(x0,y0,z0,x,y,'method')其中x0,y0分別為維和維向量，表示節(jié)點(diǎn)，z0為維矩陣，表示節(jié)點(diǎn)值，x,y為一維數(shù)組，表示插值點(diǎn)，x與y應(yīng)是方向不同的向量，即一個(gè)是行向量，另一個(gè)是列向量，z為矩陣，它的行數(shù)為的維數(shù)，列數(shù)為的維數(shù)，表示得到的插值，'method'的用法同上面的一維插值。如果是三次樣條插值，可以使用命令pp=csape(x0,y0,z0,conds,valconds)，z=fnval(pp,x,y)其中x0,y0分別為維和維向量，z0

22、為維矩陣，為矩陣，它的行數(shù)為的維數(shù)，列數(shù)為的維數(shù)，表示得到的插值，具體使用方法同一維插值。例2 在一丘陵地帶測(cè)量高程，和方向每隔100米測(cè)一個(gè)點(diǎn)，得高程如下表，試插值一曲面，確定合適的模型，并由此找出最高點(diǎn)和該點(diǎn)的高程。 100 200 300 400 500100 636 697 624 478 450 200 698 712 630 478 420300 680 674 598 412 400400 662 626 552 334 310解 編寫(xiě)程序如下：clear,clcx=100:100:500;y=100:100:400;z=636 697 624 478 450 698 712 6

23、30 478 420 680 674 598 412 400 662 626 552 334 310;pp=csape(x,y,z')xi=100:10:500;yi=100:10:400cz1=fnval(pp,xi,yi)cz2=interp2(x,y,z,xi,yi','spline')i,j=find(cz1=max(max(cz1)x=xi(i),y=yi(j),zmax=cz1(i,j)1.6.2 插值節(jié)點(diǎn)為散亂節(jié)點(diǎn)已知個(gè)節(jié)點(diǎn)：，求點(diǎn)處的插值。對(duì)上述問(wèn)題，Matlab中提供了插值函數(shù)griddata，其格式為：ZI = GRIDDATA(X,Y,Z,

24、XI,YI)其中X、Y、Z 均為n維向量，指明所給數(shù)據(jù)點(diǎn)的橫坐標(biāo)、縱坐標(biāo)和豎坐標(biāo)。向量XI、YI是給定的網(wǎng)格點(diǎn)的橫坐標(biāo)和縱坐標(biāo)，返回值ZI為網(wǎng)格（XI，YI）處的函數(shù)值。XI與YI應(yīng)是方向不同的向量，即一個(gè)是行向量，另一個(gè)是列向量。例3 在某海域測(cè)得一些點(diǎn)（x,y）處的水深z由下表給出，在矩形區(qū)域（75,200）×(-50,150) 內(nèi)畫(huà)出海底曲面的圖形。x 129 140 103.5 88 185.5 195 105 157.5 107.5 77 81 162 162 117.5y 7.5 141.5 23 147 22.5 137.5 85.5 6.5 -81 3 56.5 66

25、.5 84 -33.5z 4 8 6 8 6 8 8 9 9 8 8 9 4 9解 編寫(xiě)程序如下：x=129 140 103.5 88 185.5 195 105 157.5 107.5 77 81 162 162 117.5;y=7.5 141.5 23 147 22.5 137.5 85.5 -6.5 -81 3 56.5 -66.5 84 -33.5;z=-4 8 6 8 6 8 8 9 9 8 8 9 4 9;xi=75:1:200;yi=-50:1:150;zi=griddata(x,y,z,xi,yi','cubic')subplot(1,2,1)plot(

26、x,y,'*')subplot(1,2,2)mesh(xi,yi,zi)§2 曲線擬合的線性最小二乘法2.1 線性最小二乘法曲線擬合問(wèn)題的提法是，已知一組（二維）數(shù)據(jù)，即平面上的個(gè)點(diǎn)，互不相同，尋求一個(gè)函數(shù)（曲線），使在某種準(zhǔn)則下與所有數(shù)據(jù)點(diǎn)最為接近，即曲線擬合得最好。線性最小二乘法是解決曲線擬合最常用的方法，基本思路是，令 （7）其中是事先選定的一組線性無(wú)關(guān)的函數(shù)，是待定系數(shù)。擬合準(zhǔn)則是使，與的距離的平方和最小，稱為最小二乘準(zhǔn)則。2.1.1 系數(shù)的確定記 （8）為求使達(dá)到最小，只需利用極值的必要條件，得到關(guān)于的線性方程組即 （9）記，方程組（9）可表為 （10）當(dāng)線

27、性無(wú)關(guān)時(shí)，列滿秩，可逆，于是方程組（10）有唯一解 2.1.2 函數(shù)的選取面對(duì)一組數(shù)據(jù)，用線性最小二乘法作曲線擬合時(shí)，首要的、也是關(guān)鍵的一步是恰當(dāng)?shù)剡x取。如果通過(guò)機(jī)理分析，能夠知道與之間應(yīng)該有什么樣的函數(shù)關(guān)系，則容易確定。若無(wú)法知道與之間的關(guān)系，通?？梢詫?shù)據(jù)作圖，直觀地判斷應(yīng)該用什么樣的曲線去作擬合。人們常用的曲線有（i）直線 （ii）多項(xiàng)式 （一般，不宜太高）（iii）雙曲線（一支） （iv）指數(shù)曲線 對(duì)于指數(shù)曲線，擬合前需作變量代換，化為對(duì)的線性函數(shù)。已知一組數(shù)據(jù)，用什么樣的曲線擬合最好，可以在直觀判斷的基礎(chǔ)上，選幾種曲線分別擬合，然后比較，看哪條曲線的最小二乘指標(biāo)最小。 2.2 最小二

28、乘法的Matlab實(shí)現(xiàn)2.2.1 解方程組方法 在上面的記號(hào)下，Matlab中的線性最小二乘的標(biāo)準(zhǔn)型為 命令為 。例4 用最小二乘法求一個(gè)形如的經(jīng)驗(yàn)公式，使它與下表所示的數(shù)據(jù)擬合。 19 25 31 38 44 19.0 32.3 49.0 73.3 97.8解 編寫(xiě)程序如下x=19 25 31 38 44'y=19.0 32.3 49.0 73.3 97.8'r=ones(5,1),x.2;ab=ryx0=19:0.1:44;y0=ab(1)+ab(2)*x0.2;plot(x,y,'o',x0,y0,'r')2.2.2 多項(xiàng)式擬合方法如果取，

29、即用次多項(xiàng)式擬合給定數(shù)據(jù)，Matlab中有現(xiàn)成的函數(shù)a=polyfit(x0,y0,m)其中輸入?yún)?shù)x0,y0為要擬合的數(shù)據(jù)，m為擬合多項(xiàng)式的次數(shù)，輸出參數(shù)a為擬合多項(xiàng)式y(tǒng)=amxm+a1x+a0系數(shù)a= am, , a1, a0。多項(xiàng)式在x處的值y可用下面的函數(shù)計(jì)算y=polyval(a,x)。例5 某鄉(xiāng)鎮(zhèn)企業(yè)1990-1996年的生產(chǎn)利潤(rùn)如下表：年份 1990 1991 1992 1993 1994 1995 1996利潤(rùn)（萬(wàn)元） 70 122 144 152 174 196 202試預(yù)測(cè)1997年和1998年的利潤(rùn)。解 作已知數(shù)據(jù)的的散點(diǎn)圖，x0=1990 1991 1992 1993

30、1994 1995 1996;y0=70 122 144 152 174 196 202;plot(x0,y0,'*')發(fā)現(xiàn)該鄉(xiāng)鎮(zhèn)企業(yè)的年生產(chǎn)利潤(rùn)幾乎直線上升。因此，我們可以用作為擬合函數(shù)來(lái)預(yù)測(cè)該鄉(xiāng)鎮(zhèn)企業(yè)未來(lái)的年利潤(rùn)。編寫(xiě)程序如下：x0=1990 1991 1992 1993 1994 1995 1996;y0=70 122 144 152 174 196 202;a=polyfit(x0,y0,1)y97=polyval(a,1997)y98=polyval(a,1998)求得，1997年的生產(chǎn)利潤(rùn)y97=233.4286，1998年的生產(chǎn)利潤(rùn)y98=253.9286。

31、67;3 最小二乘優(yōu)化在無(wú)約束最優(yōu)化問(wèn)題中，有些重要的特殊情形，比如目標(biāo)函數(shù)由若干個(gè)函數(shù)的平方和構(gòu)成。這類函數(shù)一般可以寫(xiě)成： ,其中，一般假設(shè)。我們把極小化這類函數(shù)的問(wèn)題：稱為最小二乘優(yōu)化問(wèn)題。最小二乘優(yōu)化是一類比較特殊的優(yōu)化問(wèn)題，在處理這類問(wèn)題時(shí)，Matlab也提供了一些強(qiáng)大的函數(shù)。在Matlab優(yōu)化工具箱中，用于求解最小二乘優(yōu)化問(wèn)題的函數(shù)有：lsqlin、lsqcurvefit、lsqnonlin、lsqnonneg，用法介紹如下。3.1 lsqlin 函數(shù)求解 s.t. 其中為矩陣，為向量。Matlab中的函數(shù)為：X=LSQLIN(C,d,A,b,Aeq,beq,LB,UB,X0)例6

32、用lsqlin命令求解例4。解 編寫(xiě)程序如下：x=19 25 31 38 44'y=19.0 32.3 49.0 73.3 97.8'r=ones(5,1),x.2;ab=lsqlin(r,y)x0=19:0.1:44;y0=ab(1)+ab(2)*x0.2;plot(x,y,'o',x0,y0,'r')3.2 lsqcurvefit 函數(shù)給定輸入輸出數(shù)列,求參量,使得 Matlab 中的函數(shù)為X=LSQCURVEFIT(FUN,X0,XDATA,YDATA,LB,UB,OPTIONS)其中FUN是定義函數(shù)的M文件。例7 用下面一組數(shù)據(jù)擬合函數(shù)中

33、的參數(shù)。 100 200 300 400 500 600 700 800 900 1000 4.54 4.99 5.35 5.65 5.90 6.10 6.26 6.39 6.50 6.59解 該問(wèn)題即解最優(yōu)化問(wèn)題： （1）編寫(xiě)M文件fun1.m定義函數(shù)： function f=fun1(x,tdata);f=x(1)+x(2)*exp(-0.02*x(3)*tdata); %其中x(1)=a，x(2)=b，x(3)=k（2）調(diào)用函數(shù)lsqcurvefit，編寫(xiě)程序如下：td=100:100:1000;cd=4.54 4.99 5.35 5.65 5.90 6.10 6.26 6.39 6.5

34、0 6.59;x0=0.2 0.05 0.05;x=lsqcurvefit(fun1,x0,td,cd)3.3 lsqnonlin 函數(shù)已知函數(shù)向量，求使得 Matlab中的函數(shù)為X=LSQNONLIN(FUN,X0,LB,UB,OPTIONS)其中FUN是定義向量函數(shù)的M文件。例8 用lsqnonlin 函數(shù)求解例7。解 這里 （1）編寫(xiě)M文件fun2.m如下：function f=fun2(x);td=100:100:1000;cd=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)*td

35、)-cd;（2）調(diào)用函數(shù)lsqnonlin，編寫(xiě)程序如下：x0=0.2 0.05 0.05;x=lsqnonlin(fun2,x0)3.4 lsqnonneg 函數(shù)求解非負(fù)的，使得滿足 ，Matlab中的函數(shù)為X = LSQNONNEG(C,d,X0,OPTIONS)例9 已知，求滿足 最小。解 編寫(xiě)程序如下：c=0.0372 0.2869;0.6861 0.7071;0.6233 0.6245;0.6344 0.6170;d=0.8587;0.1781;0.0747;0.8405;x=lsqnonneg(c,d)§4 曲線擬合與函數(shù)逼近前面講的曲線擬合是已知一組離散數(shù)據(jù)，選擇一個(gè)較

36、簡(jiǎn)單的函數(shù)，如多項(xiàng)式，在一定準(zhǔn)則如最小二乘準(zhǔn)則下，最接近這些數(shù)據(jù)。如果已知一個(gè)較為復(fù)雜的連續(xù)函數(shù)，要求選擇一個(gè)較簡(jiǎn)單的函數(shù)，在一定準(zhǔn)則下最接近，就是所謂函數(shù)逼近。與曲線擬合的最小二乘準(zhǔn)則相對(duì)應(yīng)，函數(shù)逼近常用的一種準(zhǔn)則是最小平方逼近，即 (11)達(dá)到最小。與曲線擬合一樣，選一組函數(shù)構(gòu)造，即令代入（11）式，求使達(dá)到極小。利用極值必要條件可得 （12）這里。當(dāng)方程組（12）的系數(shù)矩陣非奇異時(shí)，有唯一解。最簡(jiǎn)單的當(dāng)然是用多項(xiàng)式逼近函數(shù)，即選，。并且如果能使，方程組（12）的系數(shù)矩陣將是對(duì)角陣，計(jì)算大大簡(jiǎn)化。滿足這種性質(zhì)的多項(xiàng)式稱正交多項(xiàng)式。勒讓得（Legendre）多項(xiàng)式是在區(qū)間上的正交多項(xiàng)式，它的表達(dá)式為可以證明常用的正交多項(xiàng)式還有第一類切比雪夫（Chebyshev）多項(xiàng)式和拉蓋爾（Laguerre）多項(xiàng)式 。例10 求在中的最佳平方逼近多項(xiàng)式。 解 編寫(xiě)程序如下：syms xbase=1,x2,x4;y1=base.'*basey2=cos(x)*base.'r1=int(y1,-pi/2,pi/2)r2=int(y2,-pi/2,pi/2)a=r1r2xishu1=double(a)digits(8),xishu2=vpa(a)求得xishu1=0.9996 -0.4964 0.0372，即所求的最佳平方逼近多項(xiàng)式為.習(xí) 題 九1.