數(shù)值分析課程報告

插值法和多項式擬合的研究摘要在科研和生產實踐中，常常需要通過一組測量數(shù)據(jù)來尋找變量X與y的函數(shù)關系近似表達式。解決這類問題的方法有兩種：一種是插值法，另一種是擬合法。插值法的原理是用一個簡單函數(shù)逼近被計算函數(shù)，然后用該簡單函數(shù)的函數(shù)值近似替代被計算函數(shù)的函數(shù)值。擬合法能夠是從給定的一組實驗數(shù)據(jù)出發(fā)，尋找函數(shù)的一個近似表達式，該近似表達式能反映數(shù)據(jù)的基本趨勢而又不一定過全部的點，即曲線擬合。本文主要介紹拉格朗日插值法、埃爾米特插值法、三次樣條插值法以及基于最小二乘法的多項式擬4□o關鍵詞：拉格朗日插值，埃爾米特插值，樣條插值，多項式擬合1方法的意義在許多實際問題及科學研究中，因素之間往往存在著函數(shù)關系，然而，這種關系經(jīng)常很難有明顯的解析表達，通常只是由觀察與測試得到一些離散數(shù)值。有時，即使給出了解析表達式，卻由于表達式過于復雜，不僅使用不便，而且不易于進行計算與理論分析。解決這類問題的方法有兩種：一種是插值法，另一種是擬合法。插值法的原理是用一個簡單函數(shù)逼近被計算函數(shù)，然后用該簡單函數(shù)的函數(shù)值近似替代被計算函數(shù)的函數(shù)值。它要求給出函數(shù)的一個函數(shù)表，然后選定一種簡單的函數(shù)形式，比如多項式、分段線性函數(shù)及三角多項式等，通過已知的函數(shù)表來確定一個簡單的函數(shù)(x)作為f(x)的近似，概括地說，就是用簡單函數(shù)為離散數(shù)組建立連續(xù)模型。插值法在實際應用中非常廣泛，但是它也有明顯的缺陷，一是測量數(shù)據(jù)常常帶有測試誤差，而插值多項式又通過所有給出的點，這樣就是插值多項式保留了這些誤差；二是如果實際得到的數(shù)據(jù)過多，則必然得到次數(shù)較高的插值多項式，這樣近似的效果并不理想。擬合法能夠很好的解決這些問題，它從給定的一組實驗數(shù)據(jù)出發(fā)，尋找函數(shù)的一個近似表達式y(tǒng)=「(x),該近似表達式能反映數(shù)據(jù)的基本趨勢而又不一定過全部的點，即曲線擬合的問題，函數(shù)的近似表達式y(tǒng)=“)稱為擬合曲線。常用最小而二乘法來確定擬合曲線。2插值法的介紹2.1插值法定義設f(x)為［a，b］上的函數(shù)，在互異點xo*…,Xn處的函數(shù)值分別為f(Xo),f(xj,...,f(Xn),構造一個簡單函數(shù)(x)作為函數(shù)f(x)的近似表達式y(tǒng)=f(x)」(x)，使g二f(XJ,i=0,1,2,...,n(1.0)則稱(x)為關于節(jié)點X0,X-!5...,Xn的插值函數(shù)；稱x0,x1,...,xn為插值節(jié)點；稱(Xif(Xi)),i=1,2，…，n為插值點；f(x)稱為被插值函數(shù)。式(1.0)稱為插值條件。這類問題稱為插值問題。插值的任務就是由已知的觀測點，為物理量(未知量)建立一個簡單的、連續(xù)的解析模型，以便能根據(jù)該模型推測該物理量在非觀測點處的特性。常用的插值函數(shù)類｛「(X)｝是代數(shù)多項式，相應插值問題是代數(shù)插值，本文主要介紹三種代數(shù)差值：拉格朗日插值，埃爾米特插值和樣條插值。2.2拉格朗日插值2.2.1兩點插值問題已知y=f(x)(=0,1)，求滿足插值條件R(x)=y(i=o,i)的插值多項式R(x)。根據(jù)解析幾何知識可知，所求的R(x)為過點區(qū)鳳),*乂)的直線，即：R(x)=yo+—(x-xo)i—xox-Xo顯然lo(Xi)=(i=O)Xi-Xoio(1=1)1(i=1)上式經(jīng)整理可改寫為：式中l(wèi)o(x)二li(x)且lo(X)，h(X)為由插值節(jié)點唯一確定的線性函數(shù)。lo(x)，I1(X)為節(jié)點Xo,X1，上的一次插值基函數(shù)。可以看出，節(jié)點Xo,X1上的插值基函數(shù)的次數(shù)為插值節(jié)點個數(shù)減一，基函數(shù)組中所含的函數(shù)個數(shù)與插值節(jié)點數(shù)相同。而滿足R(Xi)=yi(i=O,1)的插值多項式P|(x)就是節(jié)點xo，x1上插值基函數(shù)的線性組合，其組合系數(shù)分別為yo,y1。這種表示為插值基函數(shù)線性組合的一次插值多項式也就是一次拉格朗日插值多項式。當給定n+1個插值節(jié)點后，可類似定義n次插值基函數(shù)，并以此構造n次拉格朗日多項式。2.2.2n次拉格朗日插值多項式n次插值基函數(shù)：(k=O,1...n)l(x)-(X_Xo)■■-(X_XkJ(X—Xk—Xn)k_(Xk—Xo)???(XkXkJ(Xk-Xk1)..?"*)顯然lk(x)具有以下性質：性質1，lk(x)」1a—*)(k=0,1…，n)kO(i=k)性質2,lk(x)(k=0,1…，n)為由插值節(jié)點Xo,X1,…,xn唯一確定的n次函數(shù)(k=O,1...n)以n次插值基函數(shù)為基礎，可得拉格朗日插值多項式為:L(x)-l/x)\，L(X1偵版1Xk」)(x1Xk"...(x1Xn)yykl(X)lk(x)yk、nk￡7(Xk_，..(Xk_Xk/)(Xkykn.n..八(|【上乞皿k衛(wèi)i」兀-Xji#nn若記「nl(X)(X-Xj),則有51(XQ八(兀-人)。于是-(X)也可寫成：7怫5nl(X)心5—X"k)223拉格朗日插值多項式的余項Rn(X)=f(X)-Ln(X)稱為拉格朗日插值多項式的余項，也叫截斷誤差。用簡單的插值多項式Ln(x)代替復雜的函數(shù)f(x)，這種做法是否有效，取決于截斷誤差是否滿足精度要求。拉格朗日型余項：Rn(X)二mW")")其中^1&)「(x-X"'-(a,b)且依賴于Xo。一般情況下，匚e(a,b)的具體i=Q數(shù)值無法知道，但是若能夠求出maxf(n1)(x)二Mni，則可以得出插值多項式的截斷誤差限為：Rn(x)(nRn(x)(n1)33i(x)n1由此看出，Rn(x)的大小除了與Mm有關外還與插值節(jié)點有密切關系。當::m)個作為插值條件而求某點給定m個點出的函數(shù)值，但僅選用其中n1(n::m)個作為插值條件而求某點x處函數(shù)值時，n+1個節(jié)點x0，x1，…，xn的選取應該盡可能的接近x,使得計算的函數(shù)值的誤差限盡可能的小。許多實際問題不但要求插值多項式與被插值函數(shù)在節(jié)點處的函數(shù)值相當，且還需要求其導數(shù)值相等。滿足這種要求的插值多項式就是埃爾米特插值多項式。一般情形的埃爾米特插值問題一般情形的埃爾米特插值問題是指所滿足的插值條件中函數(shù)值的個數(shù)與導數(shù)值的個數(shù)相等。即當函數(shù)f(x)在區(qū)間】a,b]上n+1個節(jié)點X.(i=0,1，…,n)處的函數(shù)值V、二f(xi)及導數(shù)值mAf(xi)給定時，要求一個次數(shù)不超過2n+1的多項式H2ni(x)，使之滿足：0坤(i="n)H2n"J二m這里給出個2n+2個插值條件，可唯一確定一個形式為H2ni(x)二ao?3iX?a2niX2」1的多項式，但是要確定2n+2個系數(shù)非常復雜，因此此時可以借用構造拉格朗日插值多項式的基函數(shù)。設:“(xh/x)(j-0,1,...,n)為次數(shù)不超過2n+1的多項式，且滿足:一：仲)='『匚、㈤=。(j,j=0」，…,n)(XA0,餌■j-j(G-j則滿足上述條件的埃爾米特插值多項式可以寫成用插值基函數(shù)表示的形式：nH2n1j(x)Vi：j(x)mj]j￡上式滿足插值條件?？梢缘玫交瘮?shù):j(x),'j(x)的解析式為：-T—2(x—*)為*l2(x)(j=0,1,5,Ik=0.兀Pj(x)=(x—Xj)lj2(x)(j=0,1,..n),式中1何)為拉格朗日插值基函數(shù)因此埃爾米特插值多項式即為：nH2n〔兇八1—2(x—Xj)Elj2(x)yi、(x-X.)l.(x2m.11」4V|/|V/Ij=0j=般情形下的埃爾米特插值多項式的余項為:(2n+2)!式中：；：=(a,b),且與x有關。埃爾米特插值的幾何意義：曲線y二tn1(x)與曲線y=f(x)在插值節(jié)點處有相同的公共切線。在帶導數(shù)的插值問題中，有時插值條件中的函數(shù)值個數(shù)與導數(shù)值個數(shù)不等。這時可以以一般情況的埃爾米特插值多項式為基礎，運用待定系數(shù)法求出滿足插值條件的多項式。2.4樣條插值在上述方法中，我們根據(jù)區(qū)間［a,b］上給出的節(jié)點得以得到函數(shù)f(x)的插值多項式，但是并非插值多項式的次數(shù)越高，逼近函數(shù)f(x)的精度越好，主要原因是因為對于任意的插值節(jié)點，當n》時，插值多項式Pn(x)不一定收斂到f(x)。這種高次插值不準確的現(xiàn)象稱為龍格現(xiàn)象。為了避免高次插值的缺點，人們常常采用分段插值的方法，即將插值區(qū)間分為若干個小區(qū)間，在每個小區(qū)間上運用前面介紹的插值方法構造低次插值多項式。采用分段現(xiàn)象插值與分段二次插值，可以構造一個整個連續(xù)的函數(shù)，而采用分段三次埃爾米特插值則可以構造一個整體上具有一節(jié)連續(xù)導數(shù)的插值函數(shù)。實際問題中，很少給出插值點上的導數(shù)值。三次樣條插值就是在只給出插值點上函數(shù)值的情況下，構造一個整體上具有二階連續(xù)導數(shù)的插值函數(shù)。2.4.1三次樣條插值函數(shù)的定義設區(qū)間［a,b］上取n+1個節(jié)點a=X。；音：：：x2::...：：：xnJ<x,d=b,若函數(shù)S(x)滿足條件：在整個區(qū)間［a,b］上具有二階連續(xù)導數(shù)在每個小區(qū)間［Xy,x］(i=0,1，…,n)上是x的三次多項式S(x)=yi(i=0,1,?-n),則稱S(x)為f(x)的三次樣條插值。S(x)的邊界條件為：給定兩端點處的一階導數(shù)值，記為：S(Xo)=m°S(Xn)=mn給定兩端點處的二階導數(shù)值，記為：S(X°)=M。，s*)二M此外，對S(Xo)=S(Xn)=0的邊界條件稱為自然邊界條件2.4.2三次樣條插值函數(shù)的構造構造三次樣條插值函數(shù)，就是要寫出它在子區(qū)間［x」xj(i=0,1,...,n)上的表達式，記為S(x)(i=0,1，…,n)。一、用節(jié)點處一階導數(shù)表示的三次樣條插值函數(shù)記節(jié)點處的一階導數(shù)值為S(x)二mi(i=0,1,...,n)，若已知叫后，則S(x)在［X」X］(i=0,1，…,n)上就是滿足條件S(x」)」，S(x)=y，(i=0,1,...,n)i=yiiiS(x4)=m」，S(x)=mj(i=0,1,…,n)的三次埃爾米特插值多項式。步驟如下：根據(jù)S(x)，S(x)在內階段的連續(xù)性及插值條件，運用［Xj^Xj］上的二點三次埃爾米特插值多項式，寫出S(x)用mi(i=0,1,...,口)表示的形式。利用S(x)在內節(jié)點Xi(i=0,1,...,n-1)的連續(xù)性及邊界條件，導出含mi(i=0,1,...,n)的n+1階線性方程組。⑶求解含mi(i=0,1,...,n)的線性方程組，將得到的m：代入伙二必］的二點三次埃爾米特插值多項式具體公式為:+仆+仆2x_IX—X.LXi頑S.(x)=1+24SX—XYi4x—Xij八X、Xi-Xi+(X-X.)x一Xi!mu+(x—xjX—4mi1g—XJXi_Xi也m.的確定由S(x)得二階導數(shù)在內節(jié)點連續(xù)的條件得到，即：.m|42m.Mm1二(i=0,1,..n,1)該式即為n+1階線性方程組。其中:hi=Xi一x」hhM1i=1-hhM1i=1-■」=3(Mjh(的,1,..nr1)1mo2M1A22M2m2+++2M/mnjAn4」m勺經(jīng))州申12f2扎h.i4in-1求解該線性方程組還需兩個邊界條件。對于第一種邊界條件，將其代入，得到只含n-1個未知數(shù)的線性方程組:■2M1*-22M2+■2M1*-22M2+++》-n2對于第二種邊界條件，將其代入，其中：4n+1得到n_22m1、m2階線性方程組:/rr,、f2*lfn』-Mn/mnJ卜3坐必-九0h12fn=3上歸hn2二、用節(jié)點處二階導數(shù)表示的三次樣條插值函數(shù)記節(jié)點處的二階導數(shù)值為S(xJ=mj(i=0,1,...,n)，由于S(x)在［八4,xJ(i=0,1，…,口)上是x的線性函數(shù)，因此構造以節(jié)點處二階導數(shù)表示的三次埃爾米特插值多項式步驟如下:根據(jù)S(x)在內節(jié)點的連續(xù)性及為線性函數(shù)的特點，將S伙)表示為線性函數(shù)。在根據(jù)S(x)在內節(jié)點的連續(xù)性及插值條件，寫出S(x)用M.(i=0,1,...,可表示的形式。⑵利用S(x)在內節(jié)點xi(i=0,1,...,n-1)的連續(xù)性及邊界條件，導出

求解含Mji=0,1,...,n)的線性方程組，將得到的M/代入［Xj,Xj的S(x)的表達式，即得到二點三次埃爾米特插值多項式具體公式為：Si(x)=AAMijAXhAMi(%」MJ十、為一xMix-6hi6hi-rh)〒(yrxjhMi的確定由S(x)得二階導數(shù)在內節(jié)點連續(xù)的條件得到，即：MMij2M「：；’％Mi4=￡(i=0,1,..n,1)該式即為n+1階線性方程組。其中:hi=x1(i71,..n廠1)6hihi(i71,..n廠1)6hihi1(y［胡i1-yihi1解該方程還需要兩個邊界條件。對于第一種邊界條件，將其代入，得到只含_21【M12人+++Mn42打4式中:n-1個未知數(shù)的線性方程組:rm1irf0；(iiiM2Ifl?I■I1rpiM口fJj)Mn，X-■xy.jf6(yohihi-mo)對于第二種邊界條件，將其代入6/(mn-hn得到-yn4)hnn+1階線性方程組:，ZM1'm2

.Mn-22Mn4-嚴對于三次樣條插值函數(shù)來說，當插值節(jié)點逐漸加密時，M9n2AM、nJFfM1M0f2afn.(fn1-扎n/Mn可以證明：不但樣條插值函數(shù)收斂于函數(shù)本身，而且其導數(shù)也收斂于函數(shù)的導數(shù)。3.曲線擬合3.1用最小二乘法求解矛盾方程組工程實際中許多問題都可歸結為矛盾方程組，實際中需要尋求矛盾方程組的N未知數(shù)Xj,X2,…,xn的一組取值，它使得偏差的絕對值之和送問盡可能的小，為了7NNn2便于分析計算和應用，常采用使偏差的平方和Q=vra.X._bJ達到i=4j=4最小值，這一條件稱為最小二乘原則。按最小二乘原則來選擇X4，X2，...,Xn的一組取值的方法稱為求解矛盾方程組的最小二乘法，符合條件的Xj,X2,…,Xn的一組取值稱為矛盾方程組的最小二乘解。把Q看出n個自變量Xj,X2，…,Xn的二次函數(shù)，記為Q=f(Xj,X2，…,Xn),因此n求矛盾方程組的最小二乘解就是二乘函數(shù)的最小值點。若矛盾方程組-BjXAqj呂的系數(shù)矩陣A的秩為幾則二次函數(shù)Q=f(Xj，X2，…,Xn)一定存在最小值。即矛盾方程組的最小二乘解存在，且正則方程組AtAx二ATb有唯一解，此解就是矛盾方程組的最小二乘解。3.2多項式擬合設通過測量得到函數(shù)y二f(x)的一組數(shù)據(jù)為：y=HXJgO」,...，n),求一個次數(shù)低于N-1的多項式y(tǒng)=(x)=a0?a4X-a2X2...■amXm(m:::N-1)(其中NN2a0,a4,…,am的解a0,a4,…,amNN2a0,a4,…,am的解偏差「二「(為r的平方和Q八f「(N-yj］達到最小。i#.二4ii#即正則方程組AtAx二ATb的解。lNNNNTOC\o"1-5"\h\za°N+a.】天Xj+&2瓦x?+...+amZ小i-1i」\-1\-1NNNNN』a°送X.七1送x2+aAx;+...+am瓦x.*=送xj%i￡i）i」i=i-aNNNNNa。瓦錄七正xm++aAxm亡+...+amEx：m匹X.、-yi￡777由該正則方程組求得的唯一解代入擬合多項式y(tǒng)=「（x），即為所求。5.數(shù)值試驗表1是1971年到1990年我國總人口的統(tǒng)計數(shù)字，試根據(jù)1971年到1985年這15年人口統(tǒng)計數(shù)字用下面幾種方法預測未來20年的人口數(shù)字，并用圖示的方法比較1986年到1990年間預測人口數(shù)字與實際統(tǒng)計數(shù)字的差異，在你所使用的幾種預測方法中找出一組較為合理的預測方法。(1)指數(shù)形式y(tǒng)二aebx；拉格朗日插值、埃爾米特插值、樣條插值；三次多項式擬合；四次多項式擬合表1人口統(tǒng)計數(shù)字年份1971197219731974197519761977197819791980人口8.52298.71778.92119.08599.24209.37179.49749.62599.75429.8705年份1981198219831984198519861987198819891990人口10.007210.165410.300810.435710.585110.750710.930011.102611.270411.4333數(shù)值試驗結果:1.指數(shù)形式擬合圖形為:三次多項式擬合的效果:199019臘實際統(tǒng)計曲線實際統(tǒng)計值三次多項式擬合19861984199019臘19801978197619741972■4*aIIIIIIIIIIItIIil|197*8b1930197199019881986819761930197199019881986819761971972四初多項式擬合結論：基于最小二乘法的多項式擬合能取得良好的函數(shù)逼近效果。階數(shù)越高，擬合效果越好。參考文獻[1]王正林，龔純，等.精通MATLAB科學計算.北京：電子工業(yè)出版社，2009:135—157任玉杰.數(shù)值分析及其MATLAB實現(xiàn)]M].北京：高等教育出版社，2007:584—642李信莫車剛明，等.計算方法[M].西安：西北工業(yè)大學出版社，2010:150—167金一慶，陳越.數(shù)值方法[M].北京：機械工業(yè)出版社，2000:165—176附錄1■按指數(shù)形式擬合的程序：x=[8.5229,8.7177,8.9211,9.0859,9.2420,9.3717,9.4974,9.6259,9.7542,9.8705,10.0072,10.1654,10.3008,10.4357,10.5851,10.7507,10.9300,1;y=[1971,1972,1973,1974,1975,1976,1977,1978,1979,1980,1981,1982,19;p=polyfit(x,