6. 数値シミュレーション
この例では図1に示されるような図形をメッシュ分割し、そのメッシュを用いてHelmholtz方程式の解を計算する。
|
図1 幾何形状の構成 |
この図形は3つの部分領域に分割され、別個にメッシュ分割を施した後、D06DBFを2回コールすることによって縫い合わされる。最初の2つの図形に対する境界はD06BAFを用いて生成される。最初の図形は(D06ABFを用い)Delaunay Voronoiプロセスによって、第2の図形は(D06ACFを用い)前進フロントアルゴリズム(advancing front algorithm)によって、第3の図形は(D06DAFを用い)第2の図形をだけ回転させることによってメッシュ分割が行われる。 結果として生成された図形はD06CAFによって平滑化される。さらに頂点に対する番号の再設定、及びそのメッシュ構成に付帯する行列の疎構造(sparsity structure)がD06CCFによって生成される。
|
図2 メッシュ構造 |
生成された図形には2718個の頂点と5113個の三角形が、境界と内部の接合部には361本の辺が含まれている。ただし図2が示すのは粗いメッシュであり、それには643個の頂点、1133個の三角形、171本の辺しか含まれていない。
外部境界全体(その辺はEDGE(3,:)を介して1, 2, 3とタグ付けされている)はディリクレ型境界として、内部境界(その辺は5とタグ付けされている)はノイマン型境界として設定されている。
行列の正当性をチェックするために、下記のプログラム中では数式(5), (6)が積分公式(9)との絡みで検証されている。このテストは数値計算の面で効率的でない場合もあるが、行列の正当性に関する感触を得ることができる。次数が1以下の多項式に対しては結果は厳密値と一致するため(データファイル中でTYPE = 'P1'を選択する)、より良いチェックを行うことができる。
の典型的な値としては1が考えられるため、そのためのテストとしてが想定されている (TYPE = 'P2')。疎な線形系の解法についてはいくつかの選択肢が用意されている。
-
SOLUT = 'D'は直接法による解法機能を提供する。行列のLU分解はF11DAFを用いて行われる。その後、系はその分解を用い、F11DBFをコールする形で解が計算される。
-
SOLUT = 'I'は反復法による解法機能を提供する。行列の不完全なLU分解はF11DAFを用いて行われる。その後、系はF11DCFをコールする形で解が計算されるが、その際、次のいずれかの手法を用いる(いずれも不完全LU事前処理(incomplete LU preconditioning)を伴う)。
-
METHOD = 'RGMRES'(restarted generalised minimal residual method)
-
METHOD = 'CGS'(conjugate gradient squared method)
-
METHOD = 'BICGSTAB'(stabilised bi-conjugate gradient method)
-
METHOD = 'TFQMR'(transpose-free quasi-minimal residual method)
このプログラム中で使用されているD06, F11ルーチンの仕様については文献 [1] を参照されたい。
プログラムの最後の部分で数値解と真の値との誤差のノルム、及び勾配が計算されている(数式(7), (8), (9)を用いる)。相対誤差はオーダがであるのに対し、解の勾配に関する相対誤差はオーダがであることを示すことができる(文献 [4] 参照)。