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

1、i.常微分方程初值問(wèn)題數(shù)值解法本章討論常微分方程初值問(wèn)題數(shù)值解法，主要是差分法。解微分方程的所謂差分法的要點(diǎn)如下：首先是區(qū)域的離散，即將連續(xù)的求解區(qū)域離散化成有限個(gè)網(wǎng)格點(diǎn)。其次是方程的離散， 例如用差商代替微商，或者對(duì)微分方程積分使之變成積分方程，然后數(shù)值積分，或者最后得到網(wǎng)格點(diǎn)上的近似解所滿足的一個(gè)差分方程，解之即得差分解。i.1常微分方程差分法考慮常微分方程初值問(wèn)題：求函數(shù) u(t)滿足竺=f (t,u), 0 : t 汀(i.1a)dtu(0) =u0(i.1b)其中f(t,u)是定義在區(qū)域G : 0蘭t蘭T, u £珀上的函數(shù)，u0和T是給定的常數(shù)。我們假 設(shè)f(t,u)對(duì)u

2、滿足Lipschitz條件，即存在常數(shù) L使得f (t,uj f(t,u2)蘭 L 5 -上,g0,T;比,比亡(亠嚴(yán))(i.2)這一條件保證了( i.1)的解是適定的，即存在，唯一，而且連續(xù)依賴于初值u0。通常情況下，(i.1)的精確解不可能用簡(jiǎn)單的解析表達(dá)式給出，只能求近似解。本章討 論常微分方程最常用的近似數(shù)值解法-差分方法。先來(lái)討論最簡(jiǎn)單的Euler法。為此，首先將求解區(qū)域0,T 離散化為若干個(gè)離散點(diǎn)：0 =t° : ti ： HI ： t”:：tN =T(i.3)其中tn =hn , h0稱為步長(zhǎng)。在微積分課程中我們熟知，微商(即導(dǎo)數(shù))是差商的極限。反過(guò)來(lái)，差商就是微商的近

3、似。在t =t0處，在(i.1a)中用向前差商u(tl)曲0)代替微商du，便得hdtu(tj =u(t°) hf(t°,u(t。)p如果忽略誤差項(xiàng)；0，再換個(gè)記號(hào)，用 口代替u(ti)便得到q -u° =hf(t0，u。)般地，我們有Euler 方法：un4t =Un +hf (tn,Un), n = 0,1,1", N -1 (i.4)從(i.1b)給出的初始值u0出發(fā)，由上式可以依次算出tJiltN上的差分解u1| ,uN。F面我們用數(shù)值積分法重新導(dǎo)出Euler法以及其它幾種方法。為此，在區(qū)間tn ,tn 1上積分常微分方程（i.1a），得tn J

4、LU（tn=U（tn）f（t,U（t）dt（ i.5）用各種數(shù)值積分公式計(jì)算（i.5 ）中的積分，便導(dǎo)致各種不同的差分法。例如，若用左矩形公 式就得到Euler法（i.4）。如果用右矩形公式，便得到下面的：隱式 Euler 方法：Un4i=Un+h f（tn*,Un）, n= 0,1，卅，N1 （ i.6）類(lèi)似地，如果用梯形公式，就得到改進(jìn)的Euler方法h2f(tn1，Un1)f 心 1，®,n=0,1川，N -1(i.7)當(dāng)f（t,U）關(guān)于U是非線性函數(shù)的時(shí)候，不能由（i.6）和（i.7）從Un直接算出Un1，稱這一類(lèi)方法為隱式，通常采用某種迭代法求解。例如，將一般的隱式方法寫(xiě)成

5、Un 1=F（tn,Un,Un J （i.8）則可以利用如下的迭代法由un算出Un 4 :u： = F(tn，Un，Uk + )，k =0,1川0Un+=Un(i.9)關(guān)于k的迭代通常只需進(jìn)行很少幾步就可以滿足精度要求了。為了避免對(duì)隱式方法進(jìn)行迭代的麻煩，比如說(shuō)對(duì)于改進(jìn)的Euler方法（i.7），可以采用某種預(yù)估法近似算出f （tn 1,un，然后再用（i .7）作校正，這就導(dǎo)致所謂預(yù)估校正法。下面給出一個(gè)例子:預(yù)估:| un = f(tn,un)Un 1 = Unj 2hU nU n 1 _ f(tn-1,Un 1)(i.10)校正:Un1 f 2(Un1 Un)這是一個(gè)多步法，即計(jì)算節(jié)點(diǎn)t

6、n 1上的近似值Un 4時(shí)，除了用到前一點(diǎn)的近似值Un之外，還要用到un4，甚至可能用到un/,HI。而用前面的各種 Euler法計(jì)算節(jié)點(diǎn)tn1上的近似值 比計(jì)時(shí)，只用到un,因此稱之為單步法。F面給出另一個(gè)多步法的例子。在區(qū)間出心七上積分（i.1a），得U（tn .2） =U（tn）J f （t,U（t）dtLn用Simps on公式（即把被積函數(shù)看作二次函數(shù)）近似計(jì)算積分，便得到Mil ne 法:(i.11)un2 二 un h(fn 4fn 1 fn 2), n =0,1,”|,N -2 3用多步法（i.10 ）或（i.11）計(jì)算時(shí)，必須先用某種單步法由u0計(jì)算出u1，稱為造表頭。然后再

7、逐次算出U2,U3,|,UN。一般說(shuō)來(lái)，多步法比Euler法等簡(jiǎn)單的單步法精度要高一些。下面我們討論一類(lèi)所謂Ronge-Kutta法。他們是單步法，但是其精度可以與多步法比美。最常用的是下面的標(biāo)準(zhǔn)Ronge-Kutta 法和 Gear 法：Ki = f (tn，Un)hKi、K2 二 f (tn * Nn 2 )標(biāo)準(zhǔn) Lun ge-Kutta 法:h 2+nc-3K(i.12)9 / 7K4 = f (tnh,UK3)hUn i °(Ki 2& 2K3 K4). 6K1 = f (tn, un )K2 = f (tn 十2，山 +4Gea方法:(i.13)彳 & =

8、f (tn +£,叫 +( ft 島心 *(1-£)K2)K4 = f (J 卄心 +(-£)心 +(1+£)心) 叫出=山+£(心+2(1島心+2(1 + £)厶+ ©)從幾何上，Ronge-Kutta法可以粗略地解釋為：在區(qū)間tn,tn1中選取若干個(gè)點(diǎn)（可以重復(fù)）.1二tn2 -川-k二tn 1，僅僅利用在區(qū)間 出乙內(nèi)可以得到的所有信息，依 次給出函數(shù)f （t,u（t）在這些點(diǎn)上盡可能精確的的近似值K1,K2l,Kk，然后把它們組合起來(lái)，盡可能精確地近似計(jì)算（i.5）中的積分。Ronge-Kutta法一般的構(gòu)造方式如下。

9、選定常數(shù)k ,令K f（tn,un）K2 二 f（tn-hU亦1）K3 = f (tn，二込叫 * -31K1 -32 K2)Kk =f(tn ：kh,Un *Ki 七K2 III ",kKk)Un 1 =6 h( "Ki 上心川K)(i.14)選取這些待定常數(shù):i，'j-m的原則是：將(i.14)在(tn,Un)作Taylor展開(kāi)，然后按照h的 冪重新整理，使得1213Un 1二比 /2h3h ' ( "5)2!31與微分方程(i.1a)的解u =u(t)在t =tn處的Taylor展式1 2 1 3U(tn 1) =U(tn) 帥可時(shí)川(6)有

10、盡可能多的項(xiàng)重合，即要求1 二 fn,2 二仁 3 二 fn ,川這里花二f(tn,Un),仁df (t,u)dt等等。(t,U)蟲(chóng)“叫)(i.18)(i.17)按照(i.14)構(gòu)造出的都是顯式Runge-Kutta方法,在每一個(gè) &的表達(dá)式中只出現(xiàn) Un。如果允許在某一個(gè) Kj的表達(dá)式中出現(xiàn)un彳，則可以導(dǎo)出隱式 Runge-Kutta方法。i.2常微分方程組與高階常微分方程先來(lái)考慮下面的常微分方程組初值問(wèn)題dJ=f1(t,U1,U2l|,Um)dt如=f 2(t,U1,U2,川，Um)dtlllllldJ號(hào)=fm(t,U1,U2,川,Um)dtu1(0u；,um(0)=u：利用向量

11、記號(hào)，上式可以改寫(xiě)為u (t) =f (t,u).u(0) =Uo上節(jié)中各方法都可以直接應(yīng)用到常微分方程組(i.18)。例如，Euler方法成為UU nhf (tn, Un), n =0,1,川,N -1再來(lái)考慮高階常微分方程dmududmJu(i.19)廿二(t,u,不,川,耐)，t 0 d mJudu.mA (0) = Vm , ： 11 7 (0) = V2 , U (0) = vidtdt這時(shí)，可以令/ 1 2m、 /u = (u ,u , ,u ) =(u,dudtdmJu，川討(i.20)f = (f1, f2, II, fm) =(u2,u3, |lLum, f(t,u1,u2J

12、H,um) (i.21)U0 =(V1,V2,|),Vm)(i.22)于是可以把高階常微分方程(i.19)化成一階常微分方程組(i.18)。i.3收斂性與穩(wěn)定性截?cái)嗾`差粗略地說(shuō),截?cái)嗾`差可以定義為將微分方程解帶入到差分方程后得到的誤差，代表了微分方程與差分方程之間的誤差。例如，由Taylor展式和微分方程(i.1a)得到2!2! dt其中n是區(qū)間tn,tn 1上某個(gè)常數(shù)。與Euler法ununhf(tn,un)相比較，定義余項(xiàng)df (t,u(t) 為Euler法的截?cái)嗾`差，它關(guān)于h是2階的，記為0(f)。將上節(jié)中討論2! dt t的單步法寫(xiě)成一般形式F(un,un1,h) =0，則可以定義截?cái)?/p>

13、誤差為F(u(tn),u(tn1),h)。對(duì)1于每一個(gè)差分法在適當(dāng)?shù)狞c(diǎn)(例如 t =tn，tn-1或者(tn tn 1)作Taylor展開(kāi)，就可以得2到截?cái)嗾`差的階。對(duì)多步法可以類(lèi)似處理。上節(jié)中各方法的截?cái)嗾`差階分別為：表i.1常微分方程差分法的截?cái)嗾`差階差分法Euler隱式Euler改進(jìn)Euler預(yù)估校正(i.10)Milne標(biāo)準(zhǔn) Ronge-KuttaGear截?cái)嗾`差o(h2)o(h2)o(h3)o(h3)0(h5)0(h5)o(h5)相容性 一個(gè)差分方法稱為相容的，如果其截?cái)嗾`差至少是一階的。收斂性實(shí)用中我們更關(guān)心的是近似解的收斂性，即卩h > 0時(shí)，是否有un > u(tn

14、)。在適當(dāng)?shù)臈l件下，例如步長(zhǎng)h足夠小，右端函數(shù)f和解u足夠光滑等等，可以證明以上討論的各方法都是收斂的，并且有估計(jì)式U(tn) -U<Chp(i.23)其中常數(shù)C與u和f有關(guān)，與h和n無(wú)關(guān)；階數(shù)p等于截?cái)嗾`差的階數(shù)減去1。穩(wěn)定性收斂性考慮的是差分方程的精確解與微分方程的精確解之間的誤差。但是在計(jì)算機(jī)上求差分解時(shí)，由于計(jì)算機(jī)字長(zhǎng)的限制，不可避免地產(chǎn)生舍入誤差。另外，方程中的系數(shù)和 初值等等由于測(cè)量條件等限制，也會(huì)產(chǎn)生數(shù)據(jù)誤差。而穩(wěn)定性研究的就是舍入誤差和數(shù)據(jù)誤差對(duì)差分方程計(jì)算結(jié)果的影響，即計(jì)算過(guò)程中某一步的 “差之毫厘”會(huì)不會(huì)導(dǎo)致后面結(jié)果的“失之千里”。根據(jù)誤差來(lái)源和誤差衡量標(biāo)準(zhǔn)的不同，可

15、以定義形形色色的穩(wěn)定性。一般說(shuō)來(lái)，隱式方法的穩(wěn)定性好于同類(lèi)型的顯式方法（例如Euler方法與隱式Euler方法）。另外，為了保證穩(wěn)定性，常常要求步長(zhǎng)h足夠小。對(duì)于常微分方程數(shù)值解來(lái)說(shuō)，最簡(jiǎn)單常用的是關(guān)于初值的穩(wěn)定性。 以單步法為例。稱某差分 法關(guān)于初值穩(wěn)定，如果對(duì)同一個(gè)微分方程和所有足夠小的h，存在常數(shù)C，使得從不同初值u0和V。出發(fā)的兩個(gè)差分解 CuJ和：vn /之間的誤差滿足maxn Un v. EC u。一v°(i.24)其中N二T/h。由于任意一對(duì)Vn都可以看作是初值，關(guān)于初值的穩(wěn)定性其實(shí)是考察如下問(wèn)題：假設(shè)某一步計(jì)算有誤差， 而其后的計(jì)算不再有誤差, 的影響如何。本章所討論

16、的所有差分方法都是關(guān)于初值穩(wěn)定的。 時(shí)，通常指的是關(guān)于初值穩(wěn)定。不難證明，如果一個(gè)差分方法是穩(wěn)定且相容的，則一定收斂，即那么這一步的誤差對(duì)以后結(jié)果不加區(qū)別地說(shuō)差分方法穩(wěn)定口代）一叫T 0,當(dāng)nT嗆(i.25)du dt可以認(rèn)為這個(gè)模型方程的真解和差分解的性質(zhì)代表了非線性常微分方程dudtf (t,u)在某因而，我們進(jìn)一步考察差分方法的所在差分解的實(shí)際計(jì)算中， 每步都難以避免地舍入誤差。謂絕對(duì)穩(wěn)定性，即每一步的舍入誤差積累起來(lái)，是否會(huì)對(duì)后面的運(yùn)算結(jié)果產(chǎn)生太壞的影響?？疾炷骋粋€(gè)差分法對(duì)非線性常微分方程的絕對(duì)穩(wěn)定性是十分困難的事情，無(wú)法得到象對(duì)初值穩(wěn)定性那樣的漂亮結(jié)論。通常的做法是考慮模型方程(i.

17、26)一個(gè)局部的真解和差分解的性質(zhì)。（例如在（t, u）二（tn,un）處將f （t, u）線性化，將-看作f（tn,un），忽略其他項(xiàng)。）記h二。關(guān)于某一個(gè)差分法的絕對(duì)穩(wěn)定性的典型結(jié)果是：給:u出一個(gè)區(qū)間（，J，使得對(duì)于所有的 h （：），求解（i.25）的這個(gè)差分法都是絕對(duì)穩(wěn)定的。這個(gè)區(qū)間（）稱之為這個(gè)差分法的絕對(duì)穩(wěn)定域。絕對(duì)穩(wěn)定域越大，h的取值范圍就越大，的允許區(qū)值就越多，所代表的非線性常微分方程也就越多。表i.2常微分方程差分法的絕對(duì)穩(wěn)定域差分法Euler隱式Euler改進(jìn)EulerMilne標(biāo)準(zhǔn)R-KGear絕對(duì)穩(wěn)定域（-2，0）(皿,0)(£0)(0,0)(-2.78,0)(-2.78,0)關(guān)于絕對(duì)穩(wěn)定性的一個(gè)必要條件是：若一個(gè)多步法是相容并且穩(wěn)定的，則它絕對(duì)穩(wěn)定（即 它的絕對(duì)穩(wěn)定域非空）的一個(gè)必要條件是：0。這相當(dāng)于僅當(dāng) f : 0時(shí)，多步法才可能絕對(duì)穩(wěn)定。而當(dāng) f 0時(shí)，在差分計(jì)算中誤差增長(zhǎng)很快。那么這時(shí)候是否就無(wú)法用差分法du近似計(jì)算微分方程了呢？