克兰克-尼科尔森方法

克兰克－尼科尔森方法（Template:Lang-en）是一種数值分析的有限差分法，可用于数值求解热方程以及类似形式的偏微分方程。它在时间方向上是隐式的二阶方法，可以寫成隐式的龍格－庫塔法，数值稳定。该方法诞生于20世纪，由約翰·克蘭克与菲利斯·尼科爾森发展。

可以证明克兰克－尼科尔森方法对于扩散方程（以及许多其他方程）是无条件稳定。但是，如果时间步长ΔTemplate:Var乘以熱擴散率，再除以空间步长平方ΔTemplate:VarTemplate:Sup的值过大（根據馮諾依曼穩定性分析，以大于1/2為準），近似解中将存在虚假的振荡或衰减。基于这个原因，当要求大时间步或高空间分辨率的时候，往往会采用数值精确较差的后向欧拉法进行计算，这样即可以保证稳定，又避免了解的伪振荡。

方法

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

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

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

$\frac{u_i^{n+1}-u_i^n}{\mathrm{\Delta }t}=F_i^n(u,x,t,\frac{\partial u}{\partial x},\frac{{\partial }^{2}u}{\partial x^2})$ (前向欧拉方法)

$\frac{u_i^{n+1}-u_i^n}{\mathrm{\Delta }t}=F_i^{n+1}(u,x,t,\frac{\partial u}{\partial x},\frac{{\partial }^{2}u}{\partial x^2})$ (后向欧拉方法)

$\frac{u_i^{n+1}-u_i^n}{\mathrm{\Delta }t}=\frac{1}{2}(F_i^{n+1}(u,x,t,\frac{\partial u}{\partial x},\frac{{\partial }^{2}u}{\partial x^2})+F_i^n(u,x,t,\frac{\partial u}{\partial x},\frac{{\partial }^{2}u}{\partial x^2}))$ (克兰克－尼科尔森方法)

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

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

示例

线性扩散问题

• 线性扩散方程

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

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

$\frac{u_i^{n+1}-u_i^n}{\mathrm{\Delta }t}=\frac{a}{2(\mathrm{\Delta }x{)}^2}((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))$

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

$-r{u}_{i+1}^{n+1}+(1+2r)u_i^{n+1}-r{u}_{i-1}^{n+1}=r{u}_{i+1}^n+(1-2r)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旁流的相互作用。

$⟨0⟩\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对应时间）将以上偏微分方程中的每个部分变换为

$⟨1⟩\frac{\partial C}{\partial t}=\frac{C_i^{j+1}-C_i^j}{\mathrm{\Delta }t}$

$⟨2⟩\frac{{\partial }^{2}C}{\partial x^2}=\frac{1}{2(\mathrm{\Delta }x{)}^2}((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))$

$⟨3⟩\frac{\partial C}{\partial x}=\frac{1}{2}(\frac{(C_{i+1}^{j+1}-C_{i-1}^{j+1})}{2(\mathrm{\Delta }x)}+\frac{(C_{i+1}^j-C_{i-1}^j)}{2(\mathrm{\Delta }x)})$

$⟨4⟩C=\frac{1}{2}(C_i^{j+1}+C_i^j)$

$⟨5⟩C_N=\frac{1}{2}(C_{Ni}^{j+1}+C_{Ni}^j)$

$⟨6⟩C_M=\frac{1}{2}(C_{Mi}^{j+1}+C_{Mi}^j)$

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

$\lambda =\frac{D_x\mathrm{\Delta }t}{2\mathrm{\Delta }{x}^2}$

$\alpha =\frac{U_x\mathrm{\Delta }t}{4\mathrm{\Delta }x}$

$\beta =\frac{k\mathrm{\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\mathrm{\Delta }x}=0$

当且只当

$C_{i+1}^{j+1}=C_{i-1}^{j+1}$

时，这一条件才被满足。

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

$\left[\begin{array}{c}AA\end{array}\right]\left[\begin{array}{c}{C}^{j+1}\end{array}\right]=[BB][C^j]+[d]$

其中，

${𝐂}^{𝐣\mathbf{+}\mathrm{𝟏}}=\left[\begin{array}{c}{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{array}\right]$ ${𝐂}^{𝐣}=\left[\begin{array}{c}{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{array}\right]$

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

$𝐀𝐀=\left[\begin{array}{ccc}AA1& AA3& 0\\ AA3& AA2& AA3\\ 0& AA3& AA1\end{array}\right]$

$𝐁𝐁=\left[\begin{array}{ccc}BB1& -AA3& 0\\ -AA3& BB2& -AA3\\ 0& -AA3& BB1\end{array}\right]$

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

$𝐀𝐀\mathrm{𝟏}=\left[\begin{array}{cccc}(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{array}\right]$

$𝐀𝐀\mathrm{𝟐}=\left[\begin{array}{cccc}(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{array}\right]$

$𝐀𝐀\mathrm{𝟑}=\left[\begin{array}{cccc}-\beta & 0& 0& 0\\ 0& -\beta & 0& 0\\ 0& 0& -\beta & 0\\ 0& 0& 0& -\beta \end{array}\right]$

$𝐁𝐁\mathrm{𝟏}=\left[\begin{array}{cccc}(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{array}\right]$   &  

$𝐁𝐁\mathrm{𝟐}=\left[\begin{array}{cccc}(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{array}\right]$

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

$𝐝=\left[\begin{array}{c}(\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{array}\right]$

为了找到任意时间下污染物的聚集情况，需要对以下方程进行迭代计算：

$\left[\begin{array}{c}{C}^{j+1}\end{array}\right]=\left[\begin{array}{c}A{A}^{-1}\end{array}\right]([BB][C^j]+[d])$

二维扩散问题

將擴散問題延伸到二維的Template:Link-en，推導方程類似，但結果會是{{link-en|带形矩阵|Banded matrix||的方程式，不是三角矩陣，二維的熱方程

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

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

$\begin{array}{cc}{u}_{i,j}^{n+1}& =u_{i,j}^n+\frac{1}{2}\frac{a\mathrm{\Delta }t}{(\mathrm{\Delta }x{)}^2}[(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})\\ & +(u_{i+1,j}^n+u_{i-1,j}^n+u_{i,j+1}^n+u_{i,j-1}^n-4u_{i,j}^n)]\end{array}$

此方程可以再重組，配合Template:Link-en再進行簡化

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

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

$\begin{array}{cc}& (1+2\mu )u_{i,j}^{n+1}-\frac{\mu }{2}(u_{i+1,j}^{n+1}+u_{i-1,j}^{n+1}+u_{i,j+1}^{n+1}+u_{i,j-1}^{n+1})\\ & =(1-2\mu )u_{i,j}^n+\frac{\mu }{2}(u_{i+1,j}^n+u_{i-1,j}^n+u_{i,j+1}^n+u_{i,j-1}^n).\end{array}$

應用在金融數學上

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

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

相關條目

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