克蘭克-尼科爾森方法

克蘭克-尼科爾森方法(英語:Crank–Nicolson method)是一種數值分析有限差分法,可用於數值求解熱方程以及類似形式的偏微分方程[1]。它在時間方向上是隱式的二階方法,可以寫成隱式的龍格-庫塔法數值穩定。該方法誕生於20世紀,由約翰·克蘭克菲利斯·尼科爾森發展[2]

可以證明克蘭克-尼科爾森方法對於擴散方程(以及許多其他方程)是無條件穩定[3]。但是,如果時間步長Δt乘以熱擴散率,再除以空間步長平方Δx2的值過大(根據馮諾依曼穩定性分析,以大於1/2為準),近似解中將存在虛假的振盪或衰減。基於這個原因,當要求大時間步或高空間解像度的時候,往往會採用數值精確較差的後向歐拉法進行計算,這樣即可以保證穩定,又避免了解的偽振盪。

方法

編輯

克蘭克-尼科爾森方法在空間域上的使用中心差分;而時間域上應用梯形公式,保證了時間域上的二階收斂。例如,一維偏微分方程

 

 ,則通過克蘭克-尼科爾森方法導出的差分方程是第n步上採用前向歐拉方法與第n+1步上採用後向歐拉方法的平均值(注意,克蘭克-尼科爾森方法本身不是這兩種方法簡單地取平均,方程對解隱式依賴)。

  (前向歐拉方法)
  (後向歐拉方法)
  (克蘭克-尼科爾森方法)

對於F,通過中心差分方法使其在空間上是離散的。

注意,這是一個隱式方法,需要求解代數方程組以得到時間域上的下一個u值。如果偏微分方程是非線性的,中心差分後得到的方程依舊是非線性方程系統,因此在時間步上推進會涉及求解非線性代數方程組。許多問題中,特別是線性擴散,代數方程中的矩陣是三對角的,通過三對角矩陣算法可以高效求解,這樣,算法的時間複雜度由直接求解全矩陣的 轉化為 

示例

編輯

線性擴散問題

編輯
  • 線性擴散方程
 

通過克蘭克-尼科爾森方法將得到離散方程

 

引入變量 :

 

這是一個三對角問題,應用三對角矩陣算法(追趕法)即可得到 ,而不需要對矩陣直接求逆。

  • 準線性擴散方程
 

離散化後則會得到非線性方程系統。但是某些情況下,通過使用a的舊值,即用  替代 ,可將問題線性化。其他時候,也可能在保證穩定性的基礎上使用顯式方法估計 

一維多通道連接的擴散問題

編輯

這種模型可以用於描述水流中含穩定污染流,但只有一維信息的情況。它可以簡化為一維問題並得到有價值的信息。 可對水中污染溶質富集的問題進行建模,這種問題由三部分組成:已知的擴散方程( 為常量),平流分量(即由速度場導致的系統在空間上的變化,表示為常量Ux),以及與縱向通道k旁流的相互作用。

 

其中C表示污染物的富集水平,下標NM分別對應上一通道和下一通道。

克蘭克-尼科爾森方法(i對應位置,j對應時間)將以上偏微分方程中的每個部分變換為

 
 
 
 
 
 

現在引入以下常量用於簡化計算:

 
 
 

把 <1>, <2>, <3>, <4>, <5>, <6>, α, βλ 代入 <0>. 把新時間項(j+1)代入到左邊,當前時間項(j)代入到右邊,將得到

 

第一個通道只能與下一個通道(M)有關係,因此表達式可以簡化為:

 

同樣地, 最後一個通道只與前一個通道(N)有關聯,因此表達式可以簡化為

 

為求解此線性方程組,需要知道邊界條件在通道始端就已經給定了。

 : 當前時間步某通道的初始條件

 : 下一時間步某通道的初始條件

 : 前一通道到當前時間步下某通道的初始條件

 : 下一通道到當前時間步下某通道的初始條件

對於通道的末端最後一個節點,最方便的條件是是絕熱近似,則

 

當且只當

 

時,這一條件才被滿足。

以3個通道,5個節點為例,可以將線性系統問題表示為

 

其中,

   

需要清楚的是,AABB是由四個不同子矩陣組成的矩陣,

 
 

其中上述矩陣的的矩陣元對應於下一個矩陣和額外的4x4零矩陣。請注意,矩陣AABB的大小為12x12

 
 
 
    &  
 

這裏的d向量用於保證邊界條件成立。在此示例中為12x1的向量。

 

為了找到任意時間下污染物的聚集情況,需要對以下方程進行迭代計算:

 

二維擴散問題

編輯

將擴散問題延伸到二維的笛卡爾網格英語Cartesian grid,推導方程類似,但結果會是{{link-en|帶形矩陣|Banded matrix||的方程式,不是三角矩陣,二維的熱方程

 

假設網格滿足 的特性,即可通過克蘭克-尼科爾森方法將得到離散方程

 

此方程可以再重組,配合柯朗數英語Courant number再進行簡化

 

在克蘭克-尼科爾森方法下,不需要為了穩定性而限制柯朗數的上限,不過為了數值穩定度,柯朗數仍不能太高,可以將方程式重寫如下:

 

應用在金融數學上

編輯

許多的現象都可以用熱方程金融數學上稱為擴散方程)來建模,因此克蘭克-尼科爾森方法也可以用在這些領域中[4]。尤其金融衍生工具定價用的布萊克-休斯模型可以轉換為熱方程,因此期權定價數值解可以用克蘭克-尼科爾森方法求得。

因為期權定價若超過基本假設(例如改變股息)時,無法求得解析解,需要用上述方式求得。不過若是非平滑的最後條件(大部份的金融商品都是如此),克蘭克-尼科爾森方法會有數值的震盪,無法用濾波方式平緩。在期權定價上會反映在履約價Γ的變動。因此,一開始幾個步驟需要用其他比較不會震盪的方法(如全隱式有限差分法)。

相關條目

編輯

參考資料

編輯
  1. ^ Tuncer Cebeci. Convective Heat Transfer. Springer. 2002. ISBN 0-9668461-4-1. 
  2. ^ Crank, J.; Nicolson, P. A practical method for numerical evaluation of solutions of partial differential equations of the heat conduction type. Proc. Camb. Phil. Soc. 1947, 43 (1): 50–67. doi:10.1007/BF02127704. 
  3. ^ Thomas, J. W. Numerical Partial Differential Equations: Finite Difference Methods. Texts in Applied Mathematics 22. Berlin, New York: Springer-Verlag. 1995. ISBN 978-0-387-97999-1. . Example 3.3.2 shows that Crank–Nicolson is unconditionally stable when applied to  .
  4. ^ Wilmott, P.; Howison, S.; Dewynne, J. The Mathematics of Financial Derivatives: A Student Introduction. Cambridge Univ. Press. 1995 [2015-02-16]. ISBN 0-521-49789-2. (原始內容存檔於2012-10-19).