克兰克-尼科尔森方法

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

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

## 方法

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

　
: ${\frac {\partial u}{\partial t}}=F\left(u,x,t,{\frac {\partial u}{\partial x}},{\frac {\partial ^{2}u}{\partial x^{2}}}\right)$

令 $u(i\Delta x,n\Delta t)=u_{i}^{n}\,$ ，則通過克蘭克－尼科爾森方法導出的差分方程是第 _n_ 步上採用前向歐拉方法與第 _n+1_ 步上採用後向歐拉方法的平均值（注意，克蘭克－尼科爾森方法本身不是這兩種方法簡單地取平均，方程對解隱式依賴）。

　
: ${\frac {u_{i}^{n+1}-u_{i}^{n}}{\Delta t}}=F_{i}^{n}\left(u,x,t,{\frac {\partial u}{\partial x}},{\frac {\partial ^{2}u}{\partial x^{2}}}\right)$ (前向歐拉方法)

　
: ${\frac {u_{i}^{n+1}-u_{i}^{n}}{\Delta t}}=F_{i}^{n+1}\left(u,x,t,{\frac {\partial u}{\partial x}},{\frac {\partial ^{2}u}{\partial x^{2}}}\right)$ (後向歐拉方法)

　
: ${\frac {u_{i}^{n+1}-u_{i}^{n}}{\Delta t}}={\frac {1}{2}}\left(F_{i}^{n+1}\left(u,x,t,{\frac {\partial u}{\partial x}},{\frac {\partial ^{2}u}{\partial x^{2}}}\right)+F_{i}^{n}\left(u,x,t,{\frac {\partial u}{\partial x}},{\frac {\partial ^{2}u}{\partial x^{2}}}\right)\right)$ (克蘭克－尼科爾森方法)

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

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

## 示例

### 線性擴散問題

* 線性擴散方程

　
: ${\frac {\partial u}{\partial t}}=a{\frac {\partial ^{2}u}{\partial x^{2}}}$

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

　
: ${\frac {u_{i}^{n+1}-u_{i}^{n}}{\Delta t}}={\frac {a}{2(\Delta x)^{2}}}\left((u_{i+1}^{n+1}-2u_{i}^{n+1}+u_{i-1}^{n+1})+(u_{i+1}^{n}-2u_{i}^{n}+u_{i-1}^{n})\right)$

引入變量 $r={\frac {a\Delta t}{2(\Delta x)^{2}}}$ :

　
: $-ru_{i+1}^{n+1}+(1+2r)u_{i}^{n+1}-ru_{i-1}^{n+1}=ru_{i+1}^{n}+(1-2r)u_{i}^{n}+ru_{i-1}^{n}\,$

這是一個三對角問題，應用三對角矩陣算法（追趕法）即可得到 $u_{i}^{n+1}$ ，而不需要對矩陣直接求逆。

* 準線性擴散方程

　
: ${\frac {\partial u}{\partial t}}=a(u){\frac {\partial ^{2}u}{\partial x^{2}}}$

離散化後則會得到非線性方程系統。但是某些情況下，通過使用 _a_ 的舊值，即用 $a_{i}^{n}(u)\,$ 替代 $a_{i}^{n+1}(u)\,$ ，可將問題線性化。其他時候，也可能在保證穩定性的基礎上使用顯式方法估計 $a_{i}^{n+1}(u)\,$

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

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

　
: $\langle 0\rangle {\frac {\partial C}{\partial t}}=D_{x}{\frac {\partial ^{2}C}{\partial x^{2}}}-U_{x}{\frac {\partial C}{\partial x}}-k(C-C_{N})-k(C-C_{M})$

其中 _C_ 表示污染物的富集水平，下標 _N_ 和 _M_ 分別對應上一通道和下一通道。

克蘭克－尼科爾森方法（i 對應位置，j 對應時間）將以上偏微分方程中的每個部分變換為

　
: $\langle 1\rangle {\frac {\partial C}{\partial t}}={\frac {C_{i}^{j+1}-C_{i}^{j}}{\Delta t}}$

　
: $\langle 2\rangle {\frac {\partial ^{2}C}{\partial x^{2}}}={\frac {1}{2(\Delta x)^{2}}}\left((C_{i+1}^{j+1}-2C_{i}^{j+1}+C_{i-1}^{j+1})+(C_{i+1}^{j}-2C_{i}^{j}+C_{i-1}^{j})\right)$

　
: $\langle 3\rangle {\frac {\partial C}{\partial x}}={\frac {1}{2}}\left({\frac {(C_{i+1}^{j+1}-C_{i-1}^{j+1})}{2(\Delta x)}}+{\frac {(C_{i+1}^{j}-C_{i-1}^{j})}{2(\Delta x)}}\right)$

　
: $\langle 4\rangle C={\frac {1}{2}}(C_{i}^{j+1}+C_{i}^{j})$

　
: $\langle 5\rangle C_{N}={\frac {1}{2}}(C_{Ni}^{j+1}+C_{Ni}^{j})$

　
: $\langle 6\rangle C_{M}={\frac {1}{2}}(C_{Mi}^{j+1}+C_{Mi}^{j})$

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

　
: $\lambda ={\frac {D_{x}\Delta t}{2\Delta x^{2}}}$

　
: $\alpha ={\frac {U_{x}\Delta t}{4\Delta x}}$

　
: $\beta ={\frac {k\Delta t}{2}}$

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

　
: $-\beta C_{Ni}^{j+1}-(\lambda +\alpha )C_{i-1}^{j+1}+(1+2\lambda +2\beta )C_{i}^{j+1}-(\lambda -\alpha )C_{i+1}^{j+1}-\beta C_{Mi}^{j+1}=\beta C_{Ni}^{j}+(\lambda +\alpha )C_{i-1}^{j}+(1-2\lambda -2\beta )C_{i}^{j}+(\lambda -\alpha )C_{i+1}^{j}+\beta C_{Mi}^{j}$

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

　
: $-(\lambda +\alpha )C_{i-1}^{j+1}+(1+2\lambda +\beta )C_{i}^{j+1}-(\lambda -\alpha )C_{i+1}^{j+1}-\beta C_{Mi}^{j+1}=+(\lambda +\alpha )C_{i-1}^{j}+(1-2\lambda -\beta )C_{i}^{j}+(\lambda -\alpha )C_{i+1}^{j}+\beta C_{Mi}^{j}$

同樣地， 最後一個通道只與前一個通道（_N_）有關聯，因此表達式可以簡化為

　
: $-\beta C_{Ni}^{j+1}-(\lambda +\alpha )C_{i-1}^{j+1}+(1+2\lambda +\beta )C_{i}^{j+1}-(\lambda -\alpha )C_{i+1}^{j+1}=\beta C_{Ni}^{j}+(\lambda +\alpha )C_{i-1}^{j}+(1-2\lambda -\beta )C_{i}^{j}+(\lambda -\alpha )C_{i+1}^{j}$

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

 $C_{0}^{j}$ : 當前時間步某通道的初始條件

 $C_{0}^{j+1}$ : 下一時間步某通道的初始條件

 $C_{N0}^{j}$ : 前一通道到當前時間步下某通道的初始條件

 $C_{M0}^{j}$ : 下一通道到當前時間步下某通道的初始條件

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

　
: ${\frac {\partial C}{\partial x}}_{x=z}={\frac {(C_{i+1}-C_{i-1})}{2\Delta x}}=0$

當且只當

　
: $C_{i+1}^{j+1}=C_{i-1}^{j+1}\,$

時，這一條件才被滿足。

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

　
: ${\begin{bmatrix}AA\end{bmatrix}}{\begin{bmatrix}C^{j+1}\end{bmatrix}}=[BB][C^{j}]+[d]$

其中，

　
: $\mathbf {C^{j+1}} ={\begin{bmatrix}C_{11}^{j+1}\\C_{12}^{j+1}\\C_{13}^{j+1}\\C_{14}^{j+1}\\C_{21}^{j+1}\\C_{22}^{j+1}\\C_{23}^{j+1}\\C_{24}^{j+1}\\C_{31}^{j+1}\\C_{32}^{j+1}\\C_{33}^{j+1}\\C_{34}^{j+1}\end{bmatrix}}$ $\mathbf {C^{j}} ={\begin{bmatrix}C_{11}^{j}\\C_{12}^{j}\\C_{13}^{j}\\C_{14}^{j}\\C_{21}^{j}\\C_{22}^{j}\\C_{23}^{j}\\C_{24}^{j}\\C_{31}^{j}\\C_{32}^{j}\\C_{33}^{j}\\C_{34}^{j}\end{bmatrix}}$

需要清楚的是，_AA_ 和 _BB_ 是由四個不同子矩陣組成的矩陣，

　
: $\mathbf {AA} ={\begin{bmatrix}AA1&AA3&0\\AA3&AA2&AA3\\0&AA3&AA1\end{bmatrix}}$

　
: $\mathbf {BB} ={\begin{bmatrix}BB1&-AA3&0\\-AA3&BB2&-AA3\\0&-AA3&BB1\end{bmatrix}}$

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

　
: $\mathbf {AA1} ={\begin{bmatrix}(1+2\lambda +\beta )&-(\lambda -\alpha )&0&0\\-(\lambda +\alpha )&(1+2\lambda +\beta )&-(\lambda -\alpha )&0\\0&-(\lambda +\alpha )&(1+2\lambda +\beta )&-(\lambda -\alpha )\\0&0&-2\lambda &(1+2\lambda +\beta )\end{bmatrix}}$

　
: $\mathbf {AA2} ={\begin{bmatrix}(1+2\lambda +2\beta )&-(\lambda -\alpha )&0&0\\-(\lambda +\alpha )&(1+2\lambda +2\beta )&-(\lambda -\alpha )&0\\0&-(\lambda +\alpha )&(1+2\lambda +2\beta )&-(\lambda -\alpha )\\0&0&-2\lambda &(1+2\lambda +2\beta )\end{bmatrix}}$

　
: $\mathbf {AA3} ={\begin{bmatrix}-\beta &0&0&0\\0&-\beta &0&0\\0&0&-\beta &0\\0&0&0&-\beta \end{bmatrix}}$

　
: $\mathbf {BB1} ={\begin{bmatrix}(1-2\lambda -\beta )&(\lambda -\alpha )&0&0\\(\lambda +\alpha )&(1-2\lambda -\beta )&(\lambda -\alpha )&0\\0&(\lambda +\alpha )&(1-2\lambda -\beta )&(\lambda -\alpha )\\0&0&2\lambda &(1-2\lambda -\beta )\end{bmatrix}}$   &

　
: $\mathbf {BB2} ={\begin{bmatrix}(1-2\lambda -2\beta )&(\lambda -\alpha )&0&0\\(\lambda +\alpha )&(1-2\lambda -2\beta )&(\lambda -\alpha )&0\\0&(\lambda +\alpha )&(1-2\lambda -2\beta )&(\lambda -\alpha )\\0&0&2\lambda &(1-2\lambda -2\beta )\end{bmatrix}}$

這裡的 _d_ 矢量用於保證邊界條件成立。在此示例中為 12x1 的矢量。

　
: $\mathbf {d} ={\begin{bmatrix}(\lambda +\alpha )(C_{10}^{j+1}+C_{10}^{j})\\0\\0\\0\\(\lambda +\alpha )(C_{20}^{j+1}+C_{20}^{j})\\0\\0\\0\\(\lambda +\alpha )(C_{30}^{j+1}+C_{30}^{j})\\0\\0\\0\end{bmatrix}}$

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

　
: ${\begin{bmatrix}C^{j+1}\end{bmatrix}}={\begin{bmatrix}AA^{-1}\end{bmatrix}}([BB][C^{j}]+[d])$

### 二維擴散問題

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

　
: ${\frac {\partial u}{\partial t}}=a\left({\frac {\partial ^{2}u}{\partial x^{2}}}+{\frac {\partial ^{2}u}{\partial y^{2}}}\right)$

假設網格滿足 $\Delta x=\Delta y$ 的特性，即可通過克蘭克－尼科爾森方法將得到離散方程

　
: ${\begin{aligned}u_{i,j}^{n+1}&=u_{i,j}^{n}+{\frac {1}{2}}{\frac {a\Delta t}{(\Delta x)^{2}}}{\big [}(u_{i+1,j}^{n+1}+u_{i-1,j}^{n+1}+u_{i,j+1}^{n+1}+u_{i,j-1}^{n+1}-4u_{i,j}^{n+1})\\&\qquad {}+(u_{i+1,j}^{n}+u_{i-1,j}^{n}+u_{i,j+1}^{n}+u_{i,j-1}^{n}-4u_{i,j}^{n}){\big ]}\end{aligned}}$

此方程可以再重組，配合柯朗數再進行簡化

　
: $\mu ={\frac {a\Delta t}{(\Delta x)^{2}}}.$

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

　
: ${\begin{aligned}&(1+2\mu )u_{i,j}^{n+1}-{\frac {\mu }{2}}\left(u_{i+1,j}^{n+1}+u_{i-1,j}^{n+1}+u_{i,j+1}^{n+1}+u_{i,j-1}^{n+1}\right)\\&\quad =(1-2\mu )u_{i,j}^{n}+{\frac {\mu }{2}}\left(u_{i+1,j}^{n}+u_{i-1,j}^{n}+u_{i,j+1}^{n}+u_{i,j-1}^{n}\right).\end{aligned}}$

## 應用在金融數學上

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

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

## 相關條目

* 金融數學
* 梯形法 (微分方程)

## 參考資料