とつとつとしてろうとせずスマートフォン用の表示で見るとつとつとしてろうとせずひまつぶしにどうぞ。この広告は、90日以上更新していないブログに表示しています。以下の記述は私的な備忘録であり、内容の正確さに関しては一切保証いたしませんので、悪しからず。次の線形方程式を考える。\mathbf{f}\end{align}ここで、$\mathbf{A}$は係数行列、$\mathbf{u}$は変数ベクトル、$\mathbf{f}$は係数ベクトルである。$\mathbf{A}$が疎行列であるなら、上式の計算は反復解法によって行う必要がある。そこで、$\mathbf{A}$が次のように近似できるとする。これを元の式に代入すると、反復計算の式を得る。\mathbf{f}\end{align}ここで、$k$は反復計算のステップ数である。\mathbf{u}^{(k)}$とおけば、次式の修正方程式を得る。\mathbf{\rho}^{(k)}\end{align}上式を解いた後、次式によって解を更新する。\mathbf{\delta}^{(k+1)}\end{align}反復法では式\eqref{eq:M}を解が収束するまで繰り返し計算する。そのため、$\mathbf{M}$が$\mathbf{A}$のよい近似であり、かつ式\eqref{eq:M}の求解が容易かつ効率的に行えることが望ましい。LU分解は一つの選択肢であろう。$\mathbf{A}$が疎行列である場合、LU分解によって得られる行列は元の行列の疎構造を継承せず、0だった部分に0でない値が現れる、いわゆる「fill-in」が生じる。これを防ぐため、行列の疎構造を保存しつつ近似的に分解を行う不完全LU分解(ILU分解)を適用することができる。標準的なILU分解は解の収束が比較的遅い。SIP法)を適用することができる。例として、次式のような直交座標系におけるPoisson方程式を考える。f\end{align}ここで、$u$は求める未知関数、$f$は既知関数とする。上式を等間隔格子上で2次精度中心差分(5点スキーム)によって離散化すると、次の差分方程式を得る。y$は各方向の格子幅であり、添字$i$は$x$方向、$j$は$y$方向のインデックスであるとする。上式は次のようにまとめることができる。f_{i,j}\end{align}ここで、係数$A$の添え字は格子点上の現在位置($P$)と東西南北($EWSN$)の隣接点における量であることを表す(下図を参照)。格子点上の変数を1次元配列に格納する順番は列優先(column-wise)とする。\cdots$と並べる(配列は左下隅から始める)。このような順番付けはプログラムの中の2重ループを組むときに$i$を外側、$j$を内側のループとして回すことに対応する。して得られる係数行列$\mathbf{A}$は0でない要素が5本の対角線上に$WSPNE$の順に並んだ構造を持つ帯行列となる(行優先で変数を格納した場合は$SWPEN$の順番になる)。次式は係数行列のイメージである。A_{P}\end{array}\right)\end{align}$\mathbf{A}$が帯行列であることから、$\mathbf{M}$を$\mathbf{A}$の対角行の内$WSP$成分と等しい値を持つ下三角行列$\mathbf{L}$と、$NE$成分が$\mathbf{A}$と等しく対角成分が1である(単位)上三角行列$\mathbf{U}$に分解することを考える。$\mathbf{L}$と$\mathbf{U}$の積を計算すると、$NW$と$SE$の部分に非零要素が現れる。eq}のように表すと次式のようになる。eq}\end{align}上式は式\eqref{eq:ite}に対応するが、これをそのまま計算しても収束が遅い。これは誤差項$\mathbf{N}$が最小化されていないためである。Stoneは元々の楕円形方程式が空間的に滑らかに分布する解を持つことを利用して、誤差項を次のように近似すれば収束を加速できることを示した。u_{i,j-1}\end{align}この近似は当然、一般の行列に対しては成り立たないため、SIP法を汎用の線形計算法として利用することはできないことに注意すべきである(別の近似式を用いることもできる)。$\alpha$が$1$に近いほど収束を速めることができるが、具体的な最適値は問題による。この工夫によって反復誤差の影響が反復行列のすべての要素に反映されることにより、標準的なILU分解よりも収束性を改善することができる。ILU分解と修正方程式の反復アルゴリズム具体的な計算は以下のように行う。\mathbf{u}^{(0)}$を計算する。1$の範囲でパラメータを設定する。反復計算の前に$\mathbf{L}$と$\mathbf{U}$を一度だけ計算する。これらは反復の度に更新する必要はない。L_{P}^{l}\end{align}ここで、$l$は行列ないしベクトルのインデックスである。上式の計算は上からこの順に実行しなければならない。境界付近では境界より外側の値は0として計算する。各反復において、残差と修正量の関係を表す修正方程式を解く。\mathbf{\rho}^{(k)}\end{align}上式はLU分解によって解くことができる。後退代入によって下式を解く。\mathbf{R}^{(k)}\end{align}$\mathbf{R}^{(k)}$は$l$の昇順にしたがって次のように計算できる。L_{P}^{l}\end{align}$\mathbf{R}^{(k)}$が求まれば、$\mathbf{\delta}^{(k+1)}$は前進代入によって$l$の降順で解くことができる。\delta^{l+N_{j}}\end{align}上式を解いた後、次式によって解を更新する。\mathbf{\delta}^{(k+1)}\end{align}プログラムで計算するときには、行列要素はすべて境界の外側の位置まで含めて配列を確保する。LU分解の要素の計算は次のように書ける。);}}例題テスト計算として次の極座標系におけるPoisson方程式を解いた。\theta}\end{align}解析条件は以下のとおり。計算結果のカラーマップ図を以下に示す。収束履歴は下図のようになった。残差の値は相対値ではなく各ステップ毎の減少量としてプロットした方がわかりやすかったように思うが、面倒なのでこれでよしとする。収束に要した反復回数は下表のようになった。190反復回数で見ると、GS法に比べてSIP法の方が10倍以上速い。いえ、計算負荷や使用メモリ量はSIP法の方が大きいため、状況によってはGS法(ないしそれを加速したSOR法)で十分という場合も考えられる。問題の性質や規模に応じた解法の選択が重要である。おわりに日本語で書かれたSIP法の解説がネットにはあまりないようなので備忘録として残すことにした。基本的には参考文献(1)ないし(2)の書き写しである。サンプルプログラムが手に入るURLなども載っているので、詳しくは(1)を参照されたい(10年以上前の本だがプログラムはまだ入手できる)。LU分解部分の具体的な計算は(2)にある。5点スキーム以外の問題の行列要素を求めるときに参考になるだろう(交差微分などを含む一般の楕円形方程式には9点スキームが必要)。NAG計算ライブラリにはSIP法が収録されているようなので、それを利用できる立場の人は自分でプログラムを作る必要はないかもしれない(もっというなら普通の人は共役勾配法系のライブラリを使うんじゃなかろうか......)。MSIP法)が提案されており、それらは2次元のとき9点スキーム、3次元のとき19点スキームを用いる。MSIP法はSIP法よりも収束性が向上しているらしい。詳細な解説のある文献があまり見つからなかったことと(特に(5)が読めなかった)、プログラミングが煩雑になりそうだったので、ここでは採用しなかった。Jaxはてなブログをはじめよう!iqujack-lequinaさんは、はてなブログを使っています。ブログを報告するスターをつけました引用をストックしましたストック一覧を見る閉じる引用するにはまずログインしてくださいログイン閉じる引用をストックできませんでした。再度お試しください閉じる限定公開記事のため引用できません。