Zenn
💬

数値シミュレーションのための代数方程式計算手法まとめ

2022/07/11に公開
math
#
simulation
#
OpenFOAM
tech

数値シミュレーション（流体，構造）を想定し，代数方程式（連立一次方程式）の数値的解法についてその全体像をざっくり紹介したいと思います．

数値シミュレーションと代数方程式（連立一次方程式）

数値シミュレーション（流体解析や構造解析）では，その過程で下記の方程式（連立一次方程式）を数値的に解くことが求められます．（\bm xが変数）
例えば，静的な構造解析では\bm Aが剛性マトリックス，\bm xが変位，\bm bが外力（境界条件含む）に対応します．この方程式のことを代数方程式と呼びます．

\bm A \bm x = \bm b

有限差分法，有限要素法，有限体積法など数値シミュレーション手法では空間を離散化して扱う必要があり，大雑把には「離散化された各位置（要素，点）の数」×「自由度の数」＝「\bm xの要素数」となります．解析したい物理現象にもよりますが，求める解の次元は，構造解析で数千から数万，流体解析では数万から数百万程度にもなります．

解法は，「消去法」，「反復法」，「共役勾配法」に大別されます．（共役勾配方は反復法の一部とみなすこともできますが，ここでは特徴を際立たせるために分けています．）

それぞれについて，概要，メリット，デメリット，代表的な利用先を紹介していきたいと思います．

直接法

概要

  • 連立方程式同士の和や差をとっても解が同じであることを利用して，行列\bm Aの要素を部分的に消去していく方法．
  • 代表的な方法：
    • ガウスの消去法
    • LU分解法
    • コレスキー分解法

LU分解の場合を例に取ると，

\bm A \bm x = \bm b

を以下の形に変形します．

\bm L \bm U \bm x = \bm b

\bm Lは下三角行列，\bm Uは上三角行列を表します．その後，前進代入，後退代入を手続きを経ることで解を求めます．
このように，行列Aの形の形を「解きやすい」形に変形することで解を求めます．

メリット

  • 一意の解をもつような連立一次方程式の場合には，有限回の計算で必ず解を求めることができる．（解が存在すれば）
  • LU分解の場合，\bm bが変更になったとしても，同じL,U行列を使い回すことができる．（\bm bの変更に対して強い．）

デメリット

  • メモリ使用量や計算回数が他の方法と比較して多い．
  • 特に疎行列の場合，値が0だった場所にも値が代入され確保すべきメモリ量が増大する．（fill in 問題という．）
    • ここでは深く言及しないが，fill in問題を回避（軽減）するためには不完全LU分解などを用いる．
  • 並列計算を行いにくい．

代表的な利用例

  • 数百から数千節点程度の構造解析

計算量が多くやメモリを確保する必要があるものの，有限回の計算で必ず求解できることから，数千から数万元程度の連立一次方程式を扱うことが多い構造解析にて使われることが多いです．

幾何非線形を無視した解析の場合，「LU分解が発生する」＝「どこか部材が降伏した」ことを表します．
「部材が降伏」⇒「剛性行列が変化する」⇒「係数行列が変化する」⇒「LU分解が必要」という理屈です．

反復法

概要

  • 適当な初期値を決めた後，反復計算をしながら真の初期値に近づけていく方法
  • 代表的な方法：
    • ヤコビ法
    • ガウス・ザイデル法（Gauss-Seidel法）
\bm A \bm x = \bm b

を以下の形に変形します．

\bm x = \bm M \bm x + \bm c

ここから，適当な\bm x^{(0)}からスタートして以下の式に従い，\bm x^{(n)}を計算します．

\bm x^{(n + 1)} = \bm M \bm x^{(n)} + \bm c

以下が成立するまで（＝残差が小さくなるまで）繰り返します．

|\bm X^{(n + 1)} - \bm X^{(n)}| < \epsilon \space 又は \space \frac{|\bm X^{(n + 1)} - \bm X^{(n)}|}{|\bm X^{(n)}|} < \varepsilon

ここで\epsilon, \varepsilonは予め与えられた正数とします．

メリット

  • 直接法と比較してメモリ使用量が少ない．
  • 出発点を適切に選ぶことができれば，計算時間が短くなる．
  • fill-in 問題が発生しない．
  • 並列計算を行いやすい．

デメリット

  • 収束するまでの反復回数を見積もることができない．
  • 係数行列によっては収束しないことがある．

代表的な利用先

  • 数千から数万節点程度の流体解析

係数行列（\bm A）が非対称になる輸送方程式を解く場合，ガウス・ザイデル法などの反復法が使用されます．（OpenFOAMのチュートリアルでも頻出です．）

共役勾配法

概要

  • Aが正定値対称行列の場合に適用する反復解法．
  • 最急降下法の一種であり，CG法（Conjugate Gradient method）とも呼ばれる．
  • Aが正定値対称行列でない場合には，双対共役勾配法（BiCG法）を用いる．

以下の関数 f(\bm x)を考える．

f(\bm x) = (1/2) \bm x^{T} \bm A \bm x - \bm b \cdot \bm x

\bm Aが正定値対称行列の場合，\bm A \bm x = \bm bの解がf(\bm x)を最小にすることに注目します．
この性質を利用し，適当な出発点からスタートして，\bm z = f(\bm x)が表す曲面の最大傾斜方向に進み（谷の方向へ進む），f(\bm x)を最小にする\bm xを見つけることで，\bm A \bm x = \bm bの解を間接的に求める方法です．

メリット

  • ｎ元の方程式に対し，理論上ｎ回以下の反復計算で厳密解を求めることができる．（実用上はｎ回以上の反復回数が必要になる）
  • 直接法よりも計算回数が少ない．
  • fill-in問題を回避しやすい．(前処理によっては部分的にfill-inが発生することもある．)

デメリット

  • \bm z = f(\bm x)が表す曲面の形状が「悪い」場合，初期値によって収束する解が異なる可能性がある．
    • 形状が「悪い」とは，収束しにくい曲面のことを表す．例えば，特定の方向に扁平な曲面など．

代表的な利用先

  • 数千から数百万節点程度の流体解析

係数行列（ \bm A ）が対称になるポアソン方程式を解く場合，共役勾配法が使用されることが多いです．係数行列（ \bm A ）が非対称になる輸送方程式などでは，非対称係数行列に対応できる双対共役勾配法（BiCG法）が用いられます．
扱う次元数が多くても，直接法よりもメモリと計算量を小さくでき，反復法よりも安定して求解できるのが魅力です．

まとめ

「直接法」，「反復法」，「共役勾配法」の違いを下表にまとめました．基本的には，求める解の次元数が小さいうちは「直接法」，大きくなるにつれて，反復法，共役勾配法を利用します．

直接法 反復法 共役勾配法
計算量 大 小〜大 中
メモリ使用量 大 小 小
fill-in 問題 あり なし （基本的に）なし
応用先 構造解析 流体計算（輸送方程式） 流体計算（輸送方程式，ポアソン方程式）

参考文献

Discussion