克兰克-尼科尔森方法在空间域上的使用中心差分;而时间域上应用梯形公式,保证了时间域上的二阶收敛。例如,一维偏微分方程
∂ u ∂ t = F ( u , x , t , ∂ u ∂ x , ∂ 2 u ∂ x 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)} 令u ( i Δ x , n Δ t ) = u i n {\displaystyle u(i\Delta x,n\Delta t)=u_{i}^{n}\,} ,则通过克兰克-尼科尔森方法导出的差分方程是第n 步上采用前向欧拉方法与第n+1 步上采用后向欧拉方法的平均值(注意,克兰克-尼科尔森方法本身不是这两种方法简单地取平均,方程对解隐式依赖)。
u i n + 1 − u i n Δ t = F i n ( u , x , t , ∂ u ∂ x , ∂ 2 u ∂ x 2 ) {\displaystyle {\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)} (前向欧拉方法)
u i n + 1 − u i n Δ t = F i n + 1 ( u , x , t , ∂ u ∂ x , ∂ 2 u ∂ x 2 ) {\displaystyle {\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)} (后向欧拉方法)
u i n + 1 − u i n Δ t = 1 2 ( F i n + 1 ( u , x , t , ∂ u ∂ x , ∂ 2 u ∂ x 2 ) + F i n ( u , x , t , ∂ u ∂ x , ∂ 2 u ∂ x 2 ) ) {\displaystyle {\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 值。如果偏微分方程是非线性的,中心差分后得到的方程依旧是非线性方程系统,因此在时间步上推进会涉及求解非线性代数方程组。许多问题中,特别是线性扩散,代数方程中的矩阵是三对角的,通过三对角矩阵算法 可以高效求解,这样,算法的时间复杂度由直接求解全矩阵的O ( n 3 ) {\displaystyle {\mathcal {O}}(n^{3})} 转化为O ( n ) {\displaystyle {\mathcal {O}}(n)} 。
线性扩散问题
编辑
∂ u ∂ t = a ∂ 2 u ∂ x 2 {\displaystyle {\frac {\partial u}{\partial t}}=a{\frac {\partial ^{2}u}{\partial x^{2}}}} 通过克兰克-尼科尔森方法将得到离散方程
u i n + 1 − u i n Δ t = a 2 ( Δ x ) 2 ( ( 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 ) ) {\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)} 引入变量r = a Δ t 2 ( Δ x ) 2 {\displaystyle 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 {\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}\,} 这是一个三对角问题,应用三对角矩阵算法(追赶法)即可得到u i n + 1 {\displaystyle u_{i}^{n+1}} ,而不需要对矩阵直接求逆。
∂ u ∂ t = a ( u ) ∂ 2 u ∂ x 2 {\displaystyle {\frac {\partial u}{\partial t}}=a(u){\frac {\partial ^{2}u}{\partial x^{2}}}} 离散化后则会得到非线性方程系统。但是某些情况下,通过使用a 的旧值,即用a i n ( u ) {\displaystyle a_{i}^{n}(u)\,} 替代a i n + 1 ( u ) {\displaystyle a_{i}^{n+1}(u)\,} ,可将问题线性化。其他时候,也可能在保证稳定性的基础上使用显式方法估计a i n + 1 ( u ) {\displaystyle a_{i}^{n+1}(u)\,}
一维多通道连接的扩散问题
编辑
这种模型可以用于描述水流中含稳定污染流,但只有一维信息的情况。它可以简化为一维问题并得到有价值的信息。
可对水中污染溶质富集的问题进行建模,这种问题由三部分组成:已知的扩散方程(D x {\displaystyle D_{x}} 为常量),平流分量(即由速度场导致的系统在空间上的变化,表示为常量Ux ),以及与纵向通道k旁流的相互作用。
⟨ 0 ⟩ ∂ C ∂ t = D x ∂ 2 C ∂ x 2 − U x ∂ C ∂ x − k ( C − C N ) − k ( C − C M ) {\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})} 其中C 表示污染物的富集水平,下标N 和M 分别对应上一通道和下一通道。
克兰克-尼科尔森方法(i对应位置,j对应时间)将以上偏微分方程中的每个部分变换为
⟨ 1 ⟩ ∂ C ∂ t = C i j + 1 − C i j Δ t {\displaystyle \langle 1\rangle {\frac {\partial C}{\partial t}}={\frac {C_{i}^{j+1}-C_{i}^{j}}{\Delta t}}} ⟨ 2 ⟩ ∂ 2 C ∂ x 2 = 1 2 ( Δ x ) 2 ( ( 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 ) ) {\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)} ⟨ 3 ⟩ ∂ C ∂ x = 1 2 ( ( C i + 1 j + 1 − C i − 1 j + 1 ) 2 ( Δ x ) + ( C i + 1 j − C i − 1 j ) 2 ( Δ x ) ) {\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)} ⟨ 4 ⟩ C = 1 2 ( C i j + 1 + C i j ) {\displaystyle \langle 4\rangle C={\frac {1}{2}}(C_{i}^{j+1}+C_{i}^{j})} ⟨ 5 ⟩ C N = 1 2 ( C N i j + 1 + C N i j ) {\displaystyle \langle 5\rangle C_{N}={\frac {1}{2}}(C_{Ni}^{j+1}+C_{Ni}^{j})} ⟨ 6 ⟩ C M = 1 2 ( C M i j + 1 + C M i j ) {\displaystyle \langle 6\rangle C_{M}={\frac {1}{2}}(C_{Mi}^{j+1}+C_{Mi}^{j})} 现在引入以下常量用于简化计算:
λ = D x Δ t 2 Δ x 2 {\displaystyle \lambda ={\frac {D_{x}\Delta t}{2\Delta x^{2}}}} α = U x Δ t 4 Δ x {\displaystyle \alpha ={\frac {U_{x}\Delta t}{4\Delta x}}} β = k Δ t 2 {\displaystyle \beta ={\frac {k\Delta t}{2}}} 把 <1>, <2>, <3>, <4>, <5>, <6>, α , β 和 λ 代入 <0>. 把新时间项(j +1)代入到左边,当前时间项(j )代入到右边,将得到
− β C N i j + 1 − ( λ + α ) C i − 1 j + 1 + ( 1 + 2 λ + 2 β ) C i j + 1 − ( λ − α ) C i + 1 j + 1 − β C M i j + 1 = β C N i j + ( λ + α ) C i − 1 j + ( 1 − 2 λ − 2 β ) C i j + ( λ − α ) C i + 1 j + β C M i 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}} 第一个通道只能与下一个通道(M )有关系,因此表达式可以简化为:
− ( λ + α ) C i − 1 j + 1 + ( 1 + 2 λ + β ) C i j + 1 − ( λ − α ) C i + 1 j + 1 − β C M i j + 1 = + ( λ + α ) C i − 1 j + ( 1 − 2 λ − β ) C i j + ( λ − α ) C i + 1 j + β C M i j {\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}} 同样地, 最后一个通道只与前一个通道(N )有关联,因此表达式可以简化为
− β C N i j + 1 − ( λ + α ) C i − 1 j + 1 + ( 1 + 2 λ + β ) C i j + 1 − ( λ − α ) C i + 1 j + 1 = β C N i j + ( λ + α ) C i − 1 j + ( 1 − 2 λ − β ) C i j + ( λ − α ) C i + 1 j {\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}} 为求解此线性方程组,需要知道边界条件在通道始端就已经给定了。
C 0 j {\displaystyle C_{0}^{j}} : 当前时间步某通道的初始条件
C 0 j + 1 {\displaystyle C_{0}^{j+1}} : 下一时间步某通道的初始条件
C N 0 j {\displaystyle C_{N0}^{j}} : 前一通道到当前时间步下某通道的初始条件
C M 0 j {\displaystyle C_{M0}^{j}} : 下一通道到当前时间步下某通道的初始条件
对于通道的末端最后一个节点,最方便的条件是是绝热近似,则
∂ C ∂ x x = z = ( C i + 1 − C i − 1 ) 2 Δ x = 0 {\displaystyle {\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 {\displaystyle C_{i+1}^{j+1}=C_{i-1}^{j+1}\,} 时,这一条件才被满足。
以3个通道,5个节点为例,可以将线性系统问题表示为
[ A A ] [ C j + 1 ] = [ B B ] [ C j ] + [ d ] {\displaystyle {\begin{bmatrix}AA\end{bmatrix}}{\begin{bmatrix}C^{j+1}\end{bmatrix}}=[BB][C^{j}]+[d]} 其中,
C j + 1 = [ 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 ] {\displaystyle \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}}} C j = [ 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 ] {\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}}} 需要清楚的是,AA 和BB 是由四个不同子矩阵组成的矩阵,
A A = [ A A 1 A A 3 0 A A 3 A A 2 A A 3 0 A A 3 A A 1 ] {\displaystyle \mathbf {AA} ={\begin{bmatrix}AA1&AA3&0\\AA3&AA2&AA3\\0&AA3&AA1\end{bmatrix}}} B B = [ B B 1 − A A 3 0 − A A 3 B B 2 − A A 3 0 − A A 3 B B 1 ] {\displaystyle \mathbf {BB} ={\begin{bmatrix}BB1&-AA3&0\\-AA3&BB2&-AA3\\0&-AA3&BB1\end{bmatrix}}} 其中上述矩阵的的矩阵元对应于下一个矩阵和额外的4x4零矩阵 。请注意,矩阵AA 和BB 的大小为12x12
A A 1 = [ ( 1 + 2 λ + β ) − ( λ − α ) 0 0 − ( λ + α ) ( 1 + 2 λ + β ) − ( λ − α ) 0 0 − ( λ + α ) ( 1 + 2 λ + β ) − ( λ − α ) 0 0 − 2 λ ( 1 + 2 λ + β ) ] {\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}}} A A 2 = [ ( 1 + 2 λ + 2 β ) − ( λ − α ) 0 0 − ( λ + α ) ( 1 + 2 λ + 2 β ) − ( λ − α ) 0 0 − ( λ + α ) ( 1 + 2 λ + 2 β ) − ( λ − α ) 0 0 − 2 λ ( 1 + 2 λ + 2 β ) ] {\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}}}
A A 3 = [ − β 0 0 0 0 − β 0 0 0 0 − β 0 0 0 0 − β ] {\displaystyle \mathbf {AA3} ={\begin{bmatrix}-\beta &0&0&0\\0&-\beta &0&0\\0&0&-\beta &0\\0&0&0&-\beta \end{bmatrix}}} B B 1 = [ ( 1 − 2 λ − β ) ( λ − α ) 0 0 ( λ + α ) ( 1 − 2 λ − β ) ( λ − α ) 0 0 ( λ + α ) ( 1 − 2 λ − β ) ( λ − α ) 0 0 2 λ ( 1 − 2 λ − β ) ] {\displaystyle \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}}} & B B 2 = [ ( 1 − 2 λ − 2 β ) ( λ − α ) 0 0 ( λ + α ) ( 1 − 2 λ − 2 β ) ( λ − α ) 0 0 ( λ + α ) ( 1 − 2 λ − 2 β ) ( λ − α ) 0 0 2 λ ( 1 − 2 λ − 2 β ) ] {\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}}} 这里的d 矢量用于保证边界条件成立。在此示例中为12x1的矢量。
d = [ ( λ + α ) ( C 10 j + 1 + C 10 j ) 0 0 0 ( λ + α ) ( C 20 j + 1 + C 20 j ) 0 0 0 ( λ + α ) ( C 30 j + 1 + C 30 j ) 0 0 0 ] {\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}}} 为了找到任意时间下污染物的聚集情况,需要对以下方程进行迭代计算:
[ C j + 1 ] = [ A A − 1 ] ( [ B B ] [ C j ] + [ d ] ) {\displaystyle {\begin{bmatrix}C^{j+1}\end{bmatrix}}={\begin{bmatrix}AA^{-1}\end{bmatrix}}([BB][C^{j}]+[d])} 二维扩散问题
编辑
將擴散問題延伸到二維的笛卡爾網格 ,推導方程類似,但結果會是{{link-en|带形矩阵|Banded matrix||的方程式,不是三角矩陣 ,二維的熱方程
∂ u ∂ t = a ( ∂ 2 u ∂ x 2 + ∂ 2 u ∂ y 2 ) {\displaystyle {\frac {\partial u}{\partial t}}=a\left({\frac {\partial ^{2}u}{\partial x^{2}}}+{\frac {\partial ^{2}u}{\partial y^{2}}}\right)} 假設網格滿足Δ x = Δ y {\displaystyle \Delta x=\Delta y} 的特性,即可通過克蘭克-尼科爾森方法將得到離散方程
u i , j n + 1 = u i , j n + 1 2 a Δ t ( Δ 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 − 4 u i , j n + 1 ) + ( u i + 1 , j n + u i − 1 , j n + u i , j + 1 n + u i , j − 1 n − 4 u i , j n ) ] {\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}}} 此方程可以再重組,配合柯朗数 再進行簡化
μ = a Δ t ( Δ x ) 2 . {\displaystyle \mu ={\frac {a\Delta t}{(\Delta x)^{2}}}.} 在克蘭克-尼科爾森方法下,不需要為了穩定性而限制柯朗数的上限,不過為了數值穩定度,柯朗数仍不能太高,可以將方程式重寫如下:
( 1 + 2 μ ) u i , j n + 1 − μ 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 μ ) u i , j n + μ 2 ( u i + 1 , j n + u i − 1 , j n + u i , j + 1 n + u i , j − 1 n ) . {\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}}} 應用在金融數學上
编辑
許多的現象都可以用熱方程 (金融數學 上稱為擴散方程)來建模 ,因此克兰克-尼科尔森方法也可以用在這些領域中[4] 。尤其金融衍生工具定價用的布萊克-休斯模型 可以轉換為熱方程,因此期權定價 的數值解 可以用克兰克-尼科尔森方法求得。
因為期權定價若超過基本假設(例如改變股息)時,無法求得解析解,需要用上述方式求得。不過若是非平滑的最後條件(大部份的金融商品 都是如此),克兰克-尼科尔森方法會有數值的震盪,無法用濾波方式平緩。在期權定價 上會反映在履約價 Γ的變動。因此,一開始幾個步驟需要用其他比較不會震盪的方法(如全隱式有限差分法)。
相關條目
编辑
參考資料
编辑