化工系統(tǒng)工程基礎(chǔ)lesson5序貫?zāi)K法子系統(tǒng)收斂課件

11.基本知識2.非線性代數(shù)方程數(shù)值解法3.非線性代數(shù)方程組數(shù)值解法4.微分方程組數(shù)值解法§2.4切割物流變量的收斂（常用數(shù)學(xué)方法）11.基本知識§2.4切割物流變量的收斂（常21.基本知識過程系統(tǒng)的數(shù)學(xué)模型是一組非線性方程組，切割流線變量的收斂問題實際上是一個迭代求解非線性方程組的問題：當(dāng)假設(shè)值x與計算值y之差小于收斂容差ε時：則x為切割流線變量的收斂解。21.基本知識過程系統(tǒng)的數(shù)學(xué)模型是一組非線性方程組，切334方程的兩種表達(dá)形式方程：隱式表達(dá)形式：顯式表達(dá)形式：1、對于不同表達(dá)形式下的同一方程，可能適合采用不同的求解方法。2、兩種表達(dá)形式可以互換，且每一種表達(dá)方式下都可能存在非同等效方程形式。例：4方程的兩種表達(dá)形式方程：隱式表達(dá)形式：5兩種表達(dá)方式轉(zhuǎn)換的通式：容易實現(xiàn)轉(zhuǎn)換難以匹配恰當(dāng)5兩種表達(dá)方式轉(zhuǎn)換的通式：容易實現(xiàn)轉(zhuǎn)換難以匹配恰當(dāng)6收斂單元：數(shù)值迭代求解非線性方程組的子程序。適用于收斂單元的數(shù)值方法一般應(yīng)盡可能滿足：（1）對初值的要求不高。切割變量的初值可根據(jù)流線的實際意義給出，要求初值組數(shù)少的方法更實用。（2）數(shù)值穩(wěn)定性好。好的迭代方法應(yīng)該對各種問題都能得到收斂的解。（3）收斂速度快。三個主要影響因素：①迭代次數(shù)；②函數(shù)G(x)的計算次數(shù)，即一次流程回路的模擬計算；③矩陣求逆的次數(shù)。（4）占用計算機存儲空間少。流程模擬計算量大，數(shù)值計算的存儲空間問題也需要考慮。這些要求很難同時滿足，可以根據(jù)實際問題有所側(cè)重。6收斂單元：數(shù)值迭代求解非線性方程組的子程序。7迭代過程迭代的結(jié)果：收斂－待解變量向方程的解逐漸、無限逼近。發(fā)散－待解變量離方程的解越來越遠(yuǎn)。振蕩－待解變量的取值不斷重復(fù)同樣的變化。收斂精度，收斂容差：憑經(jīng)驗，或工程精度要求人為設(shè)定，足夠小的正數(shù)收斂判據(jù)絕對誤差：相對誤差：7迭代過程收斂精度，收斂容差：憑8一、直接迭代法（Directsubstitutionmethod）（1）顯式方程（2）迭代公式：（3）收斂特性：n=1，線性收斂2.非線性代數(shù)方程數(shù)值解法收斂y=φ(x)xkxk+1Φ(xk)Φ(xk+1)y=x8一、直接迭代法（Directsubstitutionm9舉例：x3-27=0,x*=3發(fā)散振蕩收斂直接迭代收斂的必要條件：在解的附近，有|′(x)|<1，能收斂9舉例：x3-27=0,x*=3發(fā)散振蕩收10舉例：迭代公式：取：x0=1.0，εa=0.0110舉例：迭代公式：11直接迭代法－迭代收斂性的幾何解釋發(fā)散振蕩發(fā)散φ(x)xkxk+1φ(x)x1x2φ(x)振蕩不穩(wěn)定11直接迭代法－迭代收斂性的幾何解釋發(fā)散振蕩發(fā)散φ(x)xk12解代數(shù)方程組:步驟：（1）設(shè)初猜測值x（x＝3） （2）將x（＝3）代入方程[1]，解得yc（=2/3－2/3*x＝－1.33） （3）將yc賦值給y，即y＝y(tǒng)c，代入方程[2]， 解得xc（=4/3－2/3*（－1.33）＝2.22） （4）比較xc與猜測值x，判斷是否滿足精度要求：

a、如果，則轉(zhuǎn)向（5），

b、如果，則轉(zhuǎn)向（1），將xc作為新的猜測值重新計算，直至滿足a的要求 （5）計算結(jié)束，輸出結(jié)果12解代數(shù)方程組:步驟：（1）設(shè)初猜測值x（x＝3）133y＋2x＝22y＋3x＝4（3，-1．33）xy133y＋2x＝22y＋3x＝4（3，-1．33）xy14（1）顯式形式（2）迭代公式其中：松弛因子

(relaxationfactor)

當(dāng)，為直接迭代法（3）收斂特性

n=1，線性收斂φ(x)xkxk+1φ(xk)二、部分迭代法（partialsubstitutionmethod）對振蕩不收斂的情況很有效適當(dāng)改變，可能會改進(jìn)收斂性14（1）顯式形式φ(x)xkxk+1直接迭代法單調(diào)收斂時，可取ω>1，在原有正確方向基礎(chǔ)上加大步長。即外推。振蕩收斂時，可取0<ω<1，在原有正確方向基礎(chǔ)上縮小步長。即內(nèi)插。振蕩發(fā)散時，可取0<ω<1，在原有正確方向基礎(chǔ)上縮小步長。即內(nèi)插。單調(diào)發(fā)散時，可取ω<0，沿直接迭代的相反方向進(jìn)行改進(jìn)?？煞Q“反推”。15發(fā)散，ω＝－0.05振蕩，ω＝0.5收斂，ω＝1.4直接迭代法15發(fā)散，ω＝－0.05振蕩，ω＝0.5收斂，ω＝16直接迭代法16直接迭代法17（1）顯式形式（2）迭代公式連接(k－1)點、(k)點直線：與直線y＝x的交點：三、韋格斯坦法（Wegsteinmethod）xk-1xkxk+1(k-1)(k)(k+1)φ(x)yk-1=φ(xk-1)yk=φ(xk)17（1）顯式形式三、韋格斯坦法（Wegsteinmeth18公式：初始：設(shè)x（0），用直接迭代得到x（1）（兩點法）（3）收斂特性

1<n<2，超線性收斂當(dāng)，可能導(dǎo)致迭代失敗，“有界Wegstein法”：經(jīng)驗地將ω限制在某一范圍內(nèi)。與部分迭代法有區(qū)別18與部分迭代法有區(qū)別19（1）隱式形式（2）迭代公式基本思想：設(shè)法將非線性方程轉(zhuǎn)化為某種線性方程近似求解f(x)在xk點一階展開：

得到：xxkx*f(x)f(xk)θθxk+1四、牛頓法（Newtonmethod）方向步長19（1）隱式形式xxkx*f(x)f(xk)θθxk+1四20（3）收斂特性

n=2，二階收斂缺點：a.初值點選擇影響收斂結(jié)果，對初值選擇敏感；

b.需要用分析法或數(shù)值法獲得導(dǎo)數(shù)，不適合復(fù)雜函數(shù)（4）近似牛頓法用差分代替微分：用前兩點連線的斜率近似導(dǎo)數(shù)：（割線法，與Wegstein法公式相同）20（3）收斂特性21解析法與數(shù)值法要充分小到能夠獲得導(dǎo)數(shù)的一個精確的近似值，但又不能小到產(chǎn)生顯著的舍入誤差。21解析法與數(shù)值法要充分小到能夠獲得導(dǎo)數(shù)的一個精確的近似22單變量非線性方程數(shù)值解法比較22單變量非線性方程數(shù)值解法比較23在序貫?zāi)K法中，

迭代計算不可分割子系統(tǒng)以流股5為切斷物流假設(shè)X5的初值X50；調(diào)用單元模塊計算單元，求得單元的輸出變量X2，即X2=1(X1，X50)調(diào)用單元模塊計算單元，求得單元的輸出變量X3，即X3=2(X2，X6)調(diào)用單元模塊計算單元，變量X5的計算值X^5和不可再分塊的輸出X4，即 X^5＝35(X3)，X4＝34(X3)將X5的計算值直接作為假設(shè)值（直接迭代法），令X50=X^5

，重復(fù)(2)－(5)，進(jìn)行迭代計算，直至X5和X^5之間的差值滿足收斂精度ε。最后一次迭代計算時所得X2～X5的值即為該不可再分塊的解。23在序貫?zāi)K法中，

迭代計算不可分割子系統(tǒng)以流股5為切斷物24方程組的迭代求解迭代通式：

X(k+1)=X(k)

+X(k)收斂判據(jù)絕對誤差:相對誤差：方程組的表達(dá)方式3.非線性代數(shù)方程組數(shù)值解法24方程組的迭代求解方程組的表達(dá)方式3.非線性代數(shù)方程25一、直接迭代法與部分迭代法1、顯式表達(dá)：X=(X)2、迭代公式:X(k+1)

=X(k)

+

((X(k))-X(k))

即：xi(k+1)=xi(k)+(i(X(k))-xi(k))

(i=1,…,n)當(dāng)＝1時，轉(zhuǎn)化為直接迭代公式：X(k+1)=

(X(k))

即：xi(k+1)=i(X(k))(i=1,…,n)一個參數(shù)“一個方程迭代一個變量”25一、直接迭代法與部分迭代法1、顯式表達(dá)：X=(X)263、收斂性

線性收斂，n=1

對于隱式向顯式形式轉(zhuǎn)化時，方程與變量的匹配十分重要，即(X)

構(gòu)造對收斂有影響：收斂的充分條件：每一次迭代的Jacobian矩陣的|1|均小于1。收斂的必要條件：在解x=x*處，Jacobian矩陣的|1|小于1。根據(jù)初值選擇不同，收斂結(jié)果一般為局部收斂。263、收斂性線性收斂，n=127二、韋格斯坦法（Wegsteinmethod）1、顯式形式：X=(X)2、迭代公式：也可以表示為：n個參數(shù)27二、韋格斯坦法（Wegsteinmethod）1、顯式283、收斂性1<n<2，超線性收斂。“一個方程收斂一個變量”隱式形式F(X)=0,可以用割線法283、收斂性1<n<2，超線性收斂。29例：29例：30三、牛頓－拉夫森法（Newton-Raphsonmethod）1、隱式形式F(X)＝02、迭代公式函數(shù)向量F(X)

在X=X(k)

處一階臺勞展開，可得：

F(X)=F(X(k))+J(k)(X-X(k))=0迭代公式：X(k+1)=X(k)-(J(k))-1F(X(k))或加入阻尼因子t(k):X(k+1)=X(k)-t(k)(J(k))-1F(X(k))式中，Jacobian矩陣：30三、牛頓－拉夫森法（Newton-Raphsonmet313、收斂性n=2，二次收斂“每個變量的收斂與整個方程組有關(guān)”問題：J(k)如果奇異，將導(dǎo)致求逆失敗。原因1：待解方程本身冗余或有矛盾方程存在－－檢查和修正方程組原因2：迭代過程中出現(xiàn)奇異或近似奇異Jaconbian矩陣－－用下式替換：

X(k+1)=X(k)-(J(k)＋qI)-1F(X(k))

根據(jù)經(jīng)驗，q=1.1||J(k)||1

，可滿足要求矩陣的1－范數(shù)313、收斂性n=2，二次收斂32例：a=0.0001，X*=(1,8,4)T初值：X(0)=(2,10,5)T，直接迭代，20次，收斂

Wegstein，10次，收斂

Newton-Raphson，4次，收斂

X(0)=(4,10,5)T，直接迭代，22次，收斂

Wegstein，14次，收斂

X(0)=(4,3,5)T，發(fā)散

X(0)=(8,2,1)T，發(fā)散局部收斂性32例：a=0.0001，X*=(1,8,4)T局部收斂33建立殘差函數(shù)：0.0000[0.0000，0.0000，0.0000]T[1.0000，8.0000，4.0000]T50.0018[-0.0003，0.0017，0.0003]T[0.9995，8.0002，3.9999]T40.1729[0.0014，0.1711，-0.0247]T[0.9636，8.0185，3.9932]T35.5584[1.4871，5.3527，-0.1837]T[0.7084，8.3416，4.0334]T251.4196[18.414，48.000，0.951]T[2，10，5]T1n33建立殘差函數(shù)：0.0000[0.0000，0.0000，344、擬牛頓法(Quasi-Newtonmethod)的基本思路用矩陣H(k)代替Newton-Raphson法中的(J(k))-1公式：X(k+1)=X(k)-H(k)F(X(k))H(k)構(gòu)造方法不同－－一系列Quasi-Newton法344、擬牛頓法(Quasi-Newtonmethod)的35四、布洛伊頓法（Broydenmethod）（一種典型的擬牛頓法）1、迭代公式2、迭代初始點和初始迭代矩陣的確定初始點X(0)，初始迭代矩陣H(0)＝-(J(0))-1，或單位陣I3、收斂特性：1<n<2，超線性收斂“每個變量的收斂與整個方程組有關(guān)”35四、布洛伊頓法（Broydenmethod）（一種典型36作業(yè)：1、用直接迭代法求方程的解，并用|′(x)|<1進(jìn)行收斂性的驗證。2、分別用部分迭代法（其中ω=1，0.8，0.4），Wegstein法，牛頓法，求x值；并比較各種迭代法得到的結(jié)果。36作業(yè)：1、用直接迭代法求方程的解，并用|′(x)|<37作業(yè)：3、利用直接迭代法解方程：初值（0.5，0.5），a＝0.01，比較用兩種不同的顯式表達(dá)式計算得到的結(jié)果，并分別求其在收斂點處的Jacobian矩陣。37作業(yè)：3、利用直接迭代法解方程：38微分方程數(shù)值解法

在遇到的數(shù)學(xué)問題中，如果其中的某個變量和其它變量之間的函數(shù)依賴關(guān)系是未知的，但是這個未知函數(shù)及其某些階的導(dǎo)數(shù)連同自變量一起可由一個已知的方程描述，則這樣的方程稱為微分方程。常微分方程：只有一個自變量的微分方程。在化工中，描述穩(wěn)態(tài)一維分布參數(shù)問題和非穩(wěn)態(tài)集中參數(shù)問題的微分方程均屬于常微分方程。求解方法：解析法，數(shù)值法n階常微分方程38微分方程數(shù)值解法在遇到的數(shù)學(xué)問題中，如果39一、化工流程中常見的微分方程表達(dá)1、物質(zhì)A的分解則A的分解速度：反應(yīng)速率常數(shù)反應(yīng)物料A的濃度39一、化工流程中常見的微分方程表達(dá)1、物質(zhì)A的分解反應(yīng)速率402、套管內(nèi)冷熱兩流體溫度分布內(nèi)管內(nèi)和內(nèi)管外兩流體熱平衡方程：402、套管內(nèi)冷熱兩流體溫度分布內(nèi)管內(nèi)和內(nèi)管外兩流體熱平衡方413、連續(xù)流動槽內(nèi)液面隨時間變化規(guī)律已知：液體流入量Qi(t)，流出量Qo(t)，容器內(nèi)橫截面積A

求：容器內(nèi)流體液面高度H隨時間t變化的規(guī)律z(t)＝？動態(tài)模型：容器內(nèi)積累速率＝液相流進(jìn)速率－液相流出速率413、連續(xù)流動槽內(nèi)液面隨時間變化規(guī)律已知：液體流入量Qi(424、常見的偏微分方程化工中偏微分方程很多，具有代表性的包括非穩(wěn)態(tài)傳遞過程、非穩(wěn)態(tài)反應(yīng)器模型等?？紤]軸徑向返混的非穩(wěn)態(tài)固定床反應(yīng)器模型，其主要物料衡算方程即是典型的二維偏微分方程：二維傳質(zhì)、傳熱過程的微分方程，其通用形式為：424、常見的偏微分方程化工中偏微分方程很多，具有代表性的包43

乙烯工業(yè)裂解過程的數(shù)學(xué)模型質(zhì)量平衡方程

動量平衡方程

能量平衡方程

43乙烯工業(yè)裂解過程的數(shù)學(xué)模型2023/7/2744乙烯工業(yè)裂解過程的數(shù)學(xué)模型過程模型C2H6→C2H4+H2C3H6→C2H2+CH4……反應(yīng)模型質(zhì)量衡算方程動量衡算方程熱量衡算方程C3H8→C3H6+H22023/7/2744乙烯工業(yè)裂解過程的數(shù)學(xué)模型過程模45乙烯裂解自由基反應(yīng)模型45乙烯裂解自由基反應(yīng)模型46二、常微分方程數(shù)值解法簡介初值問題：時，已知。求解進(jìn)程是“步進(jìn)式”。一般情況下，初始條件需給出所研究問題的整個系統(tǒng)的狀態(tài)，而非某一特定點的數(shù)值。只有對于均勻場問題，或集中參數(shù)問題，在給定初始條件時才可以以點代面。邊值問題：需要迭代求解。數(shù)值解法：尋求y(x)在一系列離散節(jié)點x1<x2<…<xn<xn+1<…上的近似值，y1，y2，…，yn，yn+1…

相鄰兩個節(jié)點的間距h=xn+1-xn稱為步長，h為定數(shù)。

xn=x0+nh，n=0，1，2，…解析法：只能用來求解一些特殊類型的方程。46二、常微分方程數(shù)值解法簡介初值問題：47三、一階常微分方程初值問題（Initialvalueproblem

）1、歐拉法（一階）沿點P0(x0，y0)切線方向f(x0，y0)推進(jìn)到x=x1，得到P1點。直線：y=y0+f(x0，y0)(x-x0)其中，x1=x0+h0

y1=y0+h0f(x0，y0)從P1點，沿y=y(x)的切線方向f(x1，y1)推進(jìn)，得到：y2=y1+h1f(x1，y1)

……

……yn+1=yn+hnf(xn，yn)取固定步長h，得：yn+1=yn+hf(xn，yn)xn=x0+nhy=y(x)y(x)=y0y(x1)y1x0x1P0P1Euler公式47三、一階常微分方程初值問題（Initialvalue48局部截斷誤差一般，假設(shè)yn是準(zhǔn)確的，即yn=y(xn)

把y(xn+1)-yn+1稱為局部截斷誤差。Euler公式：yn+1=yn+hf(xn，yn)

=yn+hy’(xn)臺勞展開式：y(xn+1)

=y(xn)

+hy’(xn)+(1/2)h2y’’()y(xn+1)

-yn＋1=

(1/2)h2y’’()(1/2)h2y’’(xn)48局部截斷誤差一般，假設(shè)yn是準(zhǔn)確的，即yn=y(49例：求混合槽內(nèi)濃度變化的數(shù)學(xué)模型t=0時，槽內(nèi)濃度y0=0.1進(jìn)口濃度yin=0.9槽內(nèi)液體容量不變：V＝50mol進(jìn)口、出口液體流量F＝L＝100mol/h求：t＝0.2時，y值。數(shù)學(xué)模型：49例：求混合槽內(nèi)濃度變化的數(shù)學(xué)模型t=0時，槽內(nèi)濃度y0=502、改進(jìn)的歐拉法（二階法）一階法：用每一小段初始導(dǎo)數(shù)代替整個小段上的導(dǎo)數(shù)，步長，誤差二階法：用每一小段初始和末端導(dǎo)數(shù)的平均值代替整個小段上的導(dǎo)數(shù)。

yn+1

=yn+(1/2)h[f(xn，yn)+f(xn+1，yn+1)]

初值：y0n+1=yn+hf(xn，yn)

迭代滿足：|yk+1n+1-

ykn+1|

如果只進(jìn)行一步迭代，形成二階Euler公式：

yn+1=yn+(1/2)K1+(1/2)K2K1=hf(xn，yn)K2=hf(xn+h，yn+K1)502、改進(jìn)的歐拉法（二階法）一階法：用每一小段初始導(dǎo)數(shù)代替513、四階龍格－庫塔法（Runge-Kuttamethod）基本思路：用差商代替微分，利用微分中值定理，

存在0<

<1，使得[y(xn+1)-y(xn)]/h=y’(xn+h)

得到y(tǒng)(xn+1)=y(xn)+hf[xn+h，y(xn+h)]K*，區(qū)間[xn，xn+1]上的平均斜率

K*的一種計算方法－－導(dǎo)出y(xn+1)的一種計算公式四階R-K法yn+1=yn+(1/6)h[K1+2K2+2K3+K4]K1=f(xn，yn)K2=f(xn+h/2，yn+hK1/2)K3=f(xn+h/2，yn+hK3/2)K4=f(xn+h，yn+hK3)513、四階龍格－庫塔法（Runge-Kuttamet524、一階常微分方程組初值問題數(shù)值解法一階常微分方程組：

y’1=f1(x，y1，y2，…，ym)初始條件：y1(x0)=y01

y’2=f2(x，y1，y2，…，ym)y2(x0)=y02

……

……

……y’m=fm(x，y1，y2，…，ym)ym(x0)=y0m用向量表示，可記為：Y’=F(X，Y)Y(X0)=Y0四階R-K法yi，n+1=yi，n+(1/6)h[Ki，1+2Ki，2+2Ki，3+Ki，4]Ki，1=fi(xn，y1，n，y2，n，…，ym，n)Ki，2=fi(xn+h/2，y1，n+hK1，1/2

，y2，n+hK2，1/2

…ym，n+hKm，1/2)Ki，3=fi(xn+h/2，y1，n+hK1，2/2

，y2，n+hK2，2/2

…ym，n+hKm，2/2)Ki，4=fi(xn+h，y1，n+hK1，3，y2，n+hK2，3…ym，n+hKm，3)524、一階常微分方程組初值問題數(shù)值解法一階常微分方程組：四53常用數(shù)值積分解法及比較方法名稱優(yōu)點缺點適用范圍歐拉法程序簡單、快速精度低、穩(wěn)定性差高階、低精度、實時仿真二階龍格－庫塔法程序簡單、快速精度低、穩(wěn)定性差高階、低精度、實時仿真定步長四階龍格－庫塔法程序簡單、精度較高、小步長、可解剛性方程高階較費時非剛性、不太光滑中等精度的普通方程變步長四階龍格－庫塔法可以保證滿足指定精度高階較費時同上，仿真的實時性好一些拉姆伯特五階方法精度高計算量大、高階費時，不能解剛性方程高精度計算，檢驗用默森單步法精度高計算量大，高階費時，不能解剛性方程高精度計算、檢驗用53常用數(shù)值積分解法及比較方法名稱優(yōu)點缺點適用范圍歐拉法程序54常用數(shù)值積分解法及比較方法名稱優(yōu)點缺點適用范圍外插法精度高穩(wěn)定性不高高階光滑非剛性方程阿當(dāng)姆斯預(yù)報－校正法可以估計誤差，穩(wěn)定性較好不能解剛性方程高階光滑非剛性方程小參數(shù)法簡單，變步長，速度快，可解剛性方程精度中等高階剛性方程實時仿真特雷納法可解剛性方程只適用于特殊情況特殊剛性方程吉爾法變階，變步長，適用于非線性剛性方程占用容量大，高階問題困難非線性剛性方程改進(jìn)吉爾法變階變步長，適于非線性剛性及混合方程程序復(fù)雜，須解非線性代數(shù)方程組高階非線性剛性混合方程54常用數(shù)值積分解法及比較方法名稱優(yōu)點缺點適用范圍外插法精度AspenPlus的流程收斂55切割物流：默認(rèn)方法：Wegstein撕裂物流默認(rèn)處置為零賦予初值優(yōu)先選擇設(shè)計規(guī)定：反饋控制器循環(huán)迭代求解默認(rèn)方法：Secant目的方法收斂切割物流收斂設(shè)計規(guī)定優(yōu)化Wegstein,direct,Broyden,NewtonSecant,Broyden,NewtonSQP,Compl