查看“︁克兰克-尼科尔森方法”︁的源代码

'''克兰克－尼科尔森方法'''（{{lang-en|Crank–Nicolson method}}）是一種[[数值分析]]的[[有限差分法]]，可用于数值求解[[热方程]]以及类似形式的[[偏微分方程]]<ref>{{cite book | title = Convective Heat Transfer | author = Tuncer Cebeci | publisher = Springer | year = 2002 | isbn = 0-9668461-4-1 | url = http://books.google.com/?id=xfkgT9Fd4t4C&pg=PA257&dq=%22Crank-Nicolson+method%22 }}</ref>。它在时间方向上是[[顯式和隱式方法|隐式]]的二阶方法，可以寫成隐式的[[龍格－庫塔法]]，[[数值稳定性|数值稳定]]。该方法诞生于20世纪，由[[約翰·克蘭克]]与[[菲利斯·尼科爾森]]发展<ref>{{Cite journal 
 | title = A practical method for numerical evaluation of solutions of partial differential equations of the heat conduction type
 | journal = Proc. Camb. Phil. Soc.
 | volume = 43 
 | issue = 1
 | year = 1947
 | pages = 50–67
 | doi = 10.1007/BF02127704 
 | last1 = Crank 
 | first1 = J. 
 | last2 = Nicolson 
 | first2 = P. 
 | postscript = .
}}</ref>。

可以证明克兰克－尼科尔森方法对于[[扩散方程]]（以及许多其他方程）是无条件[[数值稳定性|稳定]]<ref>{{Cite book | last1=Thomas | first1=J. W. | title=Numerical Partial Differential Equations: Finite Difference Methods | publisher=Springer-Verlag | location=Berlin, New York | series=Texts in Applied Mathematics | isbn=978-0-387-97999-1 | year=1995 |volume=22 | postscript=. }}. Example 3.3.2 shows that Crank–Nicolson is unconditionally stable when applied to <math>u_t=au_{xx}</math>.</ref>。但是，如果时间步长&Delta;{{var|t}}乘以[[熱擴散率]]，再除以空间步长平方&Delta;{{var|x}}{{sup|2}}的值过大（根據[[馮諾依曼穩定性分析]]，以大于1/2為準），近似解中将存在虚假的振荡或衰减。基于这个原因，当要求大时间步或高空间分辨率的时候，往往会采用数值精确较差的[[后向欧拉法]]进行计算，这样即可以保证稳定，又避免了解的伪振荡。

== 方法 ==
克兰克－尼科尔森方法在空间域上的使用中心差分；而时间域上应用梯形公式，保证了时间域上的二阶收敛。例如，一维偏微分方程

:<math>\frac{\partial u}{\partial t} = F\left(u, x, t, \frac{\partial u}{\partial x}, \frac{\partial^2 u}{\partial x^2}\right)</math>

令<math>u(i \Delta x, n \Delta t) = u_{i}^{n}\,</math>，则通过克兰克－尼科尔森方法导出的差分方程是第''n''步上采用前向欧拉方法与第''n+1''步上采用后向欧拉方法的平均值（注意，克兰克－尼科尔森方法本身不是这两种方法简单地取平均，方程对解隐式依赖）。

:<math>\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)</math>               (前向欧拉方法)
:<math>\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)</math>               (后向欧拉方法)
:<math>\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)</math>        (克兰克－尼科尔森方法)

对于''F''，通过中心差分方法使其在空间上是离散的。

注意，这是一个隐式方法，需要求解代数方程组以得到时间域上的下一个''u''值。如果偏微分方程是非线性的，中心差分后得到的方程依旧是非线性方程系统，因此在时间步上推进会涉及求解非线性代数方程组。许多问题中，特别是线性扩散，代数方程中的矩阵是三对角的，通过[[三对角矩阵算法]]可以高效求解，这样，算法的时间复杂度由直接求解全矩阵的<math>\mathcal{O}(n^3)</math>转化为<math>\mathcal{O}(n)</math>。

== 示例 ==
=== 线性扩散问题 ===

* 线性扩散方程

:<math>\frac{\partial u}{\partial t} = a \frac{\partial^2 u}{\partial x^2}</math>

通过克兰克－尼科尔森方法将得到离散方程

:<math>\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)</math>

引入变量<math>r = \frac{a \Delta t}{2 (\Delta x)^2}</math>:

:<math>-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}\,</math>

这是一个三对角问题，应用三对角矩阵算法（追赶法）即可得到<math>u_{i}^{n+1}</math>，而不需要对矩阵直接求逆。

* 准线性扩散方程

:<math>\frac{\partial u}{\partial t} = a(u) \frac{\partial^2 u}{\partial x^2}</math>

离散化后则会得到非线性方程系统。但是某些情况下，通过使用''a''的旧值，即用<math>a_{i}^{n}(u)\,</math> 替代<math>a_{i}^{n + 1}(u)\,</math>，可将问题线性化。其他时候，也可能在保证稳定性的基础上使用显式方法估计<math>a_{i}^{n + 1}(u)\,</math>

=== 一维多通道连接的扩散问题 ===
 
这种模型可以用于描述水流中含稳定污染流，但只有一维信息的情况。它可以简化为一维问题并得到有价值的信息。
可对水中污染溶质富集的问题进行建模，这种问题由三部分组成：已知的扩散方程（<math>D_x</math>为常量），平流分量（即由速度场导致的系统在空间上的变化，表示为常量''Ux''），以及与纵向通道k旁流的相互作用。

:<math>\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)</math>

其中''C''表示污染物的富集水平，下标''N''和''M''分别对应上一通道和下一通道。

克兰克－尼科尔森方法（i对应位置，j对应时间）将以上偏微分方程中的每个部分变换为

:<math>\langle 1 \rangle \frac{\partial C}{\partial t} = \frac{C_{i}^{j + 1} - C_{i}^{j}}{\Delta t}</math>

:<math>\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)</math>

:<math>\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)</math>

:<math> \langle 4 \rangle C= \frac{1}{2} (C_{i}^{j+1} + C_{i}^{j})</math>

:<math> \langle 5 \rangle C_N= \frac{1}{2} (C_{Ni}^{j+1} + C_{Ni}^{j})</math>

:<math> \langle 6 \rangle C_M= \frac{1}{2} (C_{Mi}^{j+1} + C_{Mi}^{j})</math>

现在引入以下常量用于简化计算：

:<math> \lambda= \frac{D_x\Delta t}{2 \Delta x^2}</math>

:<math> \alpha= \frac{U_x\Delta t}{4 \Delta x}</math>

:<math> \beta= \frac{k\Delta t}{2}</math>

把 <1>, <2>, <3>, <4>, <5>, <6>, ''α'', ''β'' 和 ''λ'' 代入 <0>. 把新时间项(''j''+1)代入到左边，当前时间项(''j'')代入到右边，将得到

:<math> -\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}</math>

第一个通道只能与下一个通道(''M'')有关系，因此表达式可以简化为：
:<math> -(\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}</math>

同样地， 最后一个通道只与前一个通道（''N''）有关联，因此表达式可以简化为
:<math> -\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}</math>

为求解此线性方程组，需要知道边界条件在通道始端就已经给定了。

<math> C_0^{j}</math>: 当前时间步某通道的初始条件 

<math> C_{0}^{j+1}</math>: 下一时间步某通道的初始条件

<math> C_{N0}^{j}</math>: 前一通道到当前时间步下某通道的初始条件

<math> C_{M0}^{j}</math>: 下一通道到当前时间步下某通道的初始条件

对于通道的末端最后一个节点，最方便的条件是是绝热近似，则

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

当且只当
:<math> C_{i + 1}^{j+1} = C_{i - 1}^{j+1} \, </math>
时，这一条件才被满足。

以3个通道，5个节点为例，可以将线性系统问题表示为

:<math> \begin{bmatrix}AA\end{bmatrix}\begin{bmatrix}C^{j+1}\end{bmatrix}=[BB][C^{j}]+[d]</math>

其中，

: <math>\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}</math>  <math>\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}</math>

需要清楚的是，''AA''和''BB''是由四个不同子矩阵组成的矩阵，

: <math>\mathbf{AA} = \begin{bmatrix}
AA1 & AA3 & 0\\
AA3 & AA2 & AA3\\
0 & AA3 & AA1\end{bmatrix}</math>

: <math>\mathbf{BB} = \begin{bmatrix}
BB1 & -AA3 & 0\\
-AA3 & BB2 & -AA3\\
0 & -AA3 & BB1\end{bmatrix}</math>

其中上述矩阵的的矩阵元对应于下一个矩阵和额外的4x4[[零矩阵]]。请注意，矩阵''AA''和''BB''的大小为12x12

: <math>\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}</math>

: <math>\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}</math>
: <math>\mathbf{AA3} = \begin{bmatrix}
-\beta & 0 & 0 & 0 \\
0 & -\beta & 0 & 0 \\
0 & 0 & -\beta & 0 \\
0 & 0 & 0 & -\beta \end{bmatrix}</math>

: <math>\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}</math> &nbsp; & &nbsp; 

: <math>\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}</math>

这里的''d''矢量用于保证边界条件成立。在此示例中为12x1的矢量。

: <math>\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}</math>

为了找到任意时间下污染物的聚集情况，需要对以下方程进行迭代计算：
:<math> \begin{bmatrix}C^{j+1}\end{bmatrix}=\begin{bmatrix}AA^{-1}\end{bmatrix}([BB][C^{j}]+[d])</math>

=== 二维扩散问题 ===
將擴散問題延伸到二維的{{link-en|笛卡爾網格|Cartesian grid}}，推導方程類似，但結果會是{{link-en|带形矩阵|Banded matrix||的方程式，不是[[三角矩陣]]，二維的熱方程

:<math>\frac{\partial u}{\partial t} = a \left(\frac{\partial^2 u}{\partial x^2} + \frac{\partial^2 u}{\partial y^2}\right)</math>

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

:<math>\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}</math>

此方程可以再重組，配合{{link-en|柯朗数|Courant number}}再進行簡化

:<math>\mu = \frac{a \Delta t}{(\Delta x)^2}.</math>

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

:<math>\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}</math>

==應用在金融數學上==
許多的現象都可以用[[熱方程]]（[[金融數學]]上稱為擴散方程）來[[數學建模|建模]]，因此克兰克－尼科尔森方法也可以用在這些領域中<ref>{{cite book | title = The Mathematics of Financial Derivatives: A Student Introduction | author1 = Wilmott, P. | author2 = Howison, S. | author3 = Dewynne, J. | publisher = Cambridge Univ. Press | year = 1995 | isbn = 0-521-49789-2 | url = http://books.google.co.in/books?hl=en&q=The%20Mathematics%20of%20Financial%20Derivatives%20Wilmott&um=1&ie=UTF-8&sa=N&tab=wp | access-date = 2015-02-16 | archive-date = 2012-10-19 | archive-url = https://web.archive.org/web/20121019072413/http://books.google.co.in/books?hl=en&q=The%20Mathematics%20of%20Financial%20Derivatives%20Wilmott&um=1&ie=UTF-8&sa=N&tab=wp | dead-url = no }}</ref>。尤其金融衍生工具定價用的[[布萊克-休斯模型]]可以轉換為熱方程，因此[[期權定價]]的[[數值方法|數值解]]可以用克兰克－尼科尔森方法求得。

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

==相關條目==
*[[金融數學]]
*[[梯形法则 (微分方程)]]

==參考資料==
{{reflist}}

{{數值偏微分方程}}

[[Category:金融数学]]
[[Category:数值微分方程]]
[[Category:有限差分]]