線性擬合方法

1、第五章第五章 實驗數(shù)據(jù)及模型參數(shù)擬合方法實驗數(shù)據(jù)及模型參數(shù)擬合方法第五章第五章 實驗數(shù)據(jù)及模型參數(shù)擬合方法實驗數(shù)據(jù)及模型參數(shù)擬合方法第一節(jié)第一節(jié) 問題的提出問題的提出第一節(jié)第一節(jié) 問題的提出問題的提出 在化工設計及化工模擬計算中，需要大量的物性參數(shù)及各種設備參數(shù)。這些參數(shù)有些可以通過計算得到，但大量的參數(shù)還是要通過實驗測量得到。實驗測量得到的常常是一組離散數(shù)據(jù)序列（xi ,yi）。 如果數(shù)據(jù)序列（xi ,yi）(為一般起見), i=1,2, ,m ,含有不可避免的誤差（或稱“噪聲” ，如圖5-1所示），或者無法同時滿足某特定的函數(shù)（如圖5-2所示）,那么，只能要求所作逼近函數(shù)（x）最優(yōu)地靠近樣

2、點，即向量Q=（(x1), (x2), , (xm)）T與Y=(y1,y2, ,ym)T的誤差或距離最小。按Q與Y之間誤差最小原則作為“最優(yōu)”標準構造的逼近函數(shù)，稱為擬合函數(shù)。第一節(jié)第一節(jié) 問題的提出問題的提出-202468101214161820050100150200Y X 024681005101520Y X 圖5-1 含有噪聲的數(shù)據(jù)圖5-2 無法同時滿足某特定函數(shù)的數(shù)據(jù)序列第一節(jié)第一節(jié) 問題的提出問題的提出 除了物性數(shù)據(jù)及設備參數(shù)需要利用數(shù)據(jù)擬合外，在化學化工中，許多模型也要利用數(shù)據(jù)擬合技術，求出最佳的模型和模型參數(shù)。如在某一反應工程實驗中，我們測得了如表5-1所示的實驗數(shù)據(jù)。 序號

3、1 2 3 4 5 6 7 8 溫度 T 10 20 30 40 50 60 70 80 轉化率 y 0.1 0.3 0.7 0.94 0.95 0.68 0.34 0.13 現(xiàn)在要確定在其他條件不變的情況下，轉化率y和溫度T的具體關系，現(xiàn)擬用兩種模型去擬合實驗數(shù)據(jù)，兩種模型分別是 2111TcTbay2222)45(Tbacy第一節(jié)第一節(jié) 問題的提出問題的提出 如何求取上述模型中的參數(shù)，并判斷兩種模型的優(yōu)劣,是化學化工工作者經(jīng)常要碰到的問題，這個問題的求解將在本章下面的有關章節(jié)中進行詳細的講解。 2111TcTbay2222)45(Tbacy第二節(jié)第二節(jié) 擬合的標準擬合的標準第二節(jié)第二節(jié) 擬

4、合的標準擬合的標準 前面已經(jīng)提到按Q與Y之間誤差最小原則作為“最優(yōu)”標準構造的逼近函數(shù)，稱為擬合函數(shù)，而向量Q與Y之間的誤差或距離有各種不同的定義方法，一般有以下幾種。 （1）用各點誤差絕對值的和表示 （2）用各點誤差按絕對值的最大值表示 （3）用各點誤差的平方和表示 miiiyxR11)(iimiyxR)(max122212Y-Q(x)R )(或imiiyxRR 式中R稱為均方誤差。由于計算均方誤差的最小值的原則容易實現(xiàn)而被廣泛采用。按均方誤差達到極小構造擬合曲按均方誤差達到極小構造擬合曲線的方法稱為最小二乘法線的方法稱為最小二乘法。同時還有許多種其他的方法構造擬合曲線，感興趣的讀者可參閱有

5、關教材。本章主要講述用最小二乘法構造擬合曲線。第二節(jié)第二節(jié) 擬合的標準擬合的標準 第二節(jié)第二節(jié) 擬合的標準擬合的標準_實例實例1 1 實驗測得二甲醇（DME）的飽和蒸氣壓和溫度的關系，見表5-2。 序號溫度 蒸氣壓 MPa1-23.70.1012-100.174300.2544100.3595200.4956300.6627400.880表5-2 DME飽和蒸氣壓和溫度的關系-30-20-10010203040500.00.20.40.60.81.0pt 由表5-2的數(shù)據(jù)觀測可得，DME的飽和蒸氣壓和溫度有正相關關系，如果以函數(shù)p=a+bt來擬合,則擬合函數(shù)是一條直線。通過計算均方誤差Q (

6、a , b )最小值而確定直線方程。 （見圖5-3 ）圖5-3 DME飽和蒸汽壓和溫度之間的線性擬合 第二節(jié)第二節(jié) 擬合的標準擬合的標準_實例實例1 1 2121)()(),(imiiimiipbtaptpbaQ 擬合得到直線方程為： 相關系數(shù)R為0.97296，平均絕對偏差SD為0.0707。 tp0.01210.30324擬合的標準擬合的標準 實例實例2 2 如果采用二次擬合，通過計算下述均方誤差 擬合得二次方程為 相關系數(shù)R為0.99972，平均絕對偏差SD為0.00815,具體擬合曲線見圖5-4。21221012210)()(),(imiiimiiiptataaptpaaaQ20001

7、50009570248450t.t.p-30-20-10010203040500.00.20.40.60.81.0y=0.24845+0.00957 x+0.00015 x2壓力, P(MPa)溫 度, t()圖5-4 DME飽和蒸氣壓和溫度之間的二次擬合 擬合的標準擬合的標準 實例實例2 2 比較圖5-3和圖5-4以及各自的相關系數(shù)和平均絕對偏差可知：n對于DME飽和蒸汽壓和溫度之間的關系，在實驗溫度范圍內用二次擬合曲線優(yōu)于線性擬合。n二次擬合曲線具有局限性，由圖5-4觀察可知，當溫度低于-30時，飽和壓力有升高的趨勢，但在擬合的溫度范圍內，二次擬合的平均絕對偏差又小于一次擬合，故對物性數(shù)據(jù)

8、進行擬合時，不僅要看在擬合條件下的擬合效果，還必須根據(jù)物性的具體性質，判斷在擬合條件之外的物性變化趨勢，以便使擬合公式在已做實驗點數(shù)據(jù)之外應用。 第三節(jié)第三節(jié) 單變量擬合和多變量擬合單變量擬合和多變量擬合 給定一組數(shù)據(jù)（xi,yi）,i=1, 2 , , m ， 作擬合直線p (x)=a + bx , 均方誤差為 2121)()(),(imiiimiiybxayxpbaQ由數(shù)學知識可知，Q (a , b)的極小值需滿足： 0)(2),(1imiiybxaabaQ0)(2),(1iimiixybxabbaQ整理得到擬合曲線滿足的方程： mimimiiiiimimiiiyxbxaxybxma111

9、211)()()(5.3.1 單變量擬合單變量擬合 1. 單變量擬合單變量擬合 線性擬合線性擬合1.單變量擬合單變量擬合線性擬合線性擬合 該方程可用消元法或克萊姆方法解出方程（如下式所示） )()()(/()( 2112111211211121121112111mimiiimimiimiiiimiimiimiiimiimiimiimiimiimiimiimiiimiimiixxmyxyxmbxxmyxxxyxxxmxyxxya單變量擬合單變量擬合 線性擬合實例線性擬合實例例：例： 下表為實驗測得的某一物性和溫度之間的關系數(shù)據(jù)，表中x為溫度數(shù)據(jù)，y為物性數(shù)據(jù)。請用線性函數(shù)擬合溫度和物性之間的關系

10、。 解：設擬合直線p(x)=a+bx ,并計算得下表： x7911131517192123252729y91215182124273033363942x313335373941434547 y454851545760636669 編號xyxyx21234521791113154756791215182169819631081652343153243267334981121169225220918389 將數(shù)據(jù)代入法方程組（1-12）中，得到： 解方程得：a = -1.5 , b = 1.5 擬合直線為： 267338191838956756721bax .p( x ) 5 .151相關系數(shù)R為

11、1。單變量擬合單變量擬合 線性擬合實例線性擬合實例線性擬合線性擬合VBVB清單清單Private Sub Command1_Click()Dim x(5),y(5),c,d,m,p,a,b,eerConst n=5For i =1 To 5 x(i)=InputBox(“x(i)=”) y(i)=InputBox(“y(i)=”)Print “x(i)=”;x(i)Print “y(i)=”;y(i)Next I c=0 d=0 m=0 p=0 For i=1 To 5 c=c+x(i) d=d+x(i)2 m=m+y(i) p=p+x(i)*y(i)Next i a=(m*d-c*p)/(n

12、*d-c2) b=(n*p-c*m)/(n*d-c2) 參數(shù)計算 a=Int(a*1000+0.5)/1000 b=Int(b*1000+0.5)/1000Text1.Text=Str(a)Text2.Text=Str(b) 參數(shù)輸出For i=1 To 5 eer=eer+(a+b*x(i)-y(i)2 誤差計算 eer=Int(eer*100000+0.5)/100000Next i eer=eer/5Text3.Text=Str(eer)End Sub有關有關線性線性擬合變型問題擬合變型問題 例如要擬合y=a+b/x2，只需在數(shù)據(jù)輸入后增加一語句x(i)=1/x(i)2，而在程序后面的誤

13、差eer 的計算中則不需要修改。2.單變量擬合單變量擬合 二次擬合函數(shù)二次擬合函數(shù) 給定數(shù)據(jù)（xi ,yi), i=1, 2 , , m ,用二次多項式函數(shù)擬合這組數(shù)據(jù)。 設 ,作出擬合函數(shù)與數(shù)據(jù)序列的均方誤差表達式2210 x ax a ap ( x ) 21221012210)()(),(imiiimiiiyxaxaayxpaaaQ由數(shù)學知識可知，Q( a0 ,a1 ,a2 )的極小值滿足 ：miiiiimiiiiimiiiixyxaxaaaQxyxaxaaaQyxaxaaaQ12221021221011221000)(20)(20)(2整理左式得二次多項式函數(shù)擬合的滿足條件方程（5-14

14、）： miiimiiimiimiimiimiimiimiimiimiimiiyxyxyaaaxxxxxxxxm121121014131213121121(5-14 ) 解此方程得到在均方誤差最小意義下的擬合函數(shù)p ( x )。式（5-14）稱為多項式擬合的法方程，法方程的系數(shù)矩陣是對稱的。當擬合多項式n 5時，法方程的系數(shù)矩陣是病態(tài)的，在用通常的迭代方法求解線性方程時會發(fā)散，在計算中要采用一些特殊算法以保護解的準確性。關于線性方程的求解方法，已在第三章中介紹。2.單變量擬合單變量擬合 二次擬合函數(shù)二次擬合函數(shù)3.二次擬合函數(shù)的拓展二次擬合函數(shù)的拓展 和一次擬合一樣，二次擬合也可以有多種變型，例

15、如 套用上面的公式，可以得到關于求解此擬合函數(shù)的法方程（5-15）。值得注意的是在此法方程的構建過程中，進行了變量的代換。首先是擬合函數(shù)中變量的代換: 。 P(x)=a0+a1x3+a2x5 253,xxxxmiiimiiimiimiimiimiimiimiimiimiimiiyxyxyaaaxxxxxxxxm1513121011018151815131513（5-15）其次是法方程的代換：將相應擬合函數(shù)中的代換引入法方程中。同時應注意法方程中x的4次冪是由兩個2次冪相乘得到，x的3次冪是由一個2次冪和一個1次冪相乘得到，而2次冪就是變量本身，而非兩個1次冪相乘得到。這個概念至關重要，在以后的

16、二次擬合的各類變型中，均需利用這個概念，千萬不要用常規(guī)的思路去進行代入計算。 3.二次擬合函數(shù)的拓展二次擬合函數(shù)的拓展如果我們需要求解是下面的擬合函數(shù)： 5 . 1110)273(273lnxbxaaymiiimiiimiimiimiimiimiimiimiimiimiiyxxyyaaaxxxxxxxxm15 . 1112101315 . 015 . 115 . 015 . 1115 . 11ln)273(273lnln)273()273()273()273()273(2731)273(2731參照上面的方法，我們很容易得到求解該擬合函數(shù)的法方程3.二次擬合函數(shù)的拓展二次擬合函數(shù)的拓展4.二次

17、擬合實例二次擬合實例請用二次多項式函數(shù)擬合下面這組數(shù)據(jù)。解解：設 并計算得下表序 號 1 2 3 4 5 6 7 x -3 -2 -1 0 1 2 3 y 4 2 3 0 -1 -2 -5 2210 x ax a ap ( x ) 序號xyxyX2X2yX3X41-34-12936-27812-22-448-8163-13-313-114000000051-1-11-11162-2-44-881673-5-159-45278101-3928-701964.4.二次擬合實例二次擬合實例將上面數(shù)據(jù)代入式 (5-14) ,相應的法方程為7196028a390280a12807a210210210aa

18、aaaa解方程得： a0 =0.66667 , a1 = -1.39286 , a2 = -0.13095所以：2130950392861666670 x.x- . - .p (x ) -3-2-10123-6-4-2024y=0.66667-1.39286 x-0.13095 x2yY X 圖圖 5-6 擬合曲線與數(shù)據(jù)序列擬合曲線與數(shù)據(jù)序列 二次擬合二次擬合-VB-VB程序清單程序清單Private Sub Command1_Click()Dim n, m As Integern = InputBox(n=方程次數(shù))m = InputBox(m=實驗次數(shù))ReDim x(m), y(m),

19、a0(n + 1), a1(n + 1), aa(m, n + 1)ReDim qq(n + 1, m), pp(n + 1, n + 1), b(n + 1), g(n + 1), tt(n + 1, n + 1)For I = 1 To m: 讀入數(shù)據(jù)READ x(I), y(I)Next IData -3, 4, -2, 2, -1, 3, 0, 0, 1, -1, 2, -2, 3, -5For i = 1 To mx(i) = InputBox(x(i)y(i) = InputBox(y(i)Next i二次擬合二次擬合-VB-VB程序清單程序清單omiga = InputBox(o

20、miga=松弛因子)For j = 1 To n + 1為計算法方程中的系數(shù)做準備 For i = 1 To m aa(i, j) = x(i) (j - 1) Next iNext jFor j = 1 To m For i = 1 To n + 1 qq(i, j) = aa(j, i) Next iNext j計算法方程中的右邊項For i = 1 To n + 1b(i) = 0For j = 1 To m b(i) = b(i) + aa(j, i) * y(j)Next jNext i開始計算法方程中的右邊項的系數(shù)For i = 1 To n + 1For j = 1 To n +

21、 1 pp(i, j) = 0 For k = 1 To m pp(i, j) = qq(i, k) * aa(k, j) + pp(i, j) Next kNext jNext i二次擬合二次擬合-VB-VB程序清單程序清單開始用松弛迭代法求解法方程中的變量For i = 1 To n + 1 a0(i) = 0 a1(i) = 0.11Next iFor i = 1 To n + 1 g(i) = b(i) / pp(i, i) For j = 1 To n + 1If j = i Then tt(i, j) = 0Else tt(i, j) = -pp(i, j) / pp(i, i)E

22、nd IfNext jNext i50 eer = 0 For i = 1 To n + 1eer = eer + Abs(a0(i) - a1(i) Next iIf eer 2 )。一般地，將含有n個未知量m個方程的線性方程組其一般形式為第四節(jié)解矛盾方程組第四節(jié)解矛盾方程組mnmnmmnnnnyxaxaxayxaxaxayxaxaxa22112222212111212111 mnmnmmnnyyyxxxaaaaaaaaa2121212222111211 寫成矩陣形為 一般情況下，當方程數(shù)m多于變量數(shù)n ，且m個方程之間線性無關, 則方程組無解，這時方程組稱為矛盾方程組。方程組在一般意義下無

23、解，也即無法找到n個變量同時滿足m個方程。這種情況和擬合曲線無法同時滿足所有的實驗數(shù)據(jù)點相仿，故可以通過求解均方誤差 極小意義下矛盾方程的解來獲取擬合曲線。 由數(shù)學知識還可證明：方程組ATAX = ATY的解就是矛盾方程組AX = Y 在最小二乘法意義下的解。這樣我們只要通過求解ATAX = ATY就可以得到矛盾方程組的解，進而得到各種擬合曲線，為擬合曲線的求解提供了另一種方法。第四節(jié)解矛盾方程組第四節(jié)解矛盾方程組minAX-B22例如，擬合直線p (x ) = a0 +a1x的矛盾方程組 ATAX = ATY的形式如下： 化簡得到與式(1-12)相同的法方程第四節(jié)解矛盾方程組第四節(jié)解矛盾方程

24、組mmmmyyyxxxaaxxxxxx2121102121 111111 111miiimiimiimiimiiyxyaaxxxm11101211 這里需要注意的是變量X和系數(shù)（a0 ， a1）之間的相互轉換關系。即 對于n次多項式曲線擬合，要計算Q ( a0 ,a1 , , an ) 的極小值問題,這與解矛盾方程組第四節(jié)解矛盾方程組第四節(jié)解矛盾方程組mxxxAaaX111 ,2110miinini - yx+ a+ x + a a1210) (mnmnmnnnnyxaxaayxaxaayxaxaa102221011110 mnyyyaaaA2110 與求 的極小問題是一回事。minaaa12

25、ini110y-xx或在這里 nmmnnxxxxxxA111 2211故對離散數(shù)據(jù)（xi,yi）,i=1, 2 , , m ；所作的n次擬合曲線y= xa + +xa +anini10，可通過解下列方程組求得： mTnTyyy A aaaAA2110 （1-21） 第四節(jié)解矛盾方程組第四節(jié)解矛盾方程組n如果擬合函數(shù)有n個自變量并進行一次擬合，則其擬合函數(shù)為： （1-22） miinimiinimiiimiiimiinnminiminiminiminiminiminiminiminiminimiimiimiiminimiimiiyxyxyxyxyaaaaaxxxxxxxxxxxxxxxm1111

26、211121012121111121111111131211121xnnnkkxaxaxaxaxaay1122110通過m（mn）次實驗，測量得到了m組 ),(, 1, 2, 1ininikiiixxxxxy的實數(shù)據(jù)，則可得到上面n個自變量擬合函數(shù)的法方程第四節(jié)解矛盾方程組第四節(jié)解矛盾方程組n只要對法方程（1-22）稍加修改，就可以得到有n個自變量的任意次方的擬合函數(shù)的法方程，通過法方程的求，就可以得到擬合函數(shù)中的各項系數(shù)。 miiinmiiinmiiimiiimiinnmininmiiinmiiinmiinmiininmiinmiiinmiinmiinimiiimiimiimiinmiimiiyxyxyxyxyaaaaxxxxxxxxxxxxxxxxxxxxxm1,1, 112,1, 1111012,1, 2,1, 1,1,1, 112, 11, 1, 11, 11, 11, 2,.112, 11, 11,1, 21, 1（1-23） 第四節(jié)解矛盾方程組第四節(jié)解矛盾方程組利用解矛