第29話  Von Neumann の安定性解析

この記事を読む助けになる数学の記事
・数値計算：1話、 2話、 4話、 5話、 6話、 7話、 22話、 23話、 24話、 26話、 27話

　微分方程式を差分近似により表した差分方程式を解くことにより数値解を求める数値積分を総称して差分法と呼びます。 差分法では実数を有限な桁数の有理数で近似して取り扱うため、四捨五入による丸め誤差が生じます。
　時刻 \( t \) を変数に持つ偏微分方程式の数値解を、ある差分法の手法により求める際に、 時間発展につれて丸め誤差が限りなく増大する場合、その手法は不安定と呼ばれます。 Von Neumann の安定性解析は、定数係数線形偏微分方程式を差分法により数値的に解く際に、差分法の手法が不安定となる条件を丸め誤差のフーリエ級数を用いて与える手法です。
　以下では、1次元移流方程式と2次元拡散方程式を陽的差分法で解く場合を例にして Von Neumann の安定性解析を解説していきます。

　1次元移流方程式は次の偏微分方程式です。 \[ \begin{align} \frac{\partial u}{\partial t} = -c \frac{\partial u}{\partial x} \ \ \ \left( c \neq 0 \right) \end{align}\] この偏微分方程式を時間の刻み幅を \( k \) 、空間の刻み幅を \( h \) として前進差分近似の差分方程式の形で書き下すと次式になります。 \[ \begin{align} u \left( x, t + k \right) = u \left( x, t \right) - \frac{ck}{h} \left\{ u \left( x + h,t \right) - u \left( x,t \right) \right\} \ \ \ \ldots (1) \end{align}\] ここで、丸め誤差を \( \tilde{u} \left( x, t \right) \) と表すと、これは \( (1) \) 式と同じ差分方程式に従って変化します。すなわち、次式が成り立ちます。 \[ \begin{align} \tilde{u} \left( x, t + k \right) = \tilde{u} \left( x, t \right) - \frac{ck}{h} \left\{ \tilde{u} \left( x + h,t \right) - \tilde{u} \left( x,t \right) \right\} \ \ \ \ldots (2) \end{align}\] ところで、時刻 \( t \) における丸め誤差はフーリエ級数を用いて表すことができます。 \[ \begin{align} \tilde{u} \left( x, t \right) = \sum ^{\infty} _{n = - \infty} c_n \left( t \right) e^{inpx} \ \ \ \ldots (3) \end{align}\] ここで、\( p \) は計算対象としている空間領域に応じた任意の実数です。\( (3) \) 式において、フーリエ係数 \( c_n \) を時刻 \( t \) の関数としていることに注意してください。 \( (3) \) 式を \( (2) \) 式に代入すると、次式になります。 \[ \begin{align} & \ \ \ \ \ \ \sum ^{\infty} _{n = - \infty} c_n \left( t + k \right) e^{inpx} \\\\ &= \sum ^{\infty} _{n = - \infty} c_n \left( t \right) e^{inpx} - \frac{ck}{h} \left\{ \sum ^{\infty} _{n = - \infty} c_n \left( t \right) e^{inp \left( x + h \right)} - \sum ^{\infty} _{n = - \infty} c_n \left( t \right) e^{inpx} \right\} \ \ \ \ldots (4) \end{align}\] \( (2) \) 式が線形差分方程式であることから、\( (4) \) 式を調べるにあたって 、ある整数 \( N \) についてのフーリエ係数 \( c_N \) が時間発展につれて限りなく大きくなるか否かを個別に調べても一般性を失いません。 そこで、\( (4) \) 式をある整数 \( N \) についてのみについて書くと次式になります。 \[ \begin{align} c_N \left( t + k \right) e^{iNpx} = c_N \left( t \right) e^{iNpx} - \frac{ck}{h} \left\{ c_N \left( t \right) e^{iNp \left( x + h \right)} - c_N \left( t \right) e^{iNpx} \right\} \ \ \ \ldots (5) \end{align}\] \( (5) \) 式を変形して、時間ステップが進む前後でのフーリエ係数の比 \( c_N \left( t + k \right) / c_N \left( t \right) \) を求めます。 \[ \begin{align} \frac{c_N \left( t + k \right)}{c_N \left( t \right)} &= 1 - \frac{ck}{h} \left( e^{iNph} - 1 \right) \\\\ &= 1 - \frac{ck}{h} \left( \cos Nph + i \sin Nph - 1 \right) \\\\ &= 1 - \frac{ck}{h} \left( \cos Nph - 1 \right) - i \frac{ck}{h} \sin Nph \end{align}\] 丸め誤差が限りなく大きくならないためには、この式の値の絶対値が1以下となればよいです。つまり、次の不等式が成り立てばよいです。 \[ \begin{align} -1 \leq \sqrt{ \left\{ 1 - \frac{ck}{h} \left( \cos Nph - 1 \right) \right\}^2 + \frac{c^2k^2}{h^2} \sin ^2 Nph } \leq 1 \end{align}\] 左側の不等号は明らかに成り立つため、右側の不等式を両辺2乗した次の不等式 \( (6) \) だけを考えます。 \[ \begin{align} \sqrt{ \left\{ 1 - \frac{ck}{h} \left( \cos Nph - 1 \right) \right\}^2 + \frac{c^2k^2}{h^2} \sin ^2 Nph } & \leq 1 \\\\ \left\{ 1 - \frac{ck}{h} \left( \cos Nph - 1 \right) \right\}^2 + \frac{c^2k^2}{h^2} \sin ^2 Nph & \leq 1 \ \ \ \ldots (6) \end{align}\] 整理していきます。 \[ \begin{align} & \ \ \ \ \ \ \left\{ 1 - \frac{ck}{h} \left( \cos Nph - 1 \right) \right\}^2 + \frac{c^2k^2}{h^2} \sin ^2 Nph \\\\ &= 1 - \frac{2ck}{h} \left( \cos Nph - 1 \right) + \frac{c^2k^2}{h^2} \left( \cos Nph - 1 \right) ^2 + \frac{c^2k^2}{h^2} \sin ^2 Nph \\\\ &= 1 - \frac{2ck}{h} \left( \cos Nph - 1 \right) \\\\ & \ \ \ \ \ \ + \frac{c^2k^2}{h^2} \cos ^2 Nph - \frac{2c^2k^2}{h^2} \cos Nph + \frac{c^2k^2}{h^2} + \frac{c^2k^2}{h^2} \sin ^2 Nph \\\\ &= 1 - \frac{2ck}{h} \left( \cos Nph - 1 \right) - \frac{2c^2k^2}{h^2} \left( \cos Nph - 1 \right) \\\\ &= 1 - \frac{2ck}{h} \left( 1 + \frac{ck}{h} \right) \left( \cos Nph - 1 \right) \leq 1 \\\\ & \ \ \ \ \ \ \ \ \ \ \ \ \ \ \frac{2ck}{h} \left( 1 + \frac{ck}{h} \right) \left( \cos Nph - 1 \right) \geq 0 \ \ \ \ldots(7) \end{align}\] ここで、\( -2 \leq \cos Nph - 1 \leq 0 \) であるため、\( c \gt 0 \) の場合、\( (7) \) の不等式は \( \cos Nph - 1 = 0 \) のときのみ成立しますが、それ以外では成り立ちません。 そのため、\( c \gt 0 \) の場合、この陽的差分法の手法は無条件に不安定であるといえます。 一方、\( c \lt 0 \) の場合、次の不等式が成り立つときに限り \( (7) \) の不等式が成り立ちます。 \[ \begin{align} 1 + \frac{ck}{h} & \geq 0 \\\\ - \frac{ck}{h} & \leq 1 \ \ \ \ldots (8) \end{align}\] 従って、\( c \lt 0 \) であり、かつ、\( (8) \) の不等式が成り立つように3つ組 \( \left( c, k, h \right) \) が定められていれば、この陽的差分法の手法は安定であるといえます。

　なお、差分方程式の空間方向の離散化を後退差分近似で行った場合は、\( c \gt 0 \) であり、かつ、 \[ \begin{align} \frac{ck}{h} \leq 1 \end{align}\] が成り立つときにこの手法は安定となります。

　2次元拡散方程式は次の偏微分方程式です。 \[ \begin{align} \frac{\partial u}{\partial t} = \kappa \left( \frac{\partial ^2 u}{\partial x^2} + \frac{\partial ^2 u}{\partial y^2} \right) \ \ \ \left( \kappa \gt 0 \right) \end{align}\] この偏微分方程式を時間の刻み幅を \( k \) 、空間の刻み幅を \( h \) として、時間について前進差分近似、空間について中心差分近似の差分方程式の形で書き下すと次式になります。 \[ \begin{align} u \left( x, y, t + k \right) = \ & u \left( x, y, t \right) + \frac{\kappa k}{h^2} \{ u \left( x + h,y,t \right) + u \left( x-h,y,t \right) \\\\ & + u \left( x,y+h,t \right) + u \left( x,y-h,t \right) - 4u \left( x,y,t \right) \} \ \ \ \ldots \left( *1 \right) \end{align}\] 丸め誤差を \( \tilde{u} \left( x, y, t \right) \) と表せば、これは \( \left( *1 \right) \) 式と同じ差分方程式に従って変化します。すなわち、次式が成り立ちます。 \[ \begin{align} \tilde{u} \left( x, y, t + k \right) = \ & \tilde{u} \left( x, y, t \right) + \frac{\kappa k}{h^2} \{ \tilde{u} \left( x + h,y,t \right) + \tilde{u} \left( x-h,y,t \right) \\\\ & + \tilde{u} \left( x,y+h,t \right) + \tilde{u} \left( x,y-h,t \right) - 4\tilde{u} \left( x,y,t \right) \} \ \ \ \ldots \left( *2 \right) \end{align}\] 時刻 \( t \) における丸め誤差は2重フーリエ級数を用いて表すことができます。 \[ \begin{align} \tilde{u} \left( x, y, t \right) = \sum ^{\infty} _{m = - \infty} \sum ^{\infty} _{n = - \infty} c_{nm} \left( t \right) e^{i \left( npx + mqy \right)} \ \ \ \ldots \left( *3 \right) \end{align}\] ここで、\( p \) と \( q \) は計算対象としている空間領域に応じた任意の実数です。 先ほど同様、ある整数の組 \( \left( N, M \right) \) についての2重フーリエ係数 \( c_{NM} \) が時間発展につれて限りなく大きくなるか否かを個別に調べても一般性を失いません。 そこで、\( \left( *3 \right) \) 式に含まれるある整数の組 \( \left( N, M \right) \) についての項のみを \( \left( *2 \right) \) 式に適用すると次式になります。 \[ \begin{align} & \ \ \ \ \ \ c_{NM} \left( t + k \right) e^{i \left( Npx + Mqy \right)} \\\\ &= c_{NM} \left( t \right) e^{i \left( Npx + Mqy \right)} \\\\ & \ \ \ \ + \frac{\kappa k}{h^2} \{ c_{NM} \left( t \right) e^{i \left( Np \left( x + h \right) + Mqy \right)} + c_{NM} \left( t \right) e^{i \left( Np \left( x - h \right) + Mqy \right)} \\\\ & \ \ \ \ + c_{NM} \left( t \right) e^{i \left( Npx + Mq \left( y + h \right) \right)} + c_{NM} \left( t \right) e^{i \left( Npx + Mq \left( y - h \right) \right)} \\\\ & \ \ \ \ - 4 c_{NM} \left( t \right) e^{i \left( Npx + Mqy \right)} \} \ \ \ \ldots \left( *4 \right) \end{align}\] \( (*4) \) 式を変形して、時間ステップが進む前後での2重フーリエ係数の比 \( c_{NM} \left( t + k \right) / c_{NM} \left( t \right) \) を求めます。 \[ \begin{align} \frac{c_{NM} \left( t + k \right)}{c_{NM} \left( t \right)} &= 1 + \frac{\kappa k}{h^2} \left( e^{iNph} + e^{-iNph} + e^{iMqh} + e^{-iMqh} - 4 \right) \\\\ &= 1 + \frac{\kappa k}{h^2} \left( 2 \cos Nph + 2 \cos Mqh - 4 \right) \\\\ &= 1 + \frac{2 \kappa k}{h^2} \left( \cos Nph + \cos Mqh - 2 \right) \end{align}\] 丸め誤差が限りなく大きくならないためには、この式の値の絶対値が1以下となればよいです。つまり、次の不等式が成り立てばよいです。 \[ \begin{align} -1 \leq 1 + \frac{2 \kappa k}{h^2} \left( \cos Nph + \cos Mqh - 2 \right) \leq 1 \\\\ -2 \leq \frac{2 \kappa k}{h^2} \left( \cos Nph + \cos Mqh - 2 \right) \leq 0 \end{align}\] \( \cos Nph \leq 1 \) かつ \( \cos Mph \leq 1 \) なので、上の式の右側の不等号は常に成り立ちます。 また、\( \cos Nph \geq -1 \) かつ \( \cos Mph \geq -1 \) なので、次の不等式 \( \left( *5 \right) \) が安定性の条件式になります。 \[ \begin{align} -2 \leq - \frac{8 \kappa k}{h^2} \\\\ \frac{\kappa k}{h^2} \leq \frac{1}{4} \ \ \ \ldots \left( *5 \right) \end{align}\]

　なお、ここまで解説してきた Von Neumann の安定性解析ですが、この方法は差分法における境界条件の影響を考慮していません。 この方法で得られる条件を満たしていても数値解が不安定な振る舞いをしたり、逆に条件を満たしていなくてもうまく計算できることもあります。 ただ、計算を安定的に行う上での実用的な指針にはなるため、Von Neumann の安定性解析で得られる条件を参考に適切な刻み幅を検討することはよい方法だと思います。

　終わりに、John von Neumann（1903-1957）はハンガリー出身のアメリカの数学者です。 数学、物理学をはじめ、経済学や気象学など様々な科学の分野に貢献した人物です。 驚異的な頭脳を持っていたと言われており、8桁×8桁のかけ算と割り算を暗算で行ったなどの様々な逸話が残っています。