克蘭克-尼科爾森方法

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

可以證明克蘭克－尼科爾森方法對於擴散方程（以及許多其他方程）是無條件穩定

$u_t=au_{xx}$。但是，如果時間步長Δ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} - 2 u_{i}^{n + 1} + u_{i - 1}^{n + 1}) + (u_{i + 1}^{n} - 2 u_{i}^{n} + u_{i - 1}^{n}) \right)$$

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

$$-r u_{i + 1}^{n + 1} + (1 + 2 r)u_{i}^{n + 1} - r u_{i - 1}^{n + 1} = r u_{i + 1}^{n} + (1 - 2 r)u_{i}^{n} + r u_{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} - 2 C_{i}^{j + 1} + C_{i - 1}^{j + 1}) + (C_{i + 1}^{j} - 2 C_{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)} + <pre> \frac{(C_{i + 1}^{j} - C_{i - 1}^{j})}{2 (\Delta x)} </pre> \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}$$

把 , , , , , , α, β 和 λ 代入 . 把新時間項(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])$$

### 二維擴散問題

方程式，不是三角矩陣，二維的熱方程

$$\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{align}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{align}$$

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

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

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

$$\begin{align}&(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{align}$$

## 應用在金融數學上

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

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