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

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

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

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

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

知識(shí)點(diǎn)：Runge-Kutta方法

線性多步法

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

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

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

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

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

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

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

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

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

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

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

在處展開

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

N階將

在處展開利用

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

（改進(jìn)的Euler方法）（中點(diǎn)方法）（Heun方法）CopyrightbyLiXinliang13Runge-Kutta方法Butcher表自動(dòng)滿足條件：新的與原來的在計(jì)算過程中有相同的增量CopyrightbyLiXinliang14Runge-Kutta方法將

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

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

§10.2線性多步法

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

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

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

4階2步方法：

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

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

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

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

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

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

和計(jì)算得到的兩組數(shù)值和，成立則稱該方法是穩(wěn)定的。誤差的增長(zhǎng)有界對(duì)初值問題若一個(gè)數(shù)值方法在某一步由計(jì)算產(chǎn)生的誤差在以后各步中逐步減少，則稱該方法關(guān)于是穩(wěn)定的。記或：當(dāng)時(shí)，，則稱用步長(zhǎng)h的這個(gè)數(shù)值積分過程是穩(wěn)定的。穩(wěn)定區(qū)域：使數(shù)值方法穩(wěn)定的的集合。CopyrightbyLiXinliang30研究微分方程數(shù)值解的計(jì)算穩(wěn)定性問題：具有廣泛的代表性，線性常系數(shù)常微分方程組通過線性變換可得到這種情形線性非齊次微分方程組數(shù)值解的計(jì)算穩(wěn)定性與齊次方程組是一樣的變系數(shù)的常微分方程組，工程上常將系數(shù)固態(tài)化后研究某個(gè)時(shí)刻鄰域的性質(zhì)非線性的常微分方程組，在研究某個(gè)解鄰域的誤差時(shí)，可經(jīng)過線化，得到誤差所滿足的線性方程組試驗(yàn)方程：CopyrightbyLiXinliang31Euler方法：從出發(fā)，假定在第m步按上式計(jì)算得到的值記為，在此基礎(chǔ)上計(jì)算下一步的值，可見，當(dāng)且僅當(dāng)時(shí)，第m步產(chǎn)生的誤差在以后各步中逐步減少，從而得到Euler方法的絕對(duì)穩(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步法的特征多項(xiàng)式線性k步法穩(wěn)定(或稱零穩(wěn)定)的充分必要條件是滿足特征根條件：的所有根均在單位圓內(nèi)，并且在單位圓周上的根只能是單重根。CopyrightbyLiXinliang34穩(wěn)定區(qū)域的計(jì)算：試驗(yàn)方程線性k步法特征方程對(duì)于給定的，若特征方程的所有根的模都小于1，則稱線性k步法關(guān)于絕對(duì)穩(wěn)定。若對(duì)所有，線性多步法都絕對(duì)穩(wěn)定，則稱為絕對(duì)穩(wěn)定區(qū)間。CopyrightbyLiXinliang35穩(wěn)定區(qū)域的計(jì)算：取，解出對(duì)應(yīng)于每個(gè)值的根，根的全體就是穩(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在可以用常微分方程來描述的許多實(shí)際的物理或化學(xué)過程中，往往包含復(fù)雜的子過程及它們之間的相互作用，其中有的子過程表現(xiàn)為快變化的，而另一些相對(duì)來說是慢變化的，并且變化速度可以相差非常大的量級(jí)。相應(yīng)地，描述這些過程的常微分方程的解中也將包含快變分量和慢變分量。如果在一個(gè)過程中的快變子過程與慢變子過程的變化速度相差非常大，在數(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是獨(dú)立變量，可以看成是時(shí)間假定矩陣A的Jordan標(biāo)準(zhǔn)型是對(duì)角陣，其特征值為（1）對(duì)應(yīng)的特征向量記為(1)的解有形式CopyrightbyLiXinliang38(1)的解有形式假定，即方程(1)是漸進(jìn)穩(wěn)定的。當(dāng)時(shí)，有，則稱為(1)的暫態(tài)解，為穩(wěn)態(tài)解暫態(tài)解可以表示成解分量的線性組合。實(shí)部確定的衰減特性，虛部確定這個(gè)量的振蕩特性對(duì)于一個(gè)穩(wěn)定系統(tǒng)，的實(shí)部一定是負(fù)的。稱為時(shí)間常數(shù)，用來表征的衰減速度。每經(jīng)過時(shí)刻，

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

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

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

甚至難以實(shí)現(xiàn)）CopyrightbyLiXinliang42另一種剛性定義（Shampine&Gear,1979）：矩陣A的所有特征值的實(shí)部不是很大的正數(shù)A至少有一個(gè)特征值的實(shí)部是很大的負(fù)數(shù)對(duì)應(yīng)于具有最大負(fù)實(shí)部的特征值的解分量變化是緩慢的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）對(duì)于任何正整數(shù)k，存在相容的A穩(wěn)定k步公式定理：A穩(wěn)定的線性多步公式的階不得超過2階CopyrightbyLiXinliang44向后Euler方法穩(wěn)定區(qū)域利用向后Euler方法計(jì)算試驗(yàn)方程穩(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)造方法時(shí)，應(yīng)考慮被求解的系統(tǒng)是否穩(wěn)定，系統(tǒng)特征值的變化等從大的負(fù)值變化到大的正值不穩(wěn)定雙時(shí)間步方法即利用了二階Gear方法CopyrightbyLiXinliang45對(duì)于剛性很大的問題，A穩(wěn)定的格式仍可能得不到正確的結(jié)果，如梯形公式，時(shí)，非常接近于1。因此，一個(gè)希望的性質(zhì)是。L穩(wěn)定(Ehle，1969)：A穩(wěn)定，且當(dāng)時(shí)，其穩(wěn)定性函數(shù)有。A穩(wěn)定條件不足如果一個(gè)隱式Runge-Kutta方法（A非奇異）滿足下列條件之一：(1)(2)則有，即該方法是L-穩(wěn)定的。對(duì)奇異攝動(dòng)問題及微分-代數(shù)方程非常重要。存在的問題(2):CopyrightbyLiXinliang46A穩(wěn)定條件太強(qiáng)，許多方法雖然不是A穩(wěn)定，但也可以對(duì)很多問題進(jìn)行計(jì)算

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

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

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