フォン・ノイマンの安定性解析

フォン・ノイマンの安定性解析（フォン・ノイマンのあんていせいかいせき、英: Von Neumann stability analysis）とは、数値解析において、線形偏微分方程式を有限差分法で解く際の数値的安定性を調べるのに使われる手法である。この手法は数値的誤差のフーリエ展開に基づいており、ジョン・クランク（英語版）とフィリス・ニコルソン（英語版） (1947) によって簡潔に述べられた後、ロスアラモス国立研究所によって発展された。後に、この方法はフォン・ノイマンとの共著により、より厳密に取り扱われた。

この手法は厳密には、等間隔格子上の線形の連立方程式に対する初期値問題にのみ適用できる。これは一見厳しい制限に見えるが、経験的にはこの解析は信頼できる結果を示し、より一般的な問題に対する指針となっている[1]。

数値的安定性

「数値的安定性」も参照

数値計算手法の安定性は数値誤差に密接にかかわっている。計算中のある時間ステップで生じた誤差が計算を続けるにあたって増大しないならば、有限差分法は安定である。計算を続けていくと誤差が一定のまま残るときは中立安定と言われる。誤差が減衰し、最終的に消失するなら、その計算手法は安定と呼ばれる。逆に、誤差が時間ステップとともに増大した場合、数値スキームは不安定であると言われる。数値スキームの安定性は、フォン・ノイマンの安定性解析によって調べることができる。時間依存の問題に対してスキームが安定であることは、厳密な微分方程式の解が有界であるならば、数値計算法も有界な解を生成することを保証する。一般に、非線形偏微分方程式である場合は特に、安定性を調べることが困難になる。

以下に挙げる特定のケースでは、フォン・ノイマンの安定性はラックス・リヒトマイヤーの意味での安定性（ラックスの等価定理で用いられる）の必要十分条件である：

  • 偏微分方程式および有限差分スキームモデルが線形である
  • 偏微分方程式は定数係数で周期的境界条件および 2 つの独立変数（時間と空間）を持っている
  • スキームは 2 つより多くの時間ステップを用いない。

フォン・ノイマンの安定性は例より多くの種類のケースで必要である。比較的単純であるために、それはしばしば、スキームで使用されるステップサイズの制限（もしあれば）についての良い推測をするために、より詳細な安定性解析の代わりに使用される。

解析手法

フォン・ノイマンの安定性解析は誤差のフーリエ分解に基づいている。ここでは1次元の熱伝導方程式：

∂ u ∂ t = α ∂ 2 u ∂ x 2 {\displaystyle {\frac {\partial u}{\partial t}}=\alpha {\frac {\partial ^{2}u}{\partial x^{2}}}}

をFTCS法（英語版）[注 1]を用い離散化した次の式の安定性を考える。

u j n + 1 = u j n + r ( u j + 1 n − 2 u j n + u j − 1 n ) {\displaystyle u_{j}^{n+1}=u_{j}^{n}+r\left(u_{j+1}^{n}-2u_{j}^{n}+u_{j-1}^{n}\right)} 　・・・(1)

ただしr は拡散数

r = α Δ t Δ x 2 {\displaystyle r={\frac {\alpha \,\Delta t}{\Delta x^{2}}}}

で、区間の長さをL とする。差分方程式の解 u j n {\displaystyle u_{j}^{n}} は格子上で偏微分方程式の解析解 u ( x , t ) {\displaystyle u(x,t)} を近似する。

丸め誤差 ϵ j n {\displaystyle \epsilon _{j}^{n}} を

ϵ j n = N j n − u j n {\displaystyle \epsilon _{j}^{n}=N_{j}^{n}-u_{j}^{n}}

と定義する。ただし u j n {\displaystyle u_{j}^{n}} は差分方程式(1)を丸め誤差なしで計算したときの解で、 N j n {\displaystyle N_{j}^{n}} は有限精度計算で得られた数値解である。厳密解 u j n {\displaystyle u_{j}^{n}} は差分方程式を厳密に満たすから、誤差 ϵ j n {\displaystyle \epsilon _{j}^{n}} もまた差分方程式を厳密に満たす。したがって誤差は次の漸化式を満たす。

ϵ j n + 1 = ϵ j n + r ( ϵ j + 1 n − 2 ϵ j n + ϵ j − 1 n ) {\displaystyle \epsilon _{j}^{n+1}=\epsilon _{j}^{n}+r\left(\epsilon _{j+1}^{n}-2\epsilon _{j}^{n}+\epsilon _{j-1}^{n}\right)} 　・・・(2)

式(1), (2)は、誤差と数値解の両方が時間ステップに応じて同じように成長または減衰することを示す。周期境界条件を持つ線形微分方程式に対し、誤差の空間的変動は区間L で次のようにフーリエ級数に展開できる：

ϵ ( x ) = ∑ m = 1 M A m e i k m x {\displaystyle \epsilon (x)=\sum _{m=1}^{M}A_{m}e^{ik_{m}x}}

ここで

k m = π m L {\displaystyle k_{m}={\frac {\pi m}{L}}} ：波数
M = L / Δ x {\displaystyle M=L/\Delta x} ：分割数

である。誤差の時間依存性は誤差の振幅 A m {\displaystyle A_{m}} が時間ステップの関数である仮定することによって考慮されている。誤差の成長・減衰は指数関数的になる傾向があるので、振幅が時間とともに指数関数的に変化すると仮定するのは妥当である。ゆえに

ϵ ( x , t ) = ∑ m = 1 M e a t e i k m x {\displaystyle \epsilon (x,t)=\sum _{m=1}^{M}e^{at}e^{ik_{m}x}}

と仮定する。ただしa は定数である。

誤差が従う差分方程式は線形なので（級数の各項の挙動は級数自体と同じである）、次の典型的な項の誤差の成長を考察すれば十分である：

ϵ m ( x , t ) = e a t e i k m x {\displaystyle \epsilon _{m}(x,t)=e^{at}e^{ik_{m}x}} 　・・・(3)

誤差に対するこの形式を使用して安定特性を調べても一般性を失わない。誤差が時間ステップを進めるごとにどのように変化するかを調べるため、式(3)を(2)に代入し整理すると

e a Δ t = 1 + r ( e i k m Δ x + e − i k m Δ x − 2 ) = 1 − 4 r sin 2 ⁡ k m Δ x 2 {\displaystyle {\begin{aligned}e^{a\Delta t}&=1+r\left(e^{ik_{m}\Delta x}+e^{-ik_{m}\Delta x}-2\right)\\&=1-4r\sin ^{2}{\frac {k_{m}\Delta x}{2}}\end{aligned}}} 　・・・(4)

を得る[注 2]。

振幅係数G を

G ≡ ϵ j n + 1 ϵ j n = e a Δ t {\displaystyle G\equiv {\frac {\epsilon _{j}^{n+1}}{\epsilon _{j}^{n}}}=e^{a\Delta t}} 　・・・(5)

と定義する。誤差が有界であるための必要十分条件は | G | ≤ 1 {\displaystyle \vert G\vert \leq 1} である。したがって式(4), (5)より、安定性の条件は

| 1 − 4 r sin 2 ⁡ k m Δ x 2 | ≤ 1 {\displaystyle \left\vert 1-4r\sin ^{2}{\frac {k_{m}\Delta x}{2}}\right\vert \leq 1}

と与えられる。この条件が任意の k m {\displaystyle k_{m}} について成り立たなければならないから、

r = α Δ t Δ x 2 ≤ 1 2 {\displaystyle r={\frac {\alpha \,\Delta t}{\Delta x^{2}}}\leq {\frac {1}{2}}} 　・・・(6)

を得る。式(6)は1次元熱伝導方程式をFTCS法で解くときの、安定性の必要条件を与える。与えられた空間ステップ幅 Δ x {\displaystyle \Delta x} に対して、時間ステップ幅 Δ t {\displaystyle \Delta t} は式(6)を満たすように十分に小さく取らなければならないことが分かる。

脚注

  1. ^ FTCS法（forward in time and central difference in space method）とは、時間微分の離散化に前進差分法を、空間微分の離散化に中心差分法を用いる離散化法である。熱伝導方程式をFTCS法で離散化すると
    u j n + 1 − u j n Δ t = α u j + 1 n − 2 u j n + u j − 1 n Δ x 2 {\displaystyle {\frac {u_{j}^{n+1}-u_{j}^{n}}{\Delta t}}=\alpha {\frac {u_{j+1}^{n}-2u_{j}^{n}+u_{j-1}^{n}}{\Delta x^{2}}}}
    となる。この式を変形すると式(1)を得る。(藤井(2014), p. 15)
  2. ^ 式(4)の導出には、関係式
    ϵ j n = e a t e i k m x ϵ j n + 1 = e a ( t + Δ t ) e i k m x = e a Δ t ϵ j n ϵ j ± 1 n = e a t e i k m ( x ± Δ x ) = e ± i k m Δ x ϵ j n {\displaystyle {\begin{aligned}\epsilon _{j}^{n}&=e^{at}e^{ik_{m}x}\\\epsilon _{j}^{n+1}&=e^{a(t+\Delta t)}e^{ik_{m}x}=e^{a\Delta t}\epsilon _{j}^{n}\\\epsilon _{j\pm 1}^{n}&=e^{at}e^{ik_{m}(x\pm \Delta x)}=e^{\pm ik_{m}\Delta x}\epsilon _{j}^{n}\end{aligned}}}
    および次の恒等式を用いる：
    cos ⁡ ( k m Δ x ) = e i k m Δ x + e − i k m Δ x 2 , sin 2 ⁡ k m Δ x 2 = 1 − cos ⁡ ( k m Δ x ) 2 {\displaystyle {\begin{aligned}\cos(k_{m}\Delta x)&={\frac {e^{ik_{m}\Delta x}+e^{-ik_{m}\Delta x}}{2}},\\\sin ^{2}{\frac {k_{m}\Delta x}{2}}&={\frac {1-\cos(k_{m}\Delta x)}{2}}\end{aligned}}}
  1. ^ 藤井(2014), p. 18.

参考文献

  • 藤井孝藏『流体力学の数値計算法』東京大学出版会、1994年。ISBN 4-13-062802-X。