クランク・ニコルソン法

数値解析において、クランク・ニコルソン法は、熱方程式や類似の偏微分方程式を数値的に解くために使用される有限差分法である。^{[ 1 ]}これは時間に関して2次法である。時間に関して陰解法であり、陰解法ルンゲ・クッタ法として記述でき、数値的に安定している。この方法は、1940年代にジョン・クランクとフィリス・ニコルソンによって開発された。 ^{[ 2 ]}

拡散方程式（および他の多くの方程式）については、クランク・ニコルソン法が無条件に安定であることが示されています。^{[ 3 ]}しかし、時間ステップの比率によっては、近似解に（減衰する）偽の振動が含まれる可能性があります。${\displaystyle \Delta t}$空間ステップの2乗に対する熱拡散率の倍数、$\displaystyle \Delta x^{2}}$は大きい（フォン・ノイマン安定性解析によれば、通常1/2より大きい）。このため、大きな時間ステップや高い空間分解能が必要な場合は、安定性と振動耐性の両方を備えた、精度は劣るものの後退オイラー法がよく用いられる。

クランク・ニコルソン法は台形公式に基づいており、時間に関して2次収束性を持つ。線形方程式の場合、台形公式は陰的中点法（ガウス・ルジャンドル陰的ルンゲ・クッタ法の最も単純な例）と等価であり、幾何学的積分器としての性質も持つ。例えば、1次元の場合、偏微分方程式が次の ようになるとしよう。$${\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\Delta x,n\Delta t)=u_{i}^{n}}$そして${\displaystyle F_{i}^{n}=F}$評価対象${\displaystyle i,n}$そして$u_i^n$クランク・ニコルソン法の式は、前進オイラー法と${\displaystyle n}$そして後退オイラー法では${\displaystyle n+1}$（ただし、この方法は単純にこれら2つの方法の平均ではなく、後退オイラー方程式は解に暗黙的に依存していることに注意が必要です。）

これは暗黙的なメソッドであることに注意してください。次の値を取得するには、${\displaystyle u}$時間経過に伴い、代数方程式系を解く必要があります。偏微分方程式が非線形の場合、離散化も非線形になるため、時間を進めるには非線形代数方程式系の解法が必要となりますが、線形化は可能です。多くの問題、特に線形拡散においては、代数問題は三重対角行列となり、三重対角行列アルゴリズムを用いて効率的に解くことができます。このアルゴリズムは高速な解法を提供します。${\displaystyle {\mathcal {O}}(N)}$通常の解決策とは対照的に、直接的な解決策${\displaystyle {\mathcal {O}}(N^{3})}$完全な行列の場合、${\displaystyle N}$行列のサイズを示します。

クランク・ニコルソン法は拡散問題によく適用される。例えば、線形拡散の場合、

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

右辺に有限差分空間離散化を適用すると、クランク・ニコルソン離散化は次のようになる。

${\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)}$

または、${\displaystyle r={\frac {a\Delta t}{2(\Delta x)^{2}}}}$、

${\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 {\frac {\partial u}{\partial t}}=a(u){\frac {\partial ^{2}u}{\partial x^{2}}},}$

これは非線形代数方程式系につながり、上記のように簡単に解くことはできません。ただし、場合によっては、古い値を使用して問題を線形化することが可能になります。${\displaystyle a}$つまり、${\displaystyle a_{i}^{n}(u)}$の代わりに${\displaystyle a_{i}^{n+1}(u)}$場合によっては、推定できるかもしれません${\displaystyle a_{i}^{n+1}(u)}$明示的な方法を用いて安定性を維持する。

これは、定常流条件下にある河川や小川で汚染問題が発生した場合に、様々な目的で一般的に用いられる解決策ですが、得られる情報は一次元のみです。多くの場合、この問題は一次元の問題に単純化しても、有用な情報を得ることができます。

ここでは、水中の溶質汚染物質の濃度をモデル化します。この問題は、次の 3 つの部分から構成されています。既知の拡散方程式 (${\displaystyle D_{x}}$（定数として選択）、移流成分（速度場によってシステムが空間的に進化していることを意味する）、定数として選択${\displaystyle U_{x}}$、縦方向チャネル間の横方向の相互作用（${\displaystyle k}$):

どこ${\displaystyle C}$は汚染物質の濃度であり、添え字は${\displaystyle N}$そして${\displaystyle M}$前後のチャンネルに対応します。

クランク・ニコルソン法（${\displaystyle i}$位置を表し、${\displaystyle j}$時間）は、偏微分方程式の各要素を次のように変換します。

代数計算を簡略化するために、以下の定数を作成します。

${\displaystyle \lambda ={\frac {D_{x}\,\Delta t}{2\,\Delta x^{2}}},}$

${\displaystyle \alpha ={\frac {U_{x}\,\Delta t}{4\,\Delta x}},}$

${\displaystyle \beta ={\frac {k\,\Delta t}{2}},}$

そして、( 2 )、( 3 )、( 4 )、( 5 )、( 6 )、( 7 )を代入する。${\displaystyle \alpha }$、${\displaystyle \beta }$そして${\displaystyle \lambda }$（１ ）に代入する。次に、新しい時間項を左辺に置く（${\displaystyle j+1}$）と右側の現在の時間項（${\displaystyle 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}={}}$

${\displaystyle \qquad \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}.}$

最初のチャネルをモデル化するために、次のチャネルとのみ接触できることに気づきます（${\displaystyle 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}={}}$

${\displaystyle \qquad {}+(\lambda +\alpha )C_{i-1}^{j}+(1-2\lambda -\beta )C_{i}^{j}+(\lambda -\alpha )C_{i+1}^{j}+\beta C_{Mi}^{j}.}$

同様に、最後のチャネルをモデル化するために、前のチャネルとのみ接触できることに気づきます（${\displaystyle 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}={}}$

${\displaystyle \qquad \beta C_{Ni}^{j}+(\lambda +\alpha )C_{i-1}^{j}+(1-2\lambda -\beta )C_{i}^{j}+(\lambda -\alpha )C_{i+1}^{j}.}$

この線形方程式系を解くには、まずチャネルの始点に境界条件を与える必要があることに注意しなければなりません。

${\displaystyle C_{0}^{j}}$：現在の時間ステップにおけるチャネルの境界条件、

${\displaystyle C_{0}^{j+1}}$: 次のタイムステップにおけるチャネルの境界条件、

${\displaystyle C_{N0}^{j}}$: 現在の時間ステップで解析されるチャネルの前のチャネルの境界条件、

${\displaystyle C_{M0}^{j}}$: 現在の時間ステップで解析されているチャネルの次のチャネルに対する境界条件。

チャネルの最後のセルについては（${\displaystyle z}$)最も都合の良い条件は断熱条件となるので、

${\displaystyle \left.{\frac {\partial C}{\partial x}}\right|_{x=z}={\frac {C_{i+1}-C_{i-1}}{2\,\Delta x}}=0.}$

この条件は、（null値に関係なく）

${\displaystyle C_{i+1}^{j+1}=C_{i-1}^{j+1}.}$

3つのチャネルと5つのノード（初期境界条件を含む）の場合について、この問題を（行列形式で）解いてみましょう。これを線形システム問題として表現します。

${\displaystyle \mathbf {AA} \,\mathbf {C^{j+1}} =\mathbf {BB} \,\mathbf {C^{j}} +\mathbf {d} ,}$

どこ

${\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}},\quad \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は4つの異なるサブアレイで構成されるアレイである必要があることを理解しなければなりません（この例では3つのチャネルのみを考慮していますが、上記で説明した主要な部分を網羅しています）。

${\displaystyle \mathbf {AA} ={\begin{bmatrix}AA1&AA3&0\\AA3&AA2&AA3\\0&AA3&AA1\end{bmatrix}},\quad \mathbf {BB} ={\begin{bmatrix}BB1&-AA3&0\\-AA3&BB2&-AA3\\0&-AA3&BB1\end{bmatrix}},}$

上記の要素は次の配列に対応し、さらに4×4のゼロで埋められた領域が追加されます。AAとBBのサイズは12×12であることに注意してください。

${\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}},}$

${\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}},}$

${\displaystyle \mathbf {AA3} ={\begin{bmatrix}-\beta &0&0&0\\0&-\beta &0&0\\0&0&-\beta &0\\0&0&0&-\beta \end{bmatrix}},}$

${\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}},}$

${\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ベクトルは境界条件を保持するために使用されます。この例では、12×1のベクトルです。

${\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}}.}$

任意の時点での濃度を求めるには、次の式を繰り返し計算する必要があります。

${\displaystyle \mathbf {C^{j+1}} =\mathbf {AA} ^{-1}(\mathbf {BB} \,\mathbf {C^{j}} +\mathbf {d} ).}$

均一な直交座標グリッド上で2次元に拡張する場合、導出は同様であり、結果は三重対角方程式ではなく帯状対角方程式のシステムにつながる可能性がある。2次元熱方程式

${\displaystyle {\frac {\partial u}{\partial t}}=a\,\nabla ^{2}u,}$

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

クランク・ニコルソン離散化を用いて解くことができる

${\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})\\&+(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}}}$

正方形のグリッドが使用されていると仮定すると、${\displaystyle \Delta x=\Delta y}$この方程式は、項を並べ替え、 CFL数を用いることでいくらか簡略化できる。

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

Crank–Nicolson数値スキームでは、安定性のために低いCFL数は必要ないが、数値精度のためには必要である。このスキームは次のように記述できる。

${\displaystyle (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)}$

${\displaystyle \qquad =(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).}$

このような線形システムを解くにはコストがかかります。そのため、数値偏微分方程式を解くために交互方向陰解法を実装できます。この方法では、割り当てられた時間ステップの半分で一方の次元を陰的に、もう一方の次元を陽的に扱い、残りの半分ではその逆を行います。この戦略の利点は、陰解法ソルバーが解く必要があるのは三重対角行列アルゴリズムだけであることです。真のクランク・ニコルソン解とADI近似解の差は、精度がオーダーで次のようになります。${\displaystyle O(\Delta t^{4})}$したがって、十分小さな時間ステップであれば無視できる。^{[ 4 ]}

クランク・ニコルソン法は陰解法であるため、一般的に正確な解を求めることは不可能です。そのため、反復法を用いて解を収束させる必要があります。一つの方法として、ニュートン法を用いて予測値を収束させる方法がありますが、これにはヤコビ行列の計算が必要です。計算流体力学や数値相対性理論のような高次元システムの場合、このヤコビ行列を計算することは現実的ではないかもしれません。

ヤコビアンを用いない代替手段として、固定点反復法がある。${\displaystyle f}$はシステムの速度であり、その場合、クランク・ニコルソン予測はマップの固定点となる。${\displaystyle \Phi (x)=x_{0}+{\frac {h}{2}}\left[f(x_{0})+f(x)\right].}$マップの反復処理${\displaystyle x^{(i+1)}=\Phi (x^{(i)})}$収束しない、パラメータ化されたマップ${\displaystyle \Theta (x,\alpha )=\alpha x+(1-\alpha )\Phi (x)}$、 と${\displaystyle \alpha \in (0,1)}$の方が挙動が良いかもしれません。展開した形式では、更新式は次のようになります。

${\displaystyle x^{i+1}=\alpha x^{i}+(1-\alpha )\left[x_{0}+{\frac {h}{2}}\left(f(x_{0})+f(x^{i})\right)\right],}$

どこ${\displaystyle x^{i}}$現在の推測と${\displaystyle x_{i-1}}$これは前のタイムステップです。

高次元システムであっても、このマップの反復処理は驚くほど速く収束する可能性がある。

熱方程式（金融数学では拡散方程式と呼ばれることが多い）でモデル化できる他の多くの現象があるため、クランク・ニコルソン法はそれらの分野にも適用されています。^{[ 5 ]}特に、ブラック・ショールズオプション価格モデルの微分方程式は熱方程式に変換できるため、 クランク・ニコルソン法でオプション価格の数値解を得ることができます。

金融におけるこの方法の重要性は、オプション価格決定問題が標準的な仮定を超えて拡張された場合（例えば、配当金の変動を組み込む場合）、閉形式で解くことができないが、この方法を用いることで解くことができる点にある。ただし、ほとんどの金融商品で発生するような、最終条件が滑らかでない場合、数値振動が抑制されないため、クランク・ニコルソン法は適切ではないことに注意が必要である。バニラオプションの場合、これは行使価格付近のガンマ値の振動につながる。したがって、特別な減衰初期化手順（例えば、完全陰解法による有限差分法）が必要となる。

• 金融数学
• 台形公式

• 科学者と技術者のための数値偏微分方程式の手法、オープンアクセス形式の講義と数値偏微分方程式のコード
• 移流方程式に対するクランク・ニコルソン法の適用および実装方法の例