插值法与最小二乘法课件

插值法與最小二乘法

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

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

的代數多項式

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

在每個插值點構造插值基函數為一個

次多項式，且滿足條件

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

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

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

的近似值。解：令用為節點，有一次和二次拉格朗日插值公式求用為節點，有拉格朗日插值多項式程式設計計算公式(3-5)的程式為二重迴圈。由內迴圈（j迴圈）通過累乘求得基函數，然後通過外迴圈（i迴圈）得到插值結果y.下圖描述了拉格朗日演算法流程，對於該框圖中的有關參數說明如下：N:

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

插值點上的函數值LN:拉格朗日插值結果T:

存放累乘積I:

外迴圈變數J:

內迴圈變數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、插值餘項對插值後，用代替必定會產生誤差，其誤差可表示為

稱上式為插值多項式

的插值餘項。

設在區間上有直到階導數，為

上個互異的節點，為滿足條件的

次插值多項式。

那麼對於任何，有

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

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

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

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

取節點為等距離，則

其中h為步長。在兩相臨節點處的函數值之差為

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

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

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

在節點上，函數值是已知的，所以

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

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

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

3.1.4牛頓插值多項式

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

依次類推

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

給出列表函數f(x)=lgx的值，如下表所示，試用牛頓插值法求lg4.01的值。表3.3解：根據給定的函數作差商表，如表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)一階差商二階差商三階差商將以上各插值點值、函數值和差商值代入插值公式得lg4.01=0.60314433.1.5

Hermite插值多項式

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

Hermite插值多項式例3-2已知函數

的函數值、導數值如下，求其插值多項式及誤差運算式。-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、三次樣條插值函數的定義對於一組已知的數據，且，若函數滿足：

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

⑶、插值條件則稱為函數f(x)關於節點的三次樣條插值函數。

3.1.6三次樣條插值多項式

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

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

3.1.6三次樣條插值多項式

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

[1]、在邊界上給定一階導數的值；

[2]、給定二階導數的值；（如果，稱其為自然邊界條件

[3]、f(x)是以b-a為週期的函數，則要求S(x)及其導數也是以b-a為週期的函數，相應邊界條件為

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

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

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

所以

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

共n+1個），就可以完全確定S(x)。在確定(3.9)式中待定係數時用了插值點的條件，下麵利用節點一階導數連續條件（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個未知數。因此，要確定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]種邊界條件的三次樣條插值函數S(x)是存在且唯一的。

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

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

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

求解後得

由(3-10)式得

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

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

3.2曲線擬合3.2.2概述實驗數據表

實驗編號

引數因變數

1x1

y12x2

y23x3

y3

…

…

…ixiyi

…

…

…mxmym3.2曲線擬合

設x和y之間為線性關係，則有：

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

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

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

a.

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

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

表示：

按最小二乘法的定義：

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

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

,,

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

經過計算得到：

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

如何衡量回歸直線與實驗數據之間的吻合程度？或回歸直線的可信度是多少？----這個問題由方差分析來完成。衡量標準：（幾個統計量）

a.殘差平方和

b.回歸平方和其中

3.2曲線擬合c.剩餘標準差其中n=1d.相關係數e.綜合檢驗值

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

3.2曲線擬合多元線性最小二乘法：引數個數大於一，如：當

Q/

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

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

如轉化為即相當於若非線性方程不能回歸轉化為線性形式時，可採用後面求根法中的迭代法。3.2曲線擬合3.2.3程式設計程式功能是：由已知的m組實測數據（xi,yi),i=1,2,3,…,m,按最小二程式原理擬合多項式的係數。其中用戶可以根據具體問題預先確定擬合多項式的次數，也可以按精度要求，使程式自動選擇多項式的次數。應該注意的是，無論何種選擇，其次數都不應該超過m。1、變數及數組使用說明M:實測數據（xi,yi)的個數；N:擬合多項式次數加1，例初值為2，即為一次多項式；EPS:允許誤差；N0:用戶直接選擇多項式次數，=0時表示按精度要求自動選擇多項式次數；3.2曲線擬合X(M):存放給定點的值；Y(M):存放給定點的實測值；F(M):存放給定點所擬合的函數值；S(2M):存放運算中形成的正規方程組係數的所有值；A(M,M):存放正規方程組係數增廣矩陣；B(M):存放正規方程組常係數列向量；Z(M):調用解方程組副程式得到的解向量；C(M):存放多項式係數列向量；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調用主元法解方程組副程式，解向量Z(I)X(I)=Z(I)I=1,…,N計算擬合曲線在給定點X(I)上函數值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輸出多項式係數C(I),I=1,…,M結束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

輸入已知數據50READN0,M60IF(N0=0)THEN70READEPS80ENDIF90DOI=1,M100READ(*,*)X(I),Y(I)110ENDDO120IF(N0=0)THEN130N=2140ELSE150N=N0+1160ENDIF170C生成係數矩陣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

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