克兰克-尼科尔森方法

维基百科,自由的百科全书
跳转至: 导航搜索

克兰克-尼科尔森方法是一種数值分析有限差分法,可用于数值求解热方程以及类似形式的偏微分方程[1]。它在时间方向上是隐式英语Explicit and implicit methods的二阶方法,可以寫成隐式的龍格-庫塔法数值稳定。该方法诞生于20世纪,由約翰·克蘭克英语John Crank菲利斯·尼科爾森英语Phyllis Nicolson发展[2]

可以证明克兰克-尼科尔森方法对于扩散方程(以及许多其他方程)是无条件稳定[3]。但是,如果时间步长Δt乘以熱擴散率,再除以空间步长平方Δx2的值过大(根據馮諾依曼穩定性分析,以大于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} - 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表示污染物的富集水平,下标NM分别对应上一通道和下一通道。

克兰克-尼科尔森方法(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)} + 
 \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}

需要清楚的是,AABB是由四个不同子矩阵组成的矩阵,

\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零矩阵。请注意,矩阵AABB的大小为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])

參考資料[编辑]

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