克兰克－尼科尔森方法

克兰克－尼科尔森方法（英语：Crank–Nicolson method）是一种数值分析的有限差分法，可用于数值求解热方程以及类似形式的偏微分方程^[1]。它在时间方向上是隐式的二阶方法，可以写成隐式的龙格－库塔法，数值稳定。该方法诞生于20世纪，由约翰·克兰克与菲利斯·尼科尔森发展^[2]。

可以证明克兰克－尼科尔森方法对于扩散方程（以及许多其他方程）是无条件稳定^[3]。但是，如果时间步长Δt乘以热扩散率，再除以空间步长平方Δx²的值过大（根据冯诺依曼稳定性分析，以大于1/2为准），近似解中将存在虚假的振荡或衰减。基于这个原因，当要求大时间步或高空间分辨率的时候，往往会采用数值精确较差的后向欧拉法（英语：Backward Euler method）进行计算，这样即可以保证稳定，又避免了解的伪振荡。

方法

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

{\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])

二维扩散问题

将扩散问题延伸到二维的笛卡尔网格（英语：Cartesian grid），推导方程类似，但结果会是{{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}}

此方程可以再重组，配合柯朗数（英语：Courant number）再进行简化

\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}}

应用在金融数学上

许多的现象都可以用热方程（金融数学上称为扩散方程）来建模，因此克兰克－尼科尔森方法也可以用在这些领域中^[4]。尤其金融衍生工具定价用的布莱克-休斯模型可以转换为热方程，因此期权定价的数值解可以用克兰克－尼科尔森方法求得。

因为期权定价若超过基本假设（例如改变股息）时，无法求得解析解，需要用上述方式求得。不过若是非平滑的最后条件（大部分的金融商品都是如此），克兰克－尼科尔森方法会有数值的震荡，无法用滤波方式平缓。在期权定价上会反映在履约价Γ的变动。因此，一开始几个步骤需要用其他比较不会震荡的方法（如全隐式有限差分法）。

参考资料

^ Tuncer Cebeci. Convective Heat Transfer. Springer. 2002. ISBN 0-9668461-4-1.
^ 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.
^ 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 $u_{t}=au_{xx}$ .
^ 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）.

[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 $u_{t}=au_{xx}$ .

[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）.

[1]

[2]

[3]

[4]