cfd2020第10講常微分方程數(shù)值解法

計算流體力學(xué)講義2020

第十講常微分方程數(shù)值解法

李新亮申義慶；力學(xué)所主樓214；82543801

；力學(xué)所主樓308；82543877

知識點：Runge-Kutta方法

線性多步法

剛性常微分方程的數(shù)值解法

1CopyrightbyLiXinliang課件下載：4階3階（TVD型）2階CopyrightbyLiXinliang2NS方程的（空間）離散推薦方法：Runge-Kutta法更高階……知識回顧：時間推進空間離散FDMFVMDGCopyrightbyLiXinliang3系統(tǒng)-----積分-微分方程組-----微分方程組------常微分方程簡單常微分方程的分析解常系數(shù)線性常微分方程、可分離變量型方程、全微分方程等常微分方程的數(shù)值計算計算機的發(fā)展、絕大多數(shù)常微分方程難以求解析解

分離變量線性方程(常數(shù)變易）全微分方程CopyrightbyLiXinliang4初值條件常微分方程（組）：m階常微分方程:引入輔助函數(shù)：……初值條件CopyrightbyLiXinliang5§10.1

單步法（Runge-Kutta方法）初值條件假定常微分方程的初值問題是唯一可解的因為f(x,y(x))恰好是待求精確解y(x)的斜率y’(x)y(x)是在x-y平面上經(jīng)過點(x0,y0)的曲線,該曲線在點(x0,y0)處的切線的斜率為f(x0,y0).Euler方法一階向前差商代替一階導(dǎo)數(shù)兩邊從

到積分,積分采用左矩形公式利用Taylor公式

CopyrightbyLiXinliang6單步法計算的值只用到前一步的值增量函數(shù)局部截斷誤差:整體截斷誤差:P階方法:局部截斷誤差顯式方法:隱式方法:右端項不包含右端項含有一階向后差商代替一階導(dǎo)數(shù)兩邊從

到積分,積分采用右矩形公式

從初值出發(fā)，計算m+1步后得到的值與準確值之差從準確值計算1步后得到的值與準確值之差CopyrightbyLiXinliang7改進的Euler方法僅迭代一次(n=1),將寫入第二式迭代求解

向前Euler方法向后Euler方法梯形法改進的Euler方法CopyrightbyLiXinliang8P階Taylor級數(shù)法……(*)式導(dǎo)數(shù)的計算比較復(fù)雜離散情形（非解析形式）設(shè)初值問題的解具有p+1階連續(xù)導(dǎo)數(shù)，利用Taylor公式(*)（很少使用）CopyrightbyLiXinliang9Runge-Kutta方法基本思想：函數(shù)在一點的導(dǎo)數(shù)值可以用該點附近若干點的函數(shù)值近似表示,因此可以把Taylor級數(shù)法中的增量函數(shù)改為

在這些點上函數(shù)值的組合,然后用Taylor展開確定待定系數(shù),使方法達到一定的階。N級Runge-Kutta方法:

是待定常數(shù)CopyrightbyLiXinliang10Runge-Kutta方法可定出:Euler公式將

在處展開

N階CopyrightbyLiXinliang11Runge-Kutta方法無窮多組解二階：

N階將

在處展開利用

CopyrightbyLiXinliang12Runge-Kutta方法無窮多組解N階2階格式滿足的條件：

（改進的Euler方法）（中點方法）（Heun方法）CopyrightbyLiXinliang13Runge-Kutta方法Butcher表自動滿足條件：新的與原來的在計算過程中有相同的增量CopyrightbyLiXinliang14Runge-Kutta方法將

在處展開CopyrightbyLiXinliang15Runge-Kutta方法CopyrightbyLiXinliang16Runge-Kutta方法三階三階Heun方法：三階Kutta方法：三階TVDRK方法(Shu-Osher)：6個系數(shù)，4個方程CopyrightbyLiXinliang17Runge-Kutta方法四階五階Kutta-Nystrom方法：經(jīng)典四階RK方法：N級RK方法能達到的最高階數(shù)顯式：10個系數(shù)，8個方程CopyrightbyLiXinliang18自適應(yīng)Runge-Kutta方法步長的選擇往往很困難一般地，要求數(shù)值解與精確解之差不超過某一誤差容限，由這個容限來確定數(shù)值方法中使用的步長（即不取固定步長）利用三階RK方法去估計二階RK方法的局部誤差，從而合理選取步長估計利用上式來選擇合理的步長以代替～假設(shè)為誤差容限，則可選取q，使得CopyrightbyLiXinliang19自適應(yīng)Runge-Kutta方法如果則如果重新計算CopyrightbyLiXinliang20高階方程的Runge-Kutta方法不把高階方程展成方程組，而是直接求解(已知)Taylor展開求待定系數(shù)一組解

CopyrightbyLiXinliang21多步法的構(gòu)造單步法：計算時只用到前一步的值

§10.2線性多步法

多步法（k步）：已知其中，是常數(shù)，：隱式理論上，有2k+1

個未知數(shù)，可獲得具有2k

階精度的格式CopyrightbyLiXinliang22多步法的構(gòu)造Taylor級數(shù)展開：p階精度

4階2步方法：

局部截斷誤差：其中由于不穩(wěn)定的原因，對于k>2，2k階的方法是沒意義的；存在k+1階及k+2(k偶數(shù))階穩(wěn)定的k步方法即可求得相應(yīng)系數(shù)如CopyrightbyLiXinliang23數(shù)值積分構(gòu)造多步法(Adams方法)以為插值節(jié)點構(gòu)造f(x,y(x))的k-1次Lagrange插值多項式（外插）余項：

Adams外插公式（顯式）Euler公式k步k階0123413-123-16555-5937-9點CopyrightbyLiXinliang24數(shù)值積分構(gòu)造多步法(Adams方法)以為插值節(jié)點構(gòu)造f(x,y(x))的k次Lagrange插值多項式（內(nèi)插）

Adams外插公式（顯式）Euler向前公式Adams內(nèi)插公式（隱式）Euler向后公式0123411158-1919-51點個點)階步(個點)步(階CopyrightbyLiXinliang25數(shù)值微分構(gòu)造多步法(Gear方法)以為插值節(jié)點構(gòu)造y(x)的k次Lagrange插值多項式余項：Gear方法（k步k階隱式）：記：

Taylor展開構(gòu)造012341124-1618-921248-3616-3CopyrightbyLiXinliang26四階Runge-Kutta方法四階Adams外插方法四階四步方法(Taylor)第一步計算需要用到。但初始條件只給出，因此需要用其它方法提供另外k-1個值。計算量穩(wěn)定性可能問題？多步法：RK方法計算4次函數(shù)值，多步法1次RK方法優(yōu)于多步法CopyrightbyLiXinliang27Admas外插方法：顯式

內(nèi)插方法：隱式預(yù)測-校正法預(yù)測公式（顯式）校正公式（隱式）Hamming

方法迭代CopyrightbyLiXinliang28Lipschitz條件：則稱該數(shù)值方法是收斂的。假設(shè)函數(shù)f(x,y)在區(qū)域D上連續(xù)，對區(qū)域D上的任意兩點，和，存在正的常數(shù)L(稱為Lipschitz常數(shù))，使得記則稱函數(shù)f(x,y)滿足Lipschitz條件。收斂：對中任意固定的點，對數(shù)值解，有相容：步長趨于0時，數(shù)值方法的截斷誤差也趨于0（至少是一階的）CopyrightbyLiXinliang29試驗方程定義穩(wěn)定：穩(wěn)定：存在正常數(shù)和，使得用任意步長以及任意兩個初值

和計算得到的兩組數(shù)值和，成立則稱該方法是穩(wěn)定的。誤差的增長有界對初值問題若一個數(shù)值方法在某一步由計算產(chǎn)生的誤差在以后各步中逐步減少，則稱該方法關(guān)于是穩(wěn)定的。記或：當時，，則稱用步長h的這個數(shù)值積分過程是穩(wěn)定的。穩(wěn)定區(qū)域：使數(shù)值方法穩(wěn)定的的集合。CopyrightbyLiXinliang30研究微分方程數(shù)值解的計算穩(wěn)定性問題：具有廣泛的代表性，線性常系數(shù)常微分方程組通過線性變換可得到這種情形線性非齊次微分方程組數(shù)值解的計算穩(wěn)定性與齊次方程組是一樣的變系數(shù)的常微分方程組，工程上常將系數(shù)固態(tài)化后研究某個時刻鄰域的性質(zhì)非線性的常微分方程組，在研究某個解鄰域的誤差時，可經(jīng)過線化，得到誤差所滿足的線性方程組試驗方程：CopyrightbyLiXinliang31Euler方法：從出發(fā)，假定在第m步按上式計算得到的值記為，在此基礎(chǔ)上計算下一步的值，可見，當且僅當時，第m步產(chǎn)生的誤差在以后各步中逐步減少，從而得到Euler方法的絕對穩(wěn)定區(qū)間為：穩(wěn)定CopyrightbyLiXinliang32p階顯式Runge-Kutta的穩(wěn)定性函數(shù)：隱式Runge-Kutta的穩(wěn)定性函數(shù)：RK方法的穩(wěn)定區(qū)域：

Euler:(-2,0)2-RK:(-2,0)3-RK:(-2.513,0)4-RK:(-2.785,0)顯式RK方法：CopyrightbyLiXinliang33一般k步多步法（改變形式）特征方程：k步多步法的穩(wěn)定區(qū)域：

令稱為線性k步法的特征多項式線性k步法穩(wěn)定(或稱零穩(wěn)定)的充分必要條件是滿足特征根條件：的所有根均在單位圓內(nèi)，并且在單位圓周上的根只能是單重根。CopyrightbyLiXinliang34穩(wěn)定區(qū)域的計算：試驗方程線性k步法特征方程對于給定的，若特征方程的所有根的模都小于1，則稱線性k步法關(guān)于絕對穩(wěn)定。若對所有，線性多步法都絕對穩(wěn)定，則稱為絕對穩(wěn)定區(qū)間。CopyrightbyLiXinliang35穩(wěn)定區(qū)域的計算：取，解出對應(yīng)于每個值的根，根的全體就是穩(wěn)定區(qū)域的邊界曲線。共軛性質(zhì)：Adams顯式格式1階(Euler):(-2,0)2階:(-1,0)3階:(-6/11,0)4階:(-3/10,0)Euler:(-2,0)2-RK:(-2,0)3-RK:(-2.513,0)4-RK:(-2.785,0)K步多步法:令Runge-Kutta方法:Adams隱式格式2階:(-∞,0)3階:(-6,0)4階:(-3,0)5階:(-90/49,0)Gear方法1-6階:(-∞,0)1-2階：A穩(wěn)定3-6階：非A穩(wěn)定穩(wěn)定區(qū)間RK格式不直接求，而是利用關(guān)系式：令①①②②②Adams顯式格式4階Gear不穩(wěn)定②CopyrightbyLiXinliang36在可以用常微分方程來描述的許多實際的物理或化學(xué)過程中，往往包含復(fù)雜的子過程及它們之間的相互作用，其中有的子過程表現(xiàn)為快變化的，而另一些相對來說是慢變化的，并且變化速度可以相差非常大的量級。相應(yīng)地，描述這些過程的常微分方程的解中也將包含快變分量和慢變分量。如果在一個過程中的快變子過程與慢變子過程的變化速度相差非常大，在數(shù)學(xué)上稱這種過程具有剛性(stiff)，而描述這種過程的常微分方程稱為剛性方程。(袁兆鼎等，剛性常微分方程初值問題的數(shù)值解法，科學(xué)出版社，1987)

§10.3剛性常微分方程的數(shù)值解法

病態(tài)方程或壞條件方程Stiffequationsareequationswherecertainimplicitmethods,inparticularBDF,performbetter,usuallytremendouslybetter,thanexplicitones.(Curtiss&Hirschfelder,1952;Hairer&Wanner,SolvingOrdinaryDifferentialEquationsII:StiffandDifferential-AlgebraticProblems,Springer,1996；科學(xué)出版社，2006)CopyrightbyLiXinliang37待求的向量函數(shù)已知的向量函數(shù)t是獨立變量，可以看成是時間假定矩陣A的Jordan標準型是對角陣，其特征值為（1）對應(yīng)的特征向量記為(1)的解有形式CopyrightbyLiXinliang38(1)的解有形式假定，即方程(1)是漸進穩(wěn)定的。當時，有，則稱為(1)的暫態(tài)解，為穩(wěn)態(tài)解暫態(tài)解可以表示成解分量的線性組合。實部確定的衰減特性，虛部確定這個量的振蕩特性對于一個穩(wěn)定系統(tǒng)，的實部一定是負的。稱為時間常數(shù)，用來表征的衰減速度。每經(jīng)過時刻，

衰減倍，約為1/3倍，越大衰減越快。振蕩頻率為

越大，振蕩越快。CopyrightbyLiXinliang39線性系統(tǒng)（1）稱作是剛性方程，如果有定義剛性方程的性質(zhì)(1)剛性方程是漸進穩(wěn)定的，解曲線從不同的初值都將趨向于它的穩(wěn)態(tài)解。各個解分量的衰減特性是不同的，衰減快的稱為快變分量，衰減慢的稱為慢變分量。剛性方程的解曲線將很快衰減到由慢變分量所確定的解上。剛性方程的解曲線可以分為二段，開始的一段為快變段，解曲線中的快變部分迅速衰減到可忽略的程度，將快變段稱作邊界層，其經(jīng)歷的時間記為，一般取快變分量衰減到原來的1/20時的時間，即。另一段是慢變段，或稱作邊界層外的段，它由較小的解分量來刻畫。(Lambert,1973)（r稱作剛性比）CopyrightbyLiXinliang40(2)剛性方程組具有奇異攝動的性質(zhì)。由于解曲線中的快變部分在邊界層內(nèi)很快衰減掉，在邊界層外，方程(1)的解曲線中所含的量的個數(shù)減少，使得解曲線的各個坐標之間不再是線性獨立的，而存在若干個代數(shù)關(guān)系式。利用這些關(guān)系式，可以用低階的方程組代替方程（1）。(3)剛性方程中矩陣A是病態(tài)的。剛性方程的性質(zhì)利用奇異攝動構(gòu)造解精確解如果J(t)在[a,b]上的變化充分小，則可用某個J(t0)來代替因此可用

的局部Jacobian矩陣的特征值來定義方程的剛性CopyrightbyLiXinliang41剛性性質(zhì)是數(shù)學(xué)問題本身的性質(zhì)，它不依賴于求解這個問題的數(shù)值方法，但正是由于這個性質(zhì)，使得傳統(tǒng)的常微分方程的數(shù)值積分方法遇到了極大的困難。由于穩(wěn)定性的要求，越大，選取的步長h越小。數(shù)值積分時，對應(yīng)于大的方程，其近似暫態(tài)解成分很快可以衰減到可以忽略的程度，即認為已達到對應(yīng)方程的穩(wěn)態(tài)解，往后的積分可以看成是求對應(yīng)于為小的穩(wěn)態(tài)解。此時，數(shù)值解的量值將由穩(wěn)態(tài)解的近似和對應(yīng)于為小的方程的近似解來確定。若要積分到穩(wěn)態(tài)解，則積分區(qū)間的長度的量級為。不管每個分量對最后的近似解是否起作用，積分過程中對每個都必須進行求解，因此所選取的步長仍應(yīng)使所有的均在穩(wěn)定區(qū)域內(nèi)，積分步長由對應(yīng)于為最大的分量來決定，步長的量級為總的積分步數(shù)的量級將不小于（剛性比，剛性大時，計算花費大，

甚至難以實現(xiàn)）CopyrightbyLiXinliang42另一種剛性定義（Shampine&Gear,1979）：矩陣A的所有特征值的實部不是很大的正數(shù)A至少有一個特征值的實部是很大的負數(shù)對應(yīng)于具有最大負實部的特征值的解分量變化是緩慢的A穩(wěn)定(Dahlquist,1963)：數(shù)值積分方法的穩(wěn)定區(qū)域含有復(fù)開左半平面，則稱該方法是A穩(wěn)定的。研究具有無限穩(wěn)定區(qū)域的方法顯式多步法（顯式Runge-Kutta方法）不可能是A穩(wěn)定的

A穩(wěn)定的隱式線性多步法不能超過2階剛性方程數(shù)值方法的發(fā)展：(1)減弱A穩(wěn)定性的限制，只需包含左半平面的一部分(2)保持A穩(wěn)定性，構(gòu)造新的數(shù)值積分方法四階向后差分（Gear）CopyrightbyLiXinliang43定理：顯式的k步公式不能是A穩(wěn)定的下面的多步公式是A穩(wěn)定的（1）向后Euler公式（2）梯形法（3）對于任何正整數(shù)k，存在相容的A穩(wěn)定k步公式定理：A穩(wěn)定的線性多步公式的階不得超過2階CopyrightbyLiXinliang44向后Euler方法穩(wěn)定區(qū)域利用向后Euler方法計算試驗方程穩(wěn)定原因:方法的穩(wěn)定區(qū)域包含了系統(tǒng)的不穩(wěn)定區(qū)域不穩(wěn)定穩(wěn)定存在的問題(1):方法的穩(wěn)定區(qū)域包含了系統(tǒng)的不穩(wěn)定區(qū)域Gear方法因此，構(gòu)造方法時，應(yīng)考慮被求解的系統(tǒng)是否穩(wěn)定，系統(tǒng)特征值的變化等從大的負值變化到大的正值不穩(wěn)定雙時間步方法即利用了二階Gear方法CopyrightbyLiXinliang45對于剛性很大的問題，A穩(wěn)定的格式仍可能得不到正確的結(jié)果，如梯形公式，時，非常接近于1。因此，一個希望的性質(zhì)是。L穩(wěn)定(Ehle，1969)：A穩(wěn)定，且當時，其穩(wěn)定性函數(shù)有。A穩(wěn)定條件不足如果一個隱式Runge-Kutta方法（A非奇異）滿足下列條件之一：(1)(2)則有，即該方法是L-穩(wěn)定的。對奇異攝動問題及微分-代數(shù)方程非常重要。存在的問題(2):CopyrightbyLiXinliang46A穩(wěn)定條件太強，許多方法雖然不是A穩(wěn)定，但也可以對很多問題進行計算

穩(wěn)定(Widlund,1967)：數(shù)值積分方法的穩(wěn)定區(qū)域含有如下復(fù)平面的集合，特別對多步方法有用。

穩(wěn)定(Cryer,1973)：數(shù)值積分方法的穩(wěn)定區(qū)域包含整個負實軸。研究穩(wěn)定的公式的性質(zhì)，可確定適合求解剛性方程的數(shù)值方法必備的一些條件許多實際的問題可能只具有實特征值A(chǔ)對任意的k，存在k步k階的穩(wěn)定的多步方法。如果高階多步法的誤差常數(shù)巨大，則不能用CopyrightbyLiXinliang47隱式Runge-Kutta方法

是待定常數(shù)Ehle(1968):s級的2s階隱式RK方法是A穩(wěn)定的。試驗方程試驗方程Cramer法則求Butcher表：求逆CopyrightbyLiXinliang48隱式Runge-Kutta的幾種構(gòu)造方法：基于Gauss求積公式（s級2s階）：s次Legendre多項式(1)求出s次Legendre多項式的s個0點，(2)由線性方程組求系數(shù)(3)計算系數(shù)矩陣基于Radau求積公式（s級2s-1階）：(1)求多項式的s個0點，，指定(2)計算系數(shù)(3)計算系數(shù)矩陣基于Lobatto求積公式（s級2s-2階）??四階三階A穩(wěn)定二階??00011/21/2??CopyrightbyLiXinliang49Runge-Kutta方法對角隱式RK方法：且至少有一個單對角隱式RK方法：11??10002/31/31/3??ImplicitEulerImplicitmidpointruleHammer&Hollingsworth（3階）??SDIRK，2級3階SDIRK，3級4階完全隱式對角隱式單對角隱式（SDIRK）顯式非A穩(wěn)定+：A穩(wěn)定-：非A穩(wěn)定CopyrightbyLiXinliang50方法的選取(1)對光滑的非剛性問題：要在大的時間區(qū)間上積分（精度要求需要的步數(shù)多），多步法可能是最好的(2)具有間斷的問題：應(yīng)選取單步方法，如RK方法(3)剛性問題：A穩(wěn)定的隱式RK方法（如DIRK）CopyrightbyLiXinliang5