ヒューリスティック解析

ヒューリスティック解析

有限差分方程式の計算結果で、本来、近似する偏微分方程式で想定されるものと大きく異なる、著しく増大して符号が頻繁に反転する解を得ることがあります。こうした解が示す挙動を「計算不安定性」と呼びます。当然ながら、こうした解は望ましくありません。定数係数を伴う線形差分方程式の計算安定性を確認する方法としては、von Neumannによるフーリエ手法を用いることができます(本シリーズ「計算安定性」を参照)。残念ながら、物理現象を表す大半の方程式は、非線形であるか、非定数係数を伴うか、もしくはその両方です。

ヒューリスティック解析手法

本書では、前述の有限差分方程式の計算安定性を調べるための、簡単なヒューリスティック解析手法について説明します。このタイプの解析は、多くのケースで、不安定性を排除する方法を示すだけでなく、近似の精度を高める方法も示すという優れた特長があります。

ここで説明するアプローチは、厳密でも完全でもないことから「ヒューリスティック」と呼ばれていますが、多くのケースに有効で、有用な情報を数多くもたらします。安定性を解析するためのヒューリスティック手法について書かれた参考文献[1]は、本書に記載する多くの情報の出典元です。

ヒューリスティック解析は、有限差分方程式を展開して各項をテイラー級数で表し、特定次数までの項のみを残すことで偏微分方程式に帰着させる、比較的簡単な考え方に基づいています。この展開は、初めは小さいものと想定される、空間増分および時間増分の累乗で表されます。

こうした展開では、元の偏微分方程式を最小次数まで再現することが不可欠です。そうでなければ、良好な近似は得られません。この必要条件は、しばしば近似の「一致性」と呼ばれます。展開された最小次数以降の項は、打ち切り誤差と呼ばれます。

ヒューリスティック解析は、テイラー展開された方程式の方が元の偏微分方程式よりも差分方程式をより高精度に表しているという基本概念に基づいています。打ち切り誤差の項をいくつか残した場合でも、項は差分方程式により近い偏微分方程式となります。この点を念頭に置きながら、ここでは計算を打ち切った式を調べることで、安定性に関する問題、必要な初期条件、深刻な不正確性など、差分方程式と共通の特性が明らかになることを示します。

まず、安定性について書かれた「計算安定性」で用いたのと同じ線形偏微分方程式について考えます。

線形方程式の例

ここでは、変数u(x,t)の1次の移流拡散方程式を用います。

(1) \frac{\partial u}{\partial t}+c\frac{\partial u}{\partial x}=\nu \frac{{{\partial }^{2}}u}{\partial {{x}^{2}}}

対流速度cおよび拡散係数νは定数であると想定します。この方程式の解は、有界であり、良好な挙動を示すことがわかっています。

ここでは、差分方程式のテイラー級数展開から生じる打ち切り誤差を調べることで、式1に対する簡潔な有限差分近似の安定性を判断できることを示します。このプロセスによって、不安定性には基本的に2つのタイプがあることが明らかになるだけでなく、ヒューリスティック手法と「計算安定性」で用いたvon Neumannタイプのフーリエ解析を直接比較できるようにもなります。こうした比較を行うことで、差分方程式の安定性を評価するうえで、テイラー展開で生じた打ち切り誤差のうち、保持すべき項と排除すべき項を決定するための有用な経験則を得ることができます。

以下の式は、「計算安定性」で説明した、式1を近似する簡潔かつ陽的な有限差分方程式です。

(2) \frac{u_{j}^{n+1}-u_{j}^{n}}{\delta t}=-\frac{c}{2\delta x}\left( u_{j+1}^{n}-u_{j-1}^{n} \right)+\frac{\nu }{\delta {{x}^{2}}}\left( u_{j+1}^{n}-2u_{j}^{n}+u_{j-1}^{n} \right)

ここで、ujnはu(jδx,nδt)を表します。これは、時間の前進差分近似と呼ばれるもので、タイムステップnにおける空間内の位置j値がすべて既知である場合に、ステップn+1におけるすべてのj値を計算することができます。つまり、元の偏微分方程式で1つの初期条件が必要であるのと同様に、単一の時間微分のみ伴うため、差分方程式でも計算を開始するにあたって初期条件が1つ必要となります。

差分方程式2では、空間位置および時間位置(jδx,nδt)ごとに、タイムステップn+1における位置j-1、j、j+1の各点に影響する特性が見られます。つまり、点(jδx,nδt)は、現在より先のある時間にて、x-t空間内で傾き±δx/δtを持つ線が境界となる影響領域を有します。これは、信号の伝搬を表す特性線と似ています。たとえば、元の式1は、擾乱の移流を示す、傾きcの特性線を有します。しかし、離散方程式における特性線は、物理特性を表すものではなく、特定点における値の変化によって差分方程式のデータ値が変化する領域を定義する、計算の特性を表します。

「計算安定性」では、フーリエ級数による手法を用いて、差分方程式2に対する3つの安定条件を導き出せることがわかりました。本書では、近似式2に関連する打ち切り誤差を調べることで得られる情報について説明します。

打ち切り誤差の評価

式2の各項は、xおよびtの連続微分可能な関数であると想定します。これにより、たとえば、uj+1,n,nはu(xj+δx,tn)となり、点(xj,tn)のまわりでδxの累乗でテイラー級数展開できます。式2のすべての項についてδxおよびδtで展開すると、以下の式を得ます。

(3) \frac{\partial u}{\partial t}+c\frac{\partial u}{\partial x}-\nu \frac{{{\partial }^{2}}u}{\partial {{x}^{2}}}=-\frac{1}{2}\delta t\frac{{{\partial }^{2}}u}{\partial {{t}^{2}}}+O\left( \delta {{x}^{2}},\delta {{t}^{2}} \right)

2次以上のδxおよびδtの項は、オーダー記号を用いてO(δx2 ,δt2)と記述されます。δxおよびδtがゼロに近づくとき、元の偏微分方程式1に帰着することから、これは一致近似であるといえます。

フーリエ解析と打ち切り誤差解析の比較

「計算安定性」では、以下の形式の典型的フーリエモード

P_{j}^{n}\propto {{r}^{n}}{{e}^{{ikxj}}}

これを、差分方程式2に代入することで、rを求める式を得ました。

(4) r=1-\left( \frac{ic\delta t}{\delta x} \right)\sin \left( k\delta x \right)-\left( \frac{2\nu \delta t}{\delta {{x}^{2}}} \right)\left[ 1-\cos \left( k\delta x \right) \right]

差分方程式の計算安定性を実現するには、rの絶対値を1.0以下にすることが必要です。

exp(i(kx+wt))の形式のフーリエモードを、計算を打ち切った式3に代入すると、r=exp(iwδt)となり、wδtの累乗で展開され、さらにsinとcosはkδxの累乗で展開されて、式4と同じ結果を得られることがわかります。式3で保持されたO(δx2 ,δt2)のように、2つの結果は同じであると確定されます。

しかし、この比較では、式4の実数部と虚数部で構成されたrの基本形式を保持するには、kδxの累乗で展開されたときに、少なくともsinおよびcosの最初の非ゼロ項を保持すべきであることも示されます。rの虚数部の最初の非ゼロ項は、sin(kδx)から導かれたもので、kδxに比例します。これは、式3のxについての1次導関数に対応します。rの実数部の最初の非ゼロ項(1を除く)は、cos(kδx)から導かれたもので、(kδx)2に比例します。これは、式3のxについての2次導関数に対応します。

これらの点から、計算を打ち切った式で振幅係数rの最小次数の実数部と虚数部を再現するには、打ち切り誤差で、独立した各変数について最小次数の偶関数および奇関数(導関数)を保持しなければならないという経験則が導かれます。式3では、δtに比例する1次の項は1つのみで、tについての2次導関数です。δxに比例する1次の項はありません。

計算を打ち切った式での安全性の評価

前述の経験則を用いると、計算を打ち切った式は以下のようになります。

(5) \frac{\delta t}{2}\frac{{{\partial }^{2}}u}{\partial {{t}^{2}}}+\frac{\partial u}{\partial t}+c\frac{\partial u}{\partial x}-\nu \frac{{{\partial }^{2}}u}{\partial {{x}^{2}}}=0

ここでまず注意すべき点は、この式は、元の偏微分方程式1と同じではないということです。ここで証明したいのは、式5は式1よりも有限差分方程式を良好に近似できる式であり、それによって差分方程式の安定性を示す特性に関する情報を得られるという点です。まさに、これが証明されます。

前述のように、差分方程式では、傾きdx/dt=±δx/δtを持つ線が境界となる影響領域に情報が伝搬されます。同様に、計算を打ち切った式5は、空間についての2次導関数および時間についての2次導関数によって双曲型(すなわち、波動)の特性を有し、有効波動速度は±(2ν/δt)½です。差分方程式で、計算を打ち切った式を近似するには、その影響領域が少なくとも計算を打ち切った式の影響領域を包含している必要があります。これにより、以下の条件が導き出されます。

(6) \frac{2\nu }{\delta t}\le {{\left( \frac{\delta x}{\delta t} \right)}^{2}} or \frac{2\nu \delta t}{\delta {{x}^{2}}}\le 1

Courant、Friedrichs、およびLewy [2]は、似たような影響領域に関する条件を用いていました。現在、これは「クーラン条件」と呼ばれており、1つの時間増分の間に波が伝搬する距離は1つの空間増分未満に制限されるというものです。クーラン条件が満たされない場合、符号の頻繁な反転や指数関数的な増大を伴う不安定性が生じます。条件式6は、まさに「計算安定性」でフーリエ解析から導き出した安定条件の1つです。

計算を打ち切った式5の2つの1次導関数項(移流項)から、以下に示す類似のクーランタイプ条件を推測できます。ここで、情報は速度cで伝搬します。

(7) \frac{c\delta t}{\delta x}\le 1

この安定条件も「計算安定性」で示したもので、満たされなかったときには、同様に符号の反転や増大を伴う不安定性が生じます。

3つ目の安定条件を導き出すには、まず、δt項を変換することで計算を打ち切った式を書き換えます。このとき、展開の1次の項が保持されるように、時間導関数の代わりに空間導関数を有するように変換します。これは、式3をtで微分し、1次以上の項をすべて無視します。

(8) \frac{{{\partial }^{2}}u}{\partial {{t}^{2}}}+c\frac{\partial }{\partial x}\frac{\partial u}{\partial t}-\nu \frac{{{\partial }^{2}}}{\partial {{x}^{2}}}\frac{\partial u}{\partial t}=O\left( \delta t \right)

次に、式1を用いて、この式のu/tの時間の1次導関数を置き換えて、以下の式を得ます。

(9) \frac{{{\partial }^{2}}u}{\partial {{t}^{2}}}={{c}^{2}}\frac{{{\partial }^{2}}u}{\partial {{x}^{2}}}-2c\nu \frac{{{\partial }^{3}}u}{\partial {{x}^{3}}}+{{\nu }^{2}}\frac{{{\partial }^{4}}u}{\partial {{x}^{4}}}+O\left( \delta t \right)

最後に、この結果を用いてδt項について、計算を打ち切った式5を書き換えます。

(10) \frac{\partial u}{\partial t}+c\frac{\partial u}{\partial x}=\left( \nu -\frac{{{c}^{2}}\delta t}{2} \right)\frac{{{\partial }^{2}}u}{\partial {{x}^{2}}}+c\nu \delta t\frac{{{\partial }^{3}}u}{\partial {{x}^{3}}}-\frac{{{\nu }^{2}}\delta t}{2}\frac{{{\partial }^{4}}u}{\partial {{x}^{4}}}

最後に得られた式は、元の有限差分方程式を点x=jδxおよびt=(n+½)δtのまわりでテイラー展開して(おそらく、より簡単に)得られるのと同じ式です。

前述の経験則から、δtに比例する右辺の最後の2つの項は、右辺の最初のδt項に含まれるものよりも高次の導関数を含むことから、破棄できます。

(11) \frac{\partial u}{\partial t}+c\frac{\partial u}{\partial x}=\left( \nu -\frac{{{c}^{2}}\delta t}{2} \right)\frac{{{\partial }^{2}}u}{\partial {{x}^{2}}}

これは、計算を打ち切った式の代替形式で、最小次数(1次)の打ち切り誤差と、独立した各変数について最小の偶関数および奇関数(導関数)を含むものだけを保持します。

式11は、変形させた拡散係数を除き、元の式1とほぼ同じです。ここで重要なのは、拡散係数は負となる可能性があることです。拡散係数が正である限り、式11の解は指数関数的な減衰挙動を示しますが、係数が負となる解では指数関数的に増大する特性が見られる、すなわち計算の不安定性が生じます。したがって、計算安定性を実現するもう1つの条件として、拡散係数が正となることを定めます。

(12) \frac{{{c}^{2}}\delta t}{2}\le \nu

このケースの不安定性は、前述の影響領域に関する2つの条件に関連した符号の反転を伴うものではなく、単に増大する特性のものです。不安定性が見受けられる場合、符号の頻繁な反転を伴うものかどうかを把握することで、影響領域に関する条件または負の拡散係数に関する条件のどちらが満たされていないかを識別することができます。こうした情報を把握できれば、不安定性を解消する方法を容易に見つけ出すことができます。

2次元の流れへの適用

実験室規模の水門の下を通過する2次元(x,z)の流れの例は、支配方程式の非線形性から生じる計算不安定性を調べるテストになります。この物理現象問題では、0.9フィートの高さまで水をせき止めている水門があります。水門の下流側(右側)の水深は0.14フィートです。重力は、-z方向(下向き)に32.2フィート/s2です。時間t=0に、水門は0.125フィート上昇し、水が下流側へ流れ込みます。図1に、ナビエストークスソルバ[3]を用いて取得した、t=0.35sにおける流れを示します。この例で使用したソルバは、不安定性を自動的に排除するよう最適化されているため、このケースでは不安定性は見受けられません。ただし、プログラムで最適化されていない設定を強制的に実行することも可能です。

図1. (左)水門の下を通過する流れ。不安定な挙動は見られない。図2. (左)タイムステップを小さくし、粘性なしで計算したときに生じる流れの不安定性

図1. (左)水門の下を通過する流れ。不安定な挙動は見られない。図2. (右)タイムステップを小さくし、粘性なしで計算したときに生じる流れの不安定性

不安定な挙動を実例で説明するために、まず、シミュレーションで使用した垂直速度式に対して実行したヒューリスティック解析を考察します。ここでは、z方向速度wに対する有効な拡散係数に焦点を当てており、他のすべての打ち切り誤差は無視します。

(13) \frac{\partial w}{\partial t}+u\frac{\partial w}{\partial x}+w\frac{\partial w}{\partial z}+\frac{\partial }{\partial z}\left( \frac{p}{\rho } \right)+g=\left( \nu +\frac{\alpha u\delta x}{2}-\frac{{{u}^{3}}\delta t}{2}-\frac{\delta {{x}^{2}}}{4}\frac{\partial u}{\partial x} \right)\frac{{{\partial }^{2}}w}{\partial {{x}^{2}}}+\left( \nu +\frac{\alpha w\delta z}{2}-\frac{{{w}^{2}}\delta t}{2}-\frac{\delta {{z}^{2}}}{2}\frac{\partial w}{\partial z} \right)\frac{{{\partial }^{2}}w}{\partial {{z}^{2}}}

xおよびｚ方向のwの拡散は、式13の右辺の2つの項で表されています。ここで、vは流体粘性、αはwのu移流を表す項(式13の左辺の第2項)の数値近似を修正するパラメータです。α=0のとき、移流の有限差分近似はwの位置を中心とした近似ですが、α=1のとき、上流側または「ドナーセル」による近似を使用します。

まず注意すべき点は、ν=0かつ中心差分近似を使用した場合(α=0)、2つの有効粘性係数の最小次数の項はδtに比例しており、負になります。これは明らかに不安定な挙動を導くもので、中心差分近似のよく知られた特性です。拡散係数を正の値に維持するために十分な粘性を加える手法も、安定性を得るうえでは確立された方法ですが、拡散が大きくなる危険性もあります。上流側での差分オプションα=1は、合理的な妥協案です。条件wδt<δxが満たされている限り、拡散係数は正であり(δx2およびδz2の項が小さい場合)、シミュレーションも安定します。

拡散係数のδx2およびδz2の項が小さくない場合、不安定な挙動が生じる可能性があります。これを示すために、粘性を0に設定し、上流側の差分量をα=0.05に減らします。負のδt項がa項より小さくなるように、非常に小さなタイムステップδt=0.00025を使用します。これらの設定で実行されたシミュレーションを図2に示します。水門の上流側でz速度における不安定性が生じています。図3に、その拡大図を示します(色はz速度の大きさを表す)。

この不安定性は、δx2項によって負となったx方向の拡散係数に起因します。水門の上流側の流れはz方向に圧縮していますが、x方向に膨張しているため負の値となります。つまり、この領域内では、δx2項のuのx方向導関数は正であり、正味では負の拡散係数となります。

この結論を確認するには、負のδx2項を補正するために、わずかに粘性を追加します(ν=0.0093)。図4では、このわずかな変化によって、流れが確かに安定したことがわかります。

この例では、元の方程式の非線形項に起因する打ち切り誤差の項が、差分方程式の計算安定性に影響することが示されました。このタイプの不安定性は、von Neumannタイプのフーリエ解析では見つけることはできません。最も重要なのは、問題となるような打ち切り誤差が存在すると判明したときに、この知識を用いて有限差分方程式を修正することで、そうした誤差を排除できることです。

図3.(左) 負のδx2項に起因する局所的に不安定な流れの拡大図。色は、ｚ速度を表す。図4. (右)負のδx2項を補正するためにわずかに粘性を追加した点を除き、図3と同様。

図3.(左) 負のδx2項に起因する局所的に不安定な流れの拡大図。色は、ｚ速度を表す。図4. (右)負のδx2項を補正するためにわずかに粘性を追加した点を除き、図3と同様。

まとめ

本書では、線形の有限差分方程式2に関連するすべての安定条件を、打ち切り誤差に対するヒューリスティックなアプローチによって特定可能であることを示しました。このアプローチは、不安定性を特定できるだけでなく、それを排除する方法も示します。たとえば、影響領域に関する条件が満たされない場合、タイムステップを小さくすることでしか問題を解決できませんが、負の拡散係数が存在する場合は、拡散を大きくして誤差を補正して安定性を取り戻す方法もあります。負の拡散誤差の原因を知ることは、この問題を回避できるよう、元の有限差分方程式をどのように修正するかという方法も示される可能性があります。

ヒューリスティックアプローチの最も重要な特長は、定数係数を伴う線形方程式に限定されない点です。これは、水門の下を通過する流れの例で示されました。計算を打ち切った式の近似式を立てるうえで、特別な仮定は必要ありませんでした。偏微分方程式を近似するための差分方程式を記述するのではなく、差分方程式を近似するための偏微分方程式を記述するという、単に逆の手順を行うことが目的でした。計算を打ち切った式を立てるための簡単な経験則についても説明しました。この近似式を使用して、解の不安定性につながる、影響領域に関する条件が満たされているか、さらに負の拡散係数が存在するかの2点を確認しました。

安定性に関するヒューリスティック解析について記述された参考文献[1]には、圧縮流れおよび非圧縮流れを伴う、さらにいくつかの流体力学シミュレーション例が示されています。また、ヒューリスティックアプローチを実際の非線形問題に適用する方法についても詳細に示されています。

参考文献

  1. C.W. Hirt、『Heuristic Stability Theory for Finite-Difference Equations』、J. Comp. Phys.、2、339 (1968)
  2. R. Courant、K.O. Friedricks、H. Lewy、『Math. Ann.』100、32 (1928)
  3. フローサイエンス(アメリカ、ニューメキシコ州サンタフェ)米国Flow Science, Inc. (Santa Fe, NM)から提供されている商用ソフトウェアパッケージFLOW-3D

^ back to top