フォン・ノイマンの安定性解析 (フォン・ノイマンのあんていせいかいせき、英 : Von Neumann stability analysis )とは、数値解析 において、線形 偏微分方程式 を有限差分法 で解く際の数値的安定性 を調べるのに使われる手法である。この手法は数値的誤差のフーリエ展開に基づいており、ジョン・クランク(英語版 )  とフィリス・ニコルソン(英語版 )   (1947) によって簡潔に述べられた後、ロスアラモス国立研究所 によって発展された。後に、この方法はフォン・ノイマン との共著により、より厳密に取り扱われた。 
  この手法は厳密には、等間隔格子上の線形の連立方程式に対する初期値問題にのみ適用できる。これは一見厳しい制限に見えるが、経験的にはこの解析は信頼できる結果を示し、より一般的な問題に対する指針となっている。 
 
  数値的安定性 数値計算手法の安定性は数値誤差に密接にかかわっている。計算中のある時間ステップで生じた誤差が計算を続けるにあたって増大しないならば、有限差分法は安定である。計算を続けていくと誤差が一定のまま残るときは中立安定 と言われる。誤差が減衰し、最終的に消失するなら、その計算手法は安定 と呼ばれる。逆に、誤差が時間ステップとともに増大した場合、数値スキームは不安定 であると言われる。数値スキームの安定性は、フォン・ノイマンの安定性解析によって調べることができる。時間依存の問題に対してスキームが安定であることは、厳密な微分方程式の解が有界であるならば、数値計算法も有界な解を生成することを保証する。一般に、非線形偏微分方程式である場合は特に、安定性を調べることが困難になる。 
  以下に挙げる特定のケースでは、フォン・ノイマンの安定性はラックス・リヒトマイヤーの意味での安定性(ラックスの等価定理 で用いられる)の必要十分条件である: 
 
  
   偏微分方程式および有限差分スキームモデルが線形である  
   偏微分方程式は定数係数で周期的境界条件および 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)}
          
         
         
    
  ただし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}}
         
        
       
       
        
         
          
          
           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)}
          
         
         
    
  式(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)を(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}}}
          
         
         
    
  を得る[注 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}}
          
         
         
    
  と定義する。誤差が有界であるための必要十分条件は 
       
        
         
         
          |
           
         
          G
           
         
          |
           
         
          ≤
           
         
          1
           
          
         
       
        {\displaystyle \vert G\vert \leq 1}
         
        
       
 
  
   
    
        
         
          
           
           
            |
             
            
            
             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)は1次元熱伝導方程式をFTCS法で解くときの、安定性の必要条件を与える。与えられた空間ステップ幅
       
        
         
         
          Δ
           
         
          x
           
          
         
       
        {\displaystyle \Delta x}
         
        
       
       
        
         
         
          Δ
           
         
          t
           
          
         
       
        {\displaystyle \Delta t}
         
        
       
 
  脚注  
    
    ^ 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)  ^ 式(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}}}
              
             
             
          
    
   
  参考文献