插值法與最小二乘法課件

插值法與最小二乘法

在實際問題中遇到的函數(shù)有些有解析表達式，但很復(fù)雜，有些只給出一些離散數(shù)據(jù)，這給我們求解函數(shù)值、導(dǎo)數(shù)值、零點、極值和積分值等帶來了諸多不便。對于這些情況，自然的想法是，設(shè)法找到某個簡單近似函數(shù)滿足。本章介紹兩種方法，即插值法和最小二乘法。

§3.1插值法3.1.1多項式的插值概念在眾多的插值函數(shù)中，多項式是最簡單最易計算的。設(shè)函數(shù)在區(qū)間[a,b]上連續(xù)，且在n+1個不同的點上分別取值為。在多項式插值中，最基本、最簡單的問題是求一個次數(shù)不超過

的代數(shù)多項式

（3-1）滿足其中均為實數(shù)（3-2）稱為被插值函數(shù)；稱插值多項式；條件（3-2）為插值條件；為插值點。其幾何意義為圖3-1插值幾何意義滿足（3-2）式的插值多項式是存在且唯一的。原因是滿足（3-2）式的多項式的未知系數(shù)行列式為著名的范德蒙德（Vardermonde）行列式故解是存在且唯一的。3.1.2拉格朗日插值多項式

在每個插值點構(gòu)造插值基函數(shù)為一個

次多項式，且滿足條件

（3-3）（3-4）由（3-4）式和多項式的定義以及（3-2）式的插值條件，我們有：再由（3-4）式知，是的零點，故按因式定理必含有因式：既是次多項式，可知

再由條件，得于是有因此拉格朗日插值多項式可寫為

（3-5）例3-1已知特殊角，用

的近似值。解：令用為節(jié)點，有一次和二次拉格朗日插值公式求用為節(jié)點，有拉格朗日插值多項式程序設(shè)計計算公式(3-5)的程序為二重循環(huán)。由內(nèi)循環(huán)（j循環(huán)）通過累乘求得基函數(shù)，然后通過外循環(huán)（i循環(huán)）得到插值結(jié)果y.下圖描述了拉格朗日算法流程，對于該框圖中的有關(guān)參數(shù)說明如下：N:

插值多項式次數(shù)（插值點個數(shù)-1）Xi：插值點Yi:

插值點上的函數(shù)值LN:拉格朗日插值結(jié)果T:

存放累乘積I:

外循環(huán)變量J:

內(nèi)循環(huán)變量STARTINPUTN,X,Xi,Yi(I=0,…N)LN=0I=0,NT=1J=0,NJ=I?NOT=(X-Xj)/(Xi-Xj)*TENDJLN=LN+Yi*TENDIOUTPUTENDYES拉格朗日算法程序DIMENSIONX(I),Y(I)REALLNREAD(*,*)NDOK=1,NREAD(*,*)X(I),Y(I)ENDDOLN=0.0DOI=0,NT=1DOJ=0,NT=T*(X-X(J))/(X(I)-X(J))ENDDOLN=LN+T*Y(I)ENDDOWRITE(*,*)‘LN=‘,LN1、插值余項對插值后，用代替必定會產(chǎn)生誤差，其誤差可表示為

稱上式為插值多項式

的插值余項。

設(shè)在區(qū)間上有直到階導(dǎo)數(shù)，為

上個互異的節(jié)點，為滿足條件的

次插值多項式。

那么對于任何，有

(3-6)其中且依賴于。（一般取最大值

3.1.3插值多項式余項3.1.4牛頓插值多項式在插值基函數(shù)中，按這種方式來求：由線性代數(shù)性質(zhì)可知，一個不高于n次的多項式可表示成的線性組合。在滿足插值條件情況下，可表示為（記為）：

（3-6）上式稱為牛頓（Newton）插值多項式。

3.1.4牛頓插值多項式1.1、向前差分與牛頓向前插值多項式

取節(jié)點為等距離，則

其中h為步長。在兩相臨節(jié)點處的函數(shù)值之差為

稱為函數(shù)f(x)在節(jié)點處以h為步長的一階向前差分（簡稱一階差分）。又稱一階差分的差分為二階差分，記為

3.1.4牛頓插值多項式為了便于計算與一目了然，用表格描述為

表3.13.1.4牛頓插值多項式在等距離節(jié)點的情況下，可以用差分來確定(3.6)式的待定系數(shù)，并將牛頓插值多項式加以簡化。

在節(jié)點上，函數(shù)值是已知的，所以

3.1.4牛頓插值多項式將上式帶入(3.6)式，并令，則式(3.6)可簡化為

該式稱為牛頓向前插值公式。其插值余項是

（3-7）（3-8）3.1.4牛頓插值多項式2、向后差分與牛頓向后插值公式a、向后差分b、中心差分

3.1.4牛頓插值多項式

3、差商與牛頓插值多項式用差商來確定(3.6)式中的系數(shù)一階差商二階差商

依次類推

3.1.4牛頓插值多項式表3.２差商表3.1.4牛頓插值多項式x4.00024.01044.02334.0294f(x)0.60208170.60318770.60458240.6052404例3-2

給出列表函數(shù)f(x)=lgx的值，如下表所示，試用牛頓插值法求lg4.01的值。表3.3解：根據(jù)給定的函數(shù)作差商表，如表3.3所示。

x=4.01,x1=4.0002,x2=4.0233,x3=4.02943.1.4牛頓插值多項式lgx=f(x0)+f(x0,x1)(x-x0)+f(x0,x1,x2)(x-x0)(x-x1)+f(x0,x1,x2,x3)(x-x0)(x-x1)(x-x2)表3.4x4.00024.01044.02334.02940.60208170.60318770.60458240.60524040.1084330.1081160.107869-0.013636-0.0130000.021781f(x)一階差商二階差商三階差商將以上各插值點值、函數(shù)值和差商值代入插值公式得lg4.01=0.60314433.1.5

Hermite插值多項式

前述的插值多項式均只要求插值點函數(shù)值相等，為進一步光滑函數(shù)曲線，有時還要求其導(dǎo)數(shù)值也相等，并等于已知值的次多項式。該多項式稱為的Hermite插值多項式，它也是一種近似式。求解Hermite插值多項式的一種最簡單方法就是直接應(yīng)用牛頓插值法。即已知節(jié)點處函數(shù)值與導(dǎo)數(shù)值時，把視為重節(jié)點，并注意個重節(jié)點的階差商3.1.5

Hermite插值多項式例3-2已知函數(shù)

的函數(shù)值、導(dǎo)數(shù)值如下，求其插值多項式及誤差表達式。-1010-4-2056解：3.1.5

Hermite插值多項式作差商表如下。f(0,0)=y’(0)/1!=0,f(0,0,0)=y”(0)/2!=3,f(1,1)=y’(1)/1!=5.按差商表3.2有：-1000110-4-4-4-2-2-400254323-1-110213.1.6三次樣條插值多項式1、三次樣條插值函數(shù)的定義對于一組已知的數(shù)據(jù)，且，若函數(shù)滿足：

(1)、在每一個子區(qū)間上都可以用最高為三次的多項式來表示。⑵、在上的二次連續(xù)可微。即，，連續(xù)。

⑶、插值條件則稱為函數(shù)f(x)關(guān)于節(jié)點的三次樣條插值函數(shù)。

3.1.6三次樣條插值多項式

2、邊界條件問題根據(jù)定義，在任意子區(qū)間上，三次樣條函數(shù)可表示為

該式中有4個待定系數(shù)，故共有4n個待定系數(shù)。由定義中的連續(xù)條件和插值條件可列出下列方程：

3.1.6三次樣條插值多項式

上式共有4n-2個方程，因此還需要2個方程才能確定4n個系數(shù)，這就要用到邊界條件。邊界條件的類型很多，常見的有三種：

[1]、在邊界上給定一階導(dǎo)數(shù)的值；

[2]、給定二階導(dǎo)數(shù)的值；（如果，稱其為自然邊界條件

[3]、f(x)是以b-a為周期的函數(shù)，則要求S(x)及其導(dǎo)數(shù)也是以b-a為周期的函數(shù)，相應(yīng)邊界條件為

3.1.6三次樣條插值多項式3、三彎矩方程設(shè)S(x)在某一節(jié)點的二階導(dǎo)數(shù)為

（在力學(xué)上解釋為：Mi是細梁在截面處的彎矩，因為該彎矩與相臨的兩個彎矩有關(guān)，故稱三彎矩方程）因為S(x)最高次數(shù)不超過三次，所以它的二階導(dǎo)數(shù)是線性函數(shù)。令，于是

3.1.6三次樣條插值多項式連續(xù)積分兩次得其中為積分常數(shù)。由以下插值條件確定

所以

（3-9）3.1.6三次樣條插值多項式代入(3.9)式得從上式可知，只要確定M值（M0,M1,…,Mn

共n+1個），就可以完全確定S(x)。在確定(3.9)式中待定系數(shù)時用了插值點的條件，下面利用節(jié)點一階導(dǎo)數(shù)連續(xù)條件（3-10）（3-11）3.1.6三次樣條插值多項式由(3.10)式得于是（3-12）（3-13）3.1.6三次樣條插值多項式由(3.11)式、(3.12)式和(3.13)式得若記則方程可簡寫為

（3-14）3.1.6三次樣條插值多項式即上方程組有n-1個方程，但有n+1個未知數(shù)。因此，要確定Mi還需用邊界條件。

[1]、第一種邊界問題邊界條件則（3-15）3.1.6三次樣條插值多項式整理得

其中

（3-16）3.1.6三次樣條插值多項式[2]、第二種邊界問題邊界條件

次邊界條件表明M0和Mn是已知的，所以(3.15)式直接改寫為

（3-17）3.1.6三次樣條插值多項式[2]、第二種邊界問題邊界條件a、

b、并注意到，得

記

3.1.6三次樣條插值多項式則方程可簡寫成

將上式和與(3.15)式合在一起

滿足第[1]或第[2]或第[3]種邊界條件的三次樣條插值函數(shù)S(x)是存在且唯一的。

（3-18）3.1.6三次樣條插值多項式例3-3：已知函數(shù)y=f(x)的函數(shù)值如下

x-1.5012y0.125-119求三次樣條函數(shù)S(x)，使其滿足邊界條件。解：這是第一種邊界條件問題，(3-14)式和(3-16)式的個參數(shù)為

3.1.6三次樣條插值多項式方程組為

求解后得

由(3-10)式得

3.2曲線擬合3.2.1問題描述實驗的目的是尋找一些物理量之間的關(guān)系，如一組實驗數(shù)據(jù)（xi,yi）,要尋找函數(shù)的一個近似表達式，這就是一個曲線擬合問題。插值方法可以在一定程度上解決曲線擬合問題，但有明顯的不足。首先，實驗數(shù)據(jù)有誤差，如果進行插值，則所求得的曲線必須嚴(yán)格通過所有實驗數(shù)據(jù)點，從而使得該曲線保留了實驗誤差；其次，實驗數(shù)據(jù)往往很多，用插值法求得的表達式因此缺乏實用性。

曲線擬合還有一個好壞的標(biāo)準(zhǔn)問題，標(biāo)準(zhǔn)不同決定著不同擬合的方法。常用的有最小二乘法，其基本思想是：使擬合曲線與實驗點之差的平方和最小。

3.2曲線擬合3.2.2概述實驗數(shù)據(jù)表

實驗編號

自變量因變量

1x1

y12x2

y23x3

y3

…

…

…ixiyi

…

…

…mxmym3.2曲線擬合

設(shè)x和y之間為線性關(guān)系，則有：

通過實驗數(shù)據(jù)來確定（3-19）式中的常數(shù)a和b，即建立x與y的關(guān)系，此為一個一元線性擬合問題（線性回歸）。

問題：如何用m組數(shù)據(jù)確定線性方程中a和b？見下圖。（3-19）3.2曲線擬合

如何確定哪一條直線是最好的？

a.

最小二乘法：使回歸的殘差平方和最小。

b.殘差：實驗數(shù)據(jù)yi與回歸方程計算的f(xi)之間的差，用qi

表示：

按最小二乘法的定義：

（3-20）（3-21）3.2曲線擬合2.`算法

由（3-21）式可知，Q是a和b的函數(shù)，即根據(jù)數(shù)學(xué)知識，要使Q最小，有：

,,

將（3-21）式代入前兩式：3.2曲線擬合

經(jīng)過計算得到：

（3-22）（3-23）3.2曲線擬合由（3-22）和（3-23）式可驗證：3.方差分析

如何衡量回歸直線與實驗數(shù)據(jù)之間的吻合程度？或回歸直線的可信度是多少？----這個問題由方差分析來完成。衡量標(biāo)準(zhǔn)：（幾個統(tǒng)計量）

a.殘差平方和

b.回歸平方和其中

3.2曲線擬合c.剩余標(biāo)準(zhǔn)差其中n=1d.相關(guān)系數(shù)e.綜合檢驗值

一般常用的是：相關(guān)系數(shù)R（R:0<=R<=1;R越接近1越好）和綜合檢驗值F（F為顯著性檢驗，它的值越大越好）。

3.2曲線擬合多元線性最小二乘法：自變量個數(shù)大于一，如：當(dāng)

Q/

bi=0時，得到一個線性方程組，用以求出bi

非線性最小二乘法：將非線性方程回歸轉(zhuǎn)化為線性形式。

如轉(zhuǎn)化為即相當(dāng)于若非線性方程不能回歸轉(zhuǎn)化為線性形式時，可采用后面求根法中的迭代法。3.2曲線擬合3.2.3程序設(shè)計程序功能是：由已知的m組實測數(shù)據(jù)（xi,yi),i=1,2,3,…,m,按最小二程序原理擬合多項式的系數(shù)。其中用戶可以根據(jù)具體問題預(yù)先確定擬合多項式的次數(shù)，也可以按精度要求，使程序自動選擇多項式的次數(shù)。應(yīng)該注意的是，無論何種選擇，其次數(shù)都不應(yīng)該超過m。1、變量及數(shù)組使用說明M:實測數(shù)據(jù)（xi,yi)的個數(shù)；N:擬合多項式次數(shù)加1，例初值為2，即為一次多項式；EPS:允許誤差；N0:用戶直接選擇多項式次數(shù)，=0時表示按精度要求自動選擇多項式次數(shù)；3.2曲線擬合X(M):存放給定點的值；Y(M):存放給定點的實測值；F(M):存放給定點所擬合的函數(shù)值；S(2M):存放運算中形成的正規(guī)方程組系數(shù)的所有值；A(M,M):存放正規(guī)方程組系數(shù)增廣矩陣；B(M):存放正規(guī)方程組常系數(shù)列向量；Z(M):調(diào)用解方程組子程序得到的解向量；C(M):存放多項式系數(shù)列向量；Q(M):存放殘量之絕對值。2、程序流程圖開始輸入N0,M,X(I),Y(I)N0=0?

N=2

N=N0+1計算S(K),K=1,2,…,2N-1生成S(I+J-1)=A(I,J)I=1,…N;J=1,…,N計算B(K)=A(K,N+1),K=1,…,N調(diào)用主元法解方程組子程序，解向量Z(I)X(I)=Z(I)I=1,…,N計算擬合曲線在給定點X(I)上函數(shù)值F(I)I=1,…,M計算殘量Q(I)=ABS(F(I)-Y(I))尋找最大值Q(I0)N0=0?Q(I0)/F(I0)<EPS?

N=N+1N<M輸出多項式系數(shù)C(I),I=1,…,M結(jié)束YNYNYNYN3.2曲線擬合3、源程序清單10

DIMENSIONX(M),Y(M),20&S(2*M),A(M,M),B(M),30&Z(M),P(M),Q(M)40C

輸入已知數(shù)據(jù)50READN0,M60IF(N0=0)THEN70READEPS80ENDIF90DOI=1,M100READ(*,*)X(I),Y(I)110ENDDO120IF(N0=0)THEN130N=2140ELSE150N=N0+1160ENDIF170C生成系數(shù)矩陣S(2N-1)180DOI=1,2*N-1190S(I)=0.0200DOJ=1,M210S(I)=S(I)+X(J)**(I-1)220ENDDO230ENDDO240DOI=1,N3.2曲線擬合250DOJ=1,N260A(I,J)=S(I+J-1)270ENDDO

280ENDDO290C

產(chǎn)生增廣矩陣300DOI=1,N310B(I)=0.0320DOJ=1,M330B(I)=B(I)+Y(J)*X(J)**(I-1)340ENDDO350A(I,N+1)=B