克兰克-尼科尔森方法(英语:Crank–Nicolson method)是一种数值分析的有限差分法,可用于数值求解热方程以及类似形式的偏微分方程[1]。它在时间方向上是隐式的二阶方法,可以写成隐式的龙格-库塔法,数值稳定。该方法诞生于20世纪,由约翰·克兰克与菲利斯·尼科尔森发展[2]。
可以证明克兰克-尼科尔森方法对于扩散方程(以及许多其他方程)是无条件稳定[3]。但是,如果时间步长Δt乘以热扩散率,再除以空间步长平方Δx2的值过大(根据冯诺依曼稳定性分析,以大于1/2为准),近似解中将存在虚假的振荡或衰减。基于这个原因,当要求大时间步或高空间分辨率的时候,往往会采用数值精确较差的后向欧拉法进行计算,这样即可以保证稳定,又避免了解的伪振荡。
克兰克-尼科尔森方法在空间域上的使用中心差分;而时间域上应用梯形公式,保证了时间域上的二阶收敛。例如,一维偏微分方程
![{\displaystyle {\frac {\partial u}{\partial t}}=F\left(u,x,t,{\frac {\partial u}{\partial x}},{\frac {\partial ^{2}u}{\partial x^{2}}}\right)}](https://wikimedia.org/api/rest_v1/media/math/render/svg/32a61d40b0f81e34bef506644f19c84ad194bfcf)
令
,则通过克兰克-尼科尔森方法导出的差分方程是第n步上采用前向欧拉方法与第n+1步上采用后向欧拉方法的平均值(注意,克兰克-尼科尔森方法本身不是这两种方法简单地取平均,方程对解隐式依赖)。
(前向欧拉方法)
(后向欧拉方法)
(克兰克-尼科尔森方法)
对于F,通过中心差分方法使其在空间上是离散的。
注意,这是一个隐式方法,需要求解代数方程组以得到时间域上的下一个u值。如果偏微分方程是非线性的,中心差分后得到的方程依旧是非线性方程系统,因此在时间步上推进会涉及求解非线性代数方程组。许多问题中,特别是线性扩散,代数方程中的矩阵是三对角的,通过三对角矩阵算法可以高效求解,这样,算法的时间复杂度由直接求解全矩阵的
转化为
。
线性扩散问题[编辑]
![{\displaystyle {\frac {\partial u}{\partial t}}=a{\frac {\partial ^{2}u}{\partial x^{2}}}}](https://wikimedia.org/api/rest_v1/media/math/render/svg/d6078eb4f0ae17369c3395b4cfe573eaf32d5fa0)
通过克兰克-尼科尔森方法将得到离散方程
![{\displaystyle {\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)}](https://wikimedia.org/api/rest_v1/media/math/render/svg/cbec5afa28c7d582f863d9c928daeec02afecf2c)
引入变量
:
![{\displaystyle -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}\,}](https://wikimedia.org/api/rest_v1/media/math/render/svg/bab440590a4a04c91b5f6285937f4583a00d66b9)
这是一个三对角问题,应用三对角矩阵算法(追赶法)即可得到
,而不需要对矩阵直接求逆。
![{\displaystyle {\frac {\partial u}{\partial t}}=a(u){\frac {\partial ^{2}u}{\partial x^{2}}}}](https://wikimedia.org/api/rest_v1/media/math/render/svg/ce58b696db770c28037513d4fa1ae5178cf94a32)
离散化后则会得到非线性方程系统。但是某些情况下,通过使用a的旧值,即用
替代
,可将问题线性化。其他时候,也可能在保证稳定性的基础上使用显式方法估计
一维多通道连接的扩散问题[编辑]
这种模型可以用于描述水流中含稳定污染流,但只有一维信息的情况。它可以简化为一维问题并得到有价值的信息。
可对水中污染溶质富集的问题进行建模,这种问题由三部分组成:已知的扩散方程(
为常量),平流分量(即由速度场导致的系统在空间上的变化,表示为常量Ux),以及与纵向通道k旁流的相互作用。
![{\displaystyle \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})}](https://wikimedia.org/api/rest_v1/media/math/render/svg/70267f956d1242abc79e64faa2d26c9df569231b)
其中C表示污染物的富集水平,下标N和M分别对应上一通道和下一通道。
克兰克-尼科尔森方法(i对应位置,j对应时间)将以上偏微分方程中的每个部分变换为
![{\displaystyle \langle 1\rangle {\frac {\partial C}{\partial t}}={\frac {C_{i}^{j+1}-C_{i}^{j}}{\Delta t}}}](https://wikimedia.org/api/rest_v1/media/math/render/svg/67397656d715a32193775427e34bc62c66f15040)
![{\displaystyle \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)}](https://wikimedia.org/api/rest_v1/media/math/render/svg/52a6266d7ffd054187c2c6c633effcff3abd24b0)
![{\displaystyle \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)}](https://wikimedia.org/api/rest_v1/media/math/render/svg/7f5049b3dc7e90707a7ca1c106d0952625818708)
![{\displaystyle \langle 4\rangle C={\frac {1}{2}}(C_{i}^{j+1}+C_{i}^{j})}](https://wikimedia.org/api/rest_v1/media/math/render/svg/e4d92fee3593a9ff97c4263ce652783107422abe)
![{\displaystyle \langle 5\rangle C_{N}={\frac {1}{2}}(C_{Ni}^{j+1}+C_{Ni}^{j})}](https://wikimedia.org/api/rest_v1/media/math/render/svg/b1d30c8d2282b0326978e4b0431f8eeeb4121fb6)
![{\displaystyle \langle 6\rangle C_{M}={\frac {1}{2}}(C_{Mi}^{j+1}+C_{Mi}^{j})}](https://wikimedia.org/api/rest_v1/media/math/render/svg/f5fd6926209312e2d8d5ae98fbf947d305784785)
现在引入以下常量用于简化计算:
![{\displaystyle \lambda ={\frac {D_{x}\Delta t}{2\Delta x^{2}}}}](https://wikimedia.org/api/rest_v1/media/math/render/svg/81d2fa40fd7488f2380ca311a9ef2630f8d62de6)
![{\displaystyle \alpha ={\frac {U_{x}\Delta t}{4\Delta x}}}](https://wikimedia.org/api/rest_v1/media/math/render/svg/4c9ec5a2753588d8106187ffc10a30dd0665e23c)
![{\displaystyle \beta ={\frac {k\Delta t}{2}}}](https://wikimedia.org/api/rest_v1/media/math/render/svg/34d3ec0b7ba2fe02cffcc51144df9515790fa67c)
把 <1>, <2>, <3>, <4>, <5>, <6>, α, β 和 λ 代入 <0>. 把新时间项(j+1)代入到左边,当前时间项(j)代入到右边,将得到
![{\displaystyle -\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}}](https://wikimedia.org/api/rest_v1/media/math/render/svg/e42b8a48bce4f2863bf5705d341656b167be8f89)
第一个通道只能与下一个通道(M)有关系,因此表达式可以简化为:
![{\displaystyle -(\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}}](https://wikimedia.org/api/rest_v1/media/math/render/svg/05204cc367a205838a93a2dfb9a8fdd89f8d6d34)
同样地, 最后一个通道只与前一个通道(N)有关联,因此表达式可以简化为
![{\displaystyle -\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}}](https://wikimedia.org/api/rest_v1/media/math/render/svg/0a18ffa14aacc3ffb3796c28093ba6a76d4d0c16)
为求解此线性方程组,需要知道边界条件在通道始端就已经给定了。
: 当前时间步某通道的初始条件
: 下一时间步某通道的初始条件
: 前一通道到当前时间步下某通道的初始条件
: 下一通道到当前时间步下某通道的初始条件
对于通道的末端最后一个节点,最方便的条件是是绝热近似,则
![{\displaystyle {\frac {\partial C}{\partial x}}_{x=z}={\frac {(C_{i+1}-C_{i-1})}{2\Delta x}}=0}](https://wikimedia.org/api/rest_v1/media/math/render/svg/4bc64870e1db2f9fbecf737f8eac1df666dcae3a)
当且只当
![{\displaystyle C_{i+1}^{j+1}=C_{i-1}^{j+1}\,}](https://wikimedia.org/api/rest_v1/media/math/render/svg/7d90098f5c25cdd8cf281cb7ce0a6d0d1271a701)
时,这一条件才被满足。
以3个通道,5个节点为例,可以将线性系统问题表示为
![{\displaystyle {\begin{bmatrix}AA\end{bmatrix}}{\begin{bmatrix}C^{j+1}\end{bmatrix}}=[BB][C^{j}]+[d]}](https://wikimedia.org/api/rest_v1/media/math/render/svg/8e2c0004cff53a3084cd60a0ff6963149780390b)
其中,
![{\displaystyle \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}}}](https://wikimedia.org/api/rest_v1/media/math/render/svg/274b425442a35a5e3df6efedd41bfd601effd696)
需要清楚的是,AA和BB是由四个不同子矩阵组成的矩阵,
![{\displaystyle \mathbf {AA} ={\begin{bmatrix}AA1&AA3&0\\AA3&AA2&AA3\\0&AA3&AA1\end{bmatrix}}}](https://wikimedia.org/api/rest_v1/media/math/render/svg/630c8627535302062431549be8c9bac1f727036c)
![{\displaystyle \mathbf {BB} ={\begin{bmatrix}BB1&-AA3&0\\-AA3&BB2&-AA3\\0&-AA3&BB1\end{bmatrix}}}](https://wikimedia.org/api/rest_v1/media/math/render/svg/47da9632045b13b1c4b894a5edb76f09299d75c8)
其中上述矩阵的的矩阵元对应于下一个矩阵和额外的4x4零矩阵。请注意,矩阵AA和BB的大小为12x12
![{\displaystyle \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}}}](https://wikimedia.org/api/rest_v1/media/math/render/svg/7151f03ea2fea8ef5546b2fbad432635b8fb49ed)
![{\displaystyle \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}}}](https://wikimedia.org/api/rest_v1/media/math/render/svg/a70e9599cd41a2226d043ed1553e95d859652f2b)
![{\displaystyle \mathbf {AA3} ={\begin{bmatrix}-\beta &0&0&0\\0&-\beta &0&0\\0&0&-\beta &0\\0&0&0&-\beta \end{bmatrix}}}](https://wikimedia.org/api/rest_v1/media/math/render/svg/4befc396376d975f6ba3deac20d3c8b2a4e3b027)
&
![{\displaystyle \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}}}](https://wikimedia.org/api/rest_v1/media/math/render/svg/aef80fe9db983b64e5b3012ead2bbd5b6a0ddb49)
这里的d矢量用于保证边界条件成立。在此示例中为12x1的矢量。
![{\displaystyle \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}}}](https://wikimedia.org/api/rest_v1/media/math/render/svg/071394782eb2eb337b52a5a95e0f054c549fa757)
为了找到任意时间下污染物的聚集情况,需要对以下方程进行迭代计算:
![{\displaystyle {\begin{bmatrix}C^{j+1}\end{bmatrix}}={\begin{bmatrix}AA^{-1}\end{bmatrix}}([BB][C^{j}]+[d])}](https://wikimedia.org/api/rest_v1/media/math/render/svg/82b9594d4023d1133ab12023a4b8a801803b0c86)
二维扩散问题[编辑]
将扩散问题延伸到二维的笛卡尔网格,推导方程类似,但结果会是{{link-en|带形矩阵|Banded matrix||的方程式,不是三角矩阵,二维的热方程
![{\displaystyle {\frac {\partial u}{\partial t}}=a\left({\frac {\partial ^{2}u}{\partial x^{2}}}+{\frac {\partial ^{2}u}{\partial y^{2}}}\right)}](https://wikimedia.org/api/rest_v1/media/math/render/svg/46663350fcf0f5ceb1615f6dd1152701f1e2e8c8)
假设网格满足
的特性,即可通过克兰克-尼科尔森方法将得到离散方程
![{\displaystyle {\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}}}](https://wikimedia.org/api/rest_v1/media/math/render/svg/5d074ece8b5df3b7f5e563a01d0fcf58768cd407)
此方程可以再重组,配合柯朗数再进行简化
![{\displaystyle \mu ={\frac {a\Delta t}{(\Delta x)^{2}}}.}](https://wikimedia.org/api/rest_v1/media/math/render/svg/8aedf0528245a91706f748bc38a000f294e31836)
在克兰克-尼科尔森方法下,不需要为了稳定性而限制柯朗数的上限,不过为了数值稳定度,柯朗数仍不能太高,可以将方程式重写如下:
![{\displaystyle {\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}}}](https://wikimedia.org/api/rest_v1/media/math/render/svg/7471ad0140bf6b69d6de2e0a1b9158a06079ba8a)
应用在金融数学上[编辑]
许多的现象都可以用热方程(金融数学上称为扩散方程)来建模,因此克兰克-尼科尔森方法也可以用在这些领域中[4]。尤其金融衍生工具定价用的布莱克-休斯模型可以转换为热方程,因此期权定价的数值解可以用克兰克-尼科尔森方法求得。
因为期权定价若超过基本假设(例如改变股息)时,无法求得解析解,需要用上述方式求得。不过若是非平滑的最后条件(大部份的金融商品都是如此),克兰克-尼科尔森方法会有数值的震荡,无法用滤波方式平缓。在期权定价上会反映在履约价Γ的变动。因此,一开始几个步骤需要用其他比较不会震荡的方法(如全隐式有限差分法)。
相关条目[编辑]
参考资料[编辑]