2.2　常微分方程

海洋數(shù)值模型總是包含著對(duì)某種形式的偏微分方程的求解，許多模擬問(wèn)題需要求解常微分方程（ODEs）。1.15節(jié)中對(duì)溫鹽環(huán)流的區(qū)塊模擬的討論和無(wú)平流和擴(kuò)散效應(yīng)的生物生態(tài)系統(tǒng)模擬是典型的例子。這些問(wèn)題由一組常微分方程表示，這些方程通常是非線性的和耦合的。無(wú)論其中所包含方程的階數(shù)是多少，該組方程總是能夠簡(jiǎn)化到N個(gè)因變量的N個(gè)一階常微分方程。

在偏微分方程的求解過(guò)程中，利用有限差分法、有限元素法或者半光譜法進(jìn)行空間離散化后，得到的方程都是常微分方程。下文中將要討論的所有形式都與這個(gè)問(wèn)題有關(guān)，要對(duì)期望精度、內(nèi)存需求和CPU需求幾方面進(jìn)行一致考慮。由于在大多數(shù)海洋數(shù)值模型問(wèn)題中，用一個(gè)已知的初始值對(duì)物理系統(tǒng)進(jìn)行的改進(jìn)受到很大關(guān)注，這些問(wèn)題稱為初始值問(wèn)題。初始值問(wèn)題包括對(duì)系統(tǒng)的未來(lái)狀態(tài)進(jìn)行求解，也就是N個(gè)變量在有限時(shí)間t內(nèi)的值，給定初始狀態(tài)下，這些變量在t=0時(shí)的值。求解過(guò)程包括將時(shí)間分成許多小的步長(zhǎng)，然后使每一步的解得到改進(jìn)。所選形式的精度、效率以及穩(wěn)定性自然與所選擇的步長(zhǎng)有關(guān)系。

對(duì)耦合的非線性一階常微分方程組的初始值問(wèn)題進(jìn)行求解的久負(fù)盛名的方法是四階龍格庫(kù)塔方法（RK）。然而，還有很多高效的方法，例如高效自適應(yīng)步長(zhǎng)嵌入式五階RK方法和BS方法，這兩種方法在所采用的步長(zhǎng)下對(duì)截?cái)嗾`差大小進(jìn)行計(jì)算，然后在當(dāng)前步長(zhǎng)下用該信息對(duì)步長(zhǎng)進(jìn)行降低，從而達(dá)到預(yù)定的精度，并在下一步中使步長(zhǎng)最優(yōu)化，這些在Press（1992）的觀點(diǎn)中有介紹和描述。我們將對(duì)這些先進(jìn)方法的原理進(jìn)行簡(jiǎn)單描述，這里只對(duì)傳統(tǒng)的RK方法和Adama-Bashforth預(yù)測(cè)校正方法進(jìn)行描述，因?yàn)樗鼈儗?duì)大多數(shù)非重復(fù)性應(yīng)用來(lái)說(shuō)已經(jīng)足夠。

2.2.1　龍格庫(kù)塔方法

首先，對(duì)單一因變量y（t）的單一一階常微分方程進(jìn)行考慮，y（t）的形式為

在f（y，t）函數(shù)形式下，方程可以是線性的，也可以是非線性的。y（0）值是已知的，要對(duì)y（T）進(jìn)行求取。時(shí)間區(qū)間t=0到t=T被分成N個(gè)子區(qū)間（n=1…N+1），那么問(wèn)題就成為：yn已知的情況下，怎樣高效而精確地確定yn+1，時(shí)間步長(zhǎng)為Δt=tn+1-tn。經(jīng)典的顯式（或前向）歐拉方法只在初始點(diǎn)tn對(duì)式（2.2.1）的導(dǎo)數(shù)進(jìn)行計(jì)算，這一一階精度方法的精度和穩(wěn)定性都不是很高。四階RK方法除了要計(jì)算初始點(diǎn)處的導(dǎo)數(shù)外，還要計(jì)算兩次導(dǎo)數(shù)，一次是中點(diǎn)處的導(dǎo)數(shù)，還有一次是終點(diǎn)tn+1處的導(dǎo)數(shù)，對(duì)4個(gè)導(dǎo)數(shù)值進(jìn)行加權(quán)平均，用它來(lái)求tn到tn+1時(shí)的導(dǎo)數(shù)值為

需要注意的是，每一次新的導(dǎo)數(shù)估計(jì)都是在之前的最佳估計(jì)基礎(chǔ)上進(jìn)行的。如果函數(shù)f只包含獨(dú)立變量t的話，那么該方法就等價(jià)于對(duì)函數(shù)f（t）積分的經(jīng)典辛普森公式。每一步的截?cái)嗾`差是o（Δt5），但整個(gè)區(qū)間上的累計(jì)截?cái)嗾`差是o（Δt4）。盡管這種方法不是最高效的，但是準(zhǔn)確度高，而且魯棒性很高。擴(kuò)展為N方程系

這是很簡(jiǎn)單的（在溫鹽災(zāi)難問(wèn)題中的應(yīng)用可以參考本章2.1節(jié)，它包含兩個(gè)耦合的一階非線性常微分方程構(gòu)成的式子）。

如果取代使用任意預(yù)定的固定步長(zhǎng)的話，可以在求解過(guò)程中對(duì)步長(zhǎng)進(jìn)行自適應(yīng)控制，RK解的效率可以得到極大改善。自適應(yīng)步長(zhǎng)控制要求每一時(shí)間步長(zhǎng)下的求解中對(duì)截?cái)嗾`差進(jìn)行估計(jì)，這可以在每一步獲得兩個(gè)解得到，其中一個(gè)解是在滿步長(zhǎng)下得到，另一解是在半步長(zhǎng)下得到。因?yàn)榻財(cái)嗾`差與步長(zhǎng)的關(guān)系是已知的，兩個(gè)解（一個(gè)的分辨率為h，而另一個(gè)的分辨率為h/2）之間的差異能夠方便地對(duì)截?cái)嗾`差進(jìn)行估計(jì)，從而使得我們可以確定用什么樣的步長(zhǎng)來(lái)獲得達(dá)到預(yù)定精度的解。如果在當(dāng)前的一步中沒(méi)有得到期望精度，便可用所確定的最優(yōu)步長(zhǎng)進(jìn)行再求解，如果超過(guò)了，則可以在下一步中將步長(zhǎng)增加到最優(yōu)值。這一局部外推法的簡(jiǎn)單原理使我們能夠根據(jù)需要使用較小的步長(zhǎng)或者在步長(zhǎng)較小時(shí)達(dá)不到期望精度的情況下使用較大步長(zhǎng)，從而極大地改善計(jì)算效率。

另一種估計(jì)截?cái)嗾`差的方法是使用嵌入的RK公式，同樣可以用兩個(gè)函數(shù)來(lái)計(jì)算精度不同的兩個(gè)解。Fehlberg的方法（Press等，1992）利用一個(gè)五階RK公式和一個(gè)嵌入的四階RK公式的解之差來(lái)估計(jì)截?cái)嗾`差，然后用它自適應(yīng)地改變步長(zhǎng)以獲得期望精度。Press等（1992）對(duì)這個(gè)高效的RK方法及與其相關(guān)的代碼進(jìn)行了描述，這種方法用于RK解的效率是一個(gè)關(guān)鍵問(wèn)題的情況。

當(dāng)解的性質(zhì)未知的情況下，RK方法具有魯棒性，它是這種情形下的最佳方法。對(duì)于在積分區(qū)間上沒(méi)有任何奇異點(diǎn)的良態(tài)函數(shù)來(lái)說(shuō)，存在一種更為有效的積分方法，這些方法被稱為Bulirsch-Stoer方法（Press等，1992），它在每一時(shí)間步長(zhǎng)中都包含對(duì)零步長(zhǎng)的極限進(jìn)行外推，這是通過(guò)在每一步中用不同的步長(zhǎng)求取多個(gè)解而實(shí)現(xiàn)的，把這些解看成步長(zhǎng)的解析函數(shù)的話，就可以用一個(gè)多項(xiàng)式函數(shù)外推或一個(gè)有理函數(shù)外推來(lái)確定步長(zhǎng)無(wú)限小的情況下的解。外推最終得到函數(shù)的解及與之相關(guān)的誤差估計(jì)，然后用它來(lái)在每一步中增加解的個(gè)數(shù)而對(duì)解進(jìn)行改進(jìn)。很顯然，每一步得到的解的個(gè)數(shù)、所采用的步長(zhǎng)以及外推方法都決定著這種方法的效率。所采用的步長(zhǎng)可以用預(yù)先確定的序列得到，但是每一步得到的解的個(gè)數(shù)則由一個(gè)迭代過(guò)程確定，起始時(shí)為兩個(gè)，在迭代過(guò)程中對(duì)外推得到的解的期望誤差特征進(jìn)行檢查。如果對(duì)得到的解不滿意，那么用下一個(gè)（減小的）步長(zhǎng)對(duì)另一個(gè)解進(jìn)行求取，求解的方法使用理查德森外推法，并對(duì)誤差進(jìn)行重新檢查，一直重復(fù)這個(gè)過(guò)程，直到滿足期望誤差，或者終止，整個(gè)過(guò)程用整體上減小的步長(zhǎng)不斷重復(fù)。如上述的嵌入RK方法那樣，自適應(yīng)步長(zhǎng)控制方法被用來(lái)控制整個(gè)步長(zhǎng)。Press等（1992）也對(duì)此方法進(jìn)行了非常詳細(xì)的描述，并介紹了與之相關(guān)的代碼，在對(duì)整個(gè)域上解的奇異點(diǎn)的存在性產(chǎn)生懷疑并且對(duì)高效性有貢獻(xiàn)。

一種被成功應(yīng)用到地球物理問(wèn)題的時(shí)間離散化中的方法是三階Adams-Bashforth公式（Durran，1991；Haidvogel等，1997）

該方法克服了海洋和大氣模型中對(duì)偏微分方程進(jìn)行時(shí)間差分的二階RK方法和其他常用方法的一些局限性，它還可以用來(lái)作為預(yù)測(cè)校正公式的一部分，預(yù)測(cè)步長(zhǎng)下基于式（2.2.4）得到y(tǒng)n+1的第一個(gè)估計(jì)值，^表示第一個(gè)估計(jì)值，然后用它來(lái)在校正步驟中對(duì)yn+1的值進(jìn)行更新

這避免了對(duì)隱式方程進(jìn)行求解。在每一步中用該方程進(jìn)行迭代求解以獲得期望精度。對(duì)于包含時(shí)間導(dǎo)數(shù)的復(fù)雜表達(dá)式的問(wèn)題來(lái)說(shuō)，例如與時(shí)間相關(guān)的多維地球物理問(wèn)題，這些預(yù)測(cè)校正方法通常比Bulirsch-Stoer方法更好（Press等，1992）。然而，這種方法代價(jià)較高，因此也很少使用。

上面討論的方法是“顯式”形式的。當(dāng)因變量為兩個(gè)或更多時(shí)，就會(huì)產(chǎn)生問(wèn)題中所包含的時(shí)間尺度不止一個(gè)的可能性，在這種情況下，盡管所考慮的精度可能會(huì)允許一個(gè)較大的時(shí)間尺度或者這種形式可能會(huì)變得不穩(wěn)定，但為了在時(shí)間尺度最小的情況下進(jìn)行求解，時(shí)間步長(zhǎng)必須要足夠小。這些方程被稱為剛性方程，當(dāng)要保持穩(wěn)定性時(shí)，它需要允許使用較大步長(zhǎng)（只與期望精度相一致）的“隱式”形式。隱式的RK和BS形式可參考Press等（1992）。在此用一個(gè)簡(jiǎn)單的線性一階方程對(duì)此進(jìn)行說(shuō)明

當(dāng)Δt＞2/α?xí)r是不穩(wěn)定的，因?yàn)閨yn|隨時(shí)間步長(zhǎng)n單調(diào)增加，這種方法就是一個(gè)隱式的形式，它不計(jì)算步長(zhǎng)n處的導(dǎo)數(shù)，而是計(jì)算步長(zhǎng)n+1處的導(dǎo)數(shù)（后向歐拉法）

yn隨n單調(diào)遞減，因?yàn)樗c原始方程一致，無(wú)論步長(zhǎng)為多少。然而，對(duì)普通的方程組執(zhí)行隱式形式是很困難的，但利用對(duì)步長(zhǎng)n+1進(jìn)行求導(dǎo)的泰勒級(jí)數(shù)展開式的第一項(xiàng)能使這種執(zhí)行變得簡(jiǎn)單。例如，可以用這種半隱式方法對(duì)上述方程進(jìn)行求解（雖然該情形下是不需要的）

在當(dāng)前討論的情形中該方程所得結(jié)果與全隱式方法所得結(jié)果一致。

有時(shí)會(huì)存在初始點(diǎn)處因變量未知的情況，但是這種未知量很少，在終點(diǎn)（或中間點(diǎn)）的其余變量已知，這種問(wèn)題被稱為兩點(diǎn)邊界值問(wèn)題。例如，我們可能知道初始點(diǎn)的浮游植物濃度，而不知道初始點(diǎn)浮游動(dòng)物的濃度，只可能知道它在終點(diǎn)的濃度，利用初始值問(wèn)題解法估計(jì)未知量在初始點(diǎn)的值并獲得終點(diǎn)處的解，便很容易對(duì)這類問(wèn)題進(jìn)行求解。對(duì)估計(jì)值進(jìn)行迭代調(diào)整，直到這些未知量的解與其在終點(diǎn)的已知值相匹配時(shí)停止迭代，從而完成求解過(guò)程。然而，這不屬于本書的范疇，對(duì)于兩點(diǎn)邊界值問(wèn)題的有效解及軟件的描述，讀者可以參考Press等（1992）。