(57)【要約】【課題】任意の格子分割数の問題に対して適用可能なマルチグリッド法を用いた連立一次方程式の数値計算方法及び装置を提供する。【解決手段】細かいグリッド上の値を粗いグリッド上の値へ変換する粗視化処理と、粗いグリッド上の値を細かいグリッド上の値へ変換する精細化処理と、において、変換先のグリッド上の格子を被覆する変換元のグリッド上の格子のグループを求め、該格子グループに属する各格子の大きさと、該格子と変換先の格子との共通部分の大きさと、の比率を重なり分率として求め、該格子グループに属する各格子における値の、変換先の格子における値への寄与を表す重み付けを設定し、該格子グループに属する全ての格子に亘って、各格子における値、重なり分率及び重み付けをの積の総和を求め、これを変換先の格子における値とする。細かいグリッドと粗いグリッドとの格子分割数比が非整数有理数の場合も粗視化及び精細化の処理が可能になる。【選択図】図１
2011154510.tif 000002
【特許請求の範囲】
【請求項１】
数値計算装置におけるマルチグリッド法を用いた連立一次方程式の数値計算方法であって、
前記数値計算装置が、
格子分割数Ｎｆの細かいグリッド上のｊ番目の格子ａｊにおける値ｘｊから、格子分割数Ｎｃ（＜Ｎｆ）の粗いグリッド上のｉ番目の格子Ａｉにおける値Ｘｉへの変換を
【数１】
2011154510.tif 000007
により行う粗視化工程と、
格子分割数Ｎｃの粗いグリッド上のｉ番目の格子Ｂｉにおける値Ｙｉから、格子分割数Ｎｆ（＞Ｎｃ）の細かいグリッド上のｊ番目の格子ｂｊにおける値ｙｊへの変換を
【数２】
2011154510.tif 000008
により行う精細化工程と、
を実行し、
前記粗視化工程は、
前記数値計算装置が、前記格子Ａｉを被覆する細かいグリッド上の格子ａｊの集合を格子グループ｛ａ｝ｉとして計算するする工程と、
前記数値計算装置が、前記格子グループ｛ａ｝ｉに属する格子ａｊ毎に、該格子ａｊと格子Ａｉとの共通部分の長さ、面積又は体積と、該格子ａｊの長さ、面積又は体積と、の比率を重なり分率ｐｉｊとして計算する工程と、
前記数値計算装置が、前記格子グループ｛ａ｝ｉに属する格子ａｊ毎に、該格子ａｊにおける値ｘｊの前記格子Ａｉにおける値Ｘｉへの寄与を重み付けｗｉｊとして所定の値に設定する工程と、
前記数値計算装置が、前記格子グループ｛ａ｝ｉに属する格子ａｊ毎に、該格子ａｊにおける値ｘｊ、重なり分率ｐｉｊ及び重み付けｗｉｊの積を計算する工程と、
前記数値計算装置が、前記積の前記格子グループ｛ａ｝ｉに属する全ての格子に亘る総和を格子Ａｉにおける値Ｘｉとして計算する工程と、
を有し、
前記精細化工程は、
前記数値計算装置が、前記格子ｂｊを被覆する粗いグリッド上の格子Ｂｉの集合を格子グループ｛Ｂ｝ｊとして計算する工程と、
前記数値計算装置が、前記格子グループ｛Ｂ｝ｊに属する格子Ｂｉ毎に、該格子Ｂｉと格子ｂｊとの共通部分の長さ、面積又は体積と、該格子Ｂｉの長さ、面積又は体積と、の比率を重なり分率ｑｊｉとして計算する工程と、
前記数値計算装置が、前記格子グループ｛Ｂ｝ｊに属する格子Ｂｉ毎に、該格子Ｂｉにおける値Ｙｉの前記格子ｂｊにおける値ｙｊへの寄与を重み付けｖｊｉとして所定の値に設定する工程と、
前記数値計算装置が、前記格子グループ｛Ｂ｝ｊに属する格子Ｂｉ毎に、該格子Ｂｉにおける値Ｙｉ、重なり分率ｑｊｉ及び重み付けｖｊｉの積を計算する工程と、
前記数値計算装置が、前記積の前記格子グループ｛Ｂ｝ｊに属する全ての格子に亘る総和を格子ｂｊにおける値ｙｊとして計算する工程と、
を有する
ことを特徴とするマルチグリッド法を用いた連立一次方程式の数値計算方法。
【請求項２】
前記数値計算装置が、前記粗視化工程において重み付けｗｉｊを格子の位置に依存しない一定値に設定し、
前記数値計算装置が、前記精細化工程において重み付けｖｊｉを格子の位置に依存しない一定値に設定する
ことを特徴とする請求項１に記載のマルチグリッド法を用いた連立一次方程式の数値計算方法。
【請求項３】
前記数値計算装置が、前記粗視化工程において重み付けｗｉｊを前記格子グループ｛ａ｝ｉに属する格子ａｊの中で前記格子Ａｉとの中心間距離が最短の格子以外は全てゼロに設定し、
前記数値計算装置が、前記精細化工程において重み付けｖｊｉを前記格子グループ｛Ｂ｝ｊに属する格子Ｂｉの中で前記格子ｂｊとの中心間距離が最短の格子以外は全てゼロに設定する
ことを特徴とする請求項１に記載のマルチグリッド法を用いた連立一次方程式の数値計算方法。
【請求項４】
前記数値計算装置は、細かいグリッドから粗いグリッドを生成するグリッド生成工程を実行するものであって、
前記グリッド生成工程は、
前記細かいグリッドの各座標軸方向の格子分割数Ｎｆが偶数の場合に、前記数値計算装置が該格子分割数Ｎｆを２で割った値を前記粗いグリッドの各座標軸方向の格子分割数Ｎｃとして決定し格子分割数比Ｎｆ÷Ｎｃが整数となる粗いグリッドを生成する工程と、
前記細かいグリッドの各座標軸方向の格子分割数Ｎｆが奇数の場合に、前記数値計算装置が該格子分割数Ｎｆに１を加算または減算してから２で割った値を粗いグリッドの各座標軸方向の格子分割数Ｎｃとして決定し格子分割数比Ｎｆ÷Ｎｃが整数ではない有理数となる粗いグリッドを生成する工程と、
を有する
ことを特徴とする請求項１から３の何れか１項に記載のマルチグリッド法を用いた連立一次方程式の数値計算方法。
【請求項５】
前記数値計算装置が、前記グリッド生成工程において粗いグリッドの格子分割数Ｎｃを（Ｎｆ＋１）÷２に決定した場合、前記数値計算装置は、
前記粗視化工程において前記格子グループ｛ａ｝ｉを、
（Ａ）ｉ＝０の場合、細かいグリッド上の格子番号がｉ×２及びｉ×２＋１の格子の集合、
（Ｂ）ｉ＝Ｎｃ−１の場合、細かいグリッド上の格子番号がｉ×２−１及びｉ×２の格子の集合、
（Ｃ）それ以外の場合、細かいグリッド上の格子番号がｉ×２−１、ｉ×２及びｉ×２＋１の格子の集合、
として計算し、
前記精細化工程において前記格子グループ｛Ｂ｝ｊを、
（Ａ）ｊが偶数の場合、粗いグリッド上の格子番号がｊ÷２の格子の集合、
（Ｂ）ｊが奇数の場合、粗いグリッド上の格子番号が（ｊ−１）÷２及び（ｊ−１）÷２＋１の格子の集合、
として計算し、
前記数値計算装置が、前記グリッド生成工程において粗いグリッドの格子分割数Ｎｃを（Ｎｆ−１）÷２に決定した場合、前記数値計算装置は、
前記粗視化工程において前記格子グループ｛ａ｝ｉを、細かいグリッド上の格子番号がｉ×２、ｉ×２＋１及びｉ×２＋２の格子の集合として計算し、
前記精細化工程において前記格子グループ｛Ｂ｝ｊを、
（Ａ）ｊ＝０の場合、粗いグリッド上の格子番号がｊ÷２の格子の集合、
（Ｂ）ｊ＝Ｎｆ−１の場合、粗いグリッド上の格子番号がｊ÷２−１の格子の集合、
（Ｃ）ｊが０またはＮｆ−１以外の偶数の場合、粗いグリッド上の格子番号がｊ÷２及びｊ÷２−１の格子の集合、
（Ｄ）ｊが奇数の場合、粗いグリッド上の格子番号が（ｊ−１）÷２の格子の集合、として計算する
ことを特徴とする請求項４に記載のマルチグリッド法を用いた連立一次方程式の数値計算方法。
【請求項６】
前記数値計算装置は、細かいグリッドから粗いグリッドを生成するグリッド生成工程を実行するものであって、
前記グリッド生成工程は、
前記細かいグリッドの各座標軸方向の格子分割数Ｎｆが偶数の場合に、前記数値計算装置が該格子分割数Ｎｆを２で割った値を前記粗いグリッドの各座標軸方向の格子分割数Ｎｃとして決定し格子分割数比Ｎｆ÷Ｎｃが整数となる粗いグリッドを生成する工程と、
前記細かいグリッドの各座標軸方向の格子分割数Ｎｆが奇数の場合に、半開区間［（Ｎｆ−１）÷２，Ｎｆ）に含まれる２のべき乗、２のべき乗と３の積、２のべき乗と５の積または８の倍数のいずれかを粗いグリッドの各座標軸方向の格子分割数Ｎｃとして決定し格子分割数比Ｎｆ÷Ｎｃが整数ではない有理数となる粗いグリッドを生成する工程と、を有する
ことを特徴とする請求項１から３の何れか１項に記載のマルチグリッド法を用いた連立一次方程式の数値計算方法。
【請求項７】
マルチグリッド法を用いて連立一次方程式を数値計算する数値計算装置であって、
格子分割数Ｎｆの細かいグリッド上のｊ番目の格子ａｊにおける値ｘｊから、格子分割数Ｎｃ（＜Ｎｆ）の粗いグリッド上のｉ番目の格子Ａｉにおける値Ｘｉへの変換を
【数３】
2011154510.tif 000009
により行う粗視化手段と、
格子分割数Ｎｃの粗いグリッド上のｉ番目の格子Ｂｉにおける値Ｙｉから、格子分割数Ｎｆ（＞Ｎｃ）の細かいグリッド上のｊ番目の格子ｂｊにおける値ｙｊへの変換を
【数４】
2011154510.tif 000010
により行う精細化手段と、
を備え、
前記粗視化手段は、
前記格子Ａｉを被覆する細かいグリッド上の格子ａｊの集合を格子グループ｛ａ｝ｉとして計算するする手段と、
前記格子グループ｛ａ｝ｉに属する格子ａｊ毎に、該格子ａｊと格子Ａｉとの共通部分の長さ、面積又は体積と、該格子ａｊの長さ、面積又は体積と、の比率を重なり分率ｐｉｊとして計算する手段と、
前記格子グループ｛ａ｝ｉに属する格子ａｊ毎に、該格子ａｊにおける値ｘｊの前記格子Ａｉにおける値Ｘｉへの寄与を重み付けｗｉｊとして所定の値に設定する手段と、
前記格子グループ｛ａ｝ｉに属する格子ａｊ毎に、該格子ａｊにおける値ｘｊ、重な
り分率ｐｉｊ及び重み付けｗｉｊの積を計算する手段と、
前記積の前記格子グループ｛ａ｝ｉに属する全ての格子に亘る総和を格子Ａｉにおける値Ｘｉとして計算する手段と、
を備え、
前記精細化手段は、
前記格子ｂｊを被覆する粗いグリッド上の格子Ｂｉの集合を格子グループ｛Ｂ｝ｊとして計算する手段と、
前記格子グループ｛Ｂ｝ｊに属する格子Ｂｉ毎に、該格子Ｂｉと格子ｂｊとの共通部分の長さ、面積又は体積と、該格子Ｂｉの長さ、面積又は体積と、の比率を重なり分率ｑｊｉとして計算する手段と、
前記格子グループ｛Ｂ｝ｊに属する格子Ｂｉ毎に、該格子Ｂｉにおける値Ｙｉの前記格子ｂｊにおける値ｙｊへの寄与を重み付けｖｊｉとして所定の値に設定する手段と、
前記格子グループ｛Ｂ｝ｊに属する格子Ｂｉ毎に、該格子Ｂｉにおける値Ｙｉ、重なり分率ｑｊｉ及び重み付けｖｊｉの積を計算する手段と、
前記積の前記格子グループ｛Ｂ｝ｊに属する全ての格子に亘る総和を格子ｂｊにおける値ｙｊとして計算する手段と、
を備えることを特徴とする数値計算装置。
【発明の詳細な説明】
【技術分野】
【０００１】
本発明は、連立一次方程式の数値計算方法及び装置に関し、特にマルチグリッド法を用いた連立一次方程式の数値解法に関する。
【背景技術】
【０００２】
一般に、ある物理系に対して物理現象の問題を解く場合は、系を離散格子からなるグリッド上に離散化する。そして、物理系を表す基礎方程式は、グリッド上の差分方程式になり、大規模な連立一次方程式Ａ×ｘ＝ｂを解く問題に帰着する。Ａは係数行列、ｂは定数項ベクトル、ｘは解ベクトルである。
このような大規模な連立一時方程式を解く数値計算手法としては、直接法であるガウス消去法や、反復法であるヤコビ法、ガウス−ザイデル法、ＳＯＲ法、共役勾配法などが知られている。特に連立一次方程式の係数行列が零成分の多い疎行列の場合には反復解法が有効である。
しかし、通常用いられるＳＯＲ法や共役勾配法は、主に隣接格子を介して誤差情報が伝播するので、誤差の短波長成分の減衰は速いが、誤差の長波長成分の減衰は遅い。そのため、結果的に十分満足できる解へ収束するまでの反復回数が多くなり、長い計算時間を要する。
【０００３】
そこで、誤差の長波長成分の減衰を速くして反復回数を少なくする方法として、マルチグルッド法（多重格子法、ＭＧ法）が提案されている。この方法は、与えられたグリッド上の連立一次方程式と、それをより粗いグリッド上に変数変換したものをそれぞれに解くことによって、誤差の短波長成分と長波長成分の両方を効率よく減衰させることができる（非特許文献１、非特許文献２を参照）。
図２はマルチグリッド法において、粗さの異なるグリッド間で変数を変換する処理について説明するための図である。解くべき連立一次方程式Ａ×ｘ＝ｂが与えられている最も細かいグリッドを第１グリッドとし、第１グリッドよりも粗いグリッドを順番に第２グリッド、第３グリッドと階層化する。細かいグリッド上の変数を粗いグリッド上に変換する操作が粗視化処理であり、リストリクション（ｒｅｓｔｒｉｃｔｉｏｎ）とも呼ばれている。粗いグリッド上の変数を細かいグリッド上に変換する操作が精細化処理であり、プロロンゲーション（ｐｒｏｌｏｎｇａｔｉｏｎ）とも呼ばれている。
【０００４】
マルチグリッド法の典型的な処理例を図３を用いて説明する。
ステップＳ３１０では、与えられた第１グリッド（細かいグリッド）の格子分割数を２で割った値が格子分割数となる第２グリッド（粗いグリッド）を生成する。格子分割数が２で割れない場合は処理を終了する。
ステップＳ３２０では、解くべき連立一次方程式の係数行列（細かいグリッド上の係数行列）を粗視化して、粗いグリッド上の係数行列を生成する。ここでは、細かいグリッドの格子分割数と粗いグリッドの格子分割数の比が２になる場合の粗視化作用素と精細化作用素を用いる。
細かいグリッド上の係数行列Ａｆから粗いグリッド上の係数行列Ａｃを生成する方法としては、ＤＣＡとＧＣＡがある。ＤＣＡ（Ｄｉｒｅｃｔ  Ｃｏａｒｓｅ  ｇｒｉｄ  Ａｐｐｒｏｘｉｍａｔｉｏｎ）は、物理系を表す基礎方程式を粗いグリッド上で直接離散化して連立一次方程式を導出する方法である。ＧＣＡ（Ｇａｌｅｒｋｉｎ  Ｃｏａｒｓｅ  ｇｒｉｄ  Ａｐｐｒｏｘｉｍａｔｉｏｎ）は、細かいグリッド上の係数行列Ａｆを粗視化作用素の行列形式Ｒと精細化作用素の行列形式Ｐで挟んで行列積Ｒ×Ａｆ×Ｐを計算して導出する方法である。細かいグリッドの分割数をＮｆ、粗いグリッドの分割数をＮｃとすると、ＡｆはＮｆ行Ｎｆ列の正方行列、ＲはＮｃ行Ｎｆ列の行列、ＰはＮｆ行Ｎｃ列の行列
であり、Ｒ×Ａｆ×Ｐ＝ＡｃはＮｃ行Ｎｃ列の正方行列になる。
【０００５】
ステップＳ３３０では、解くべき連立一次方程式Ａ×ｘ＝ｂの解ベクトルｘと定数項ベクトルｂについて、細かいグリッドと粗いグリッドとの格子分割数比が２になる場合の粗視化作用素を作用させて粗いグリッド上の連立一次方程式を準備する。
ステップＳ３４０では、粗いグリッド上の連立一次方程式を解くことによって、誤差の長波長成分を減衰させ易くする。ここで用いる求解法はガウス−ザイデル法やＳＯＲ法等、適当に選んで良い。
ステップＳ３５０では、粗いグリッド上で解いた連立一次方程式の解ベクトルについて、細かいグリッドの格子分割数と粗いグリッドの格子分割数の比が２になる場合の精細化作用素を作用させて、元の第１グリッド上の変数に戻す。
ステップＳ３６０では、初めに与えられた連立一次方程式を解くことによって、誤差の短波長成分を減衰させ易くする。このときステップＳ３５０で精細化された解ベクトルを初期値として用いることによって、誤差の長波長成分と短波長成分を同時に減衰させ易くすることができる。
図３の処理例は第１グリッドと第２グリッドのみ計算しているが、実用上はＳ３３０からＳ３６０の処理を再帰的に実行してより粗いグリッド上でも計算することによって、誤差の様々な波長成分を同時に減衰させ易くしている。
【０００６】
Ｓ３２０、Ｓ３３０、Ｓ３５０で使用する粗視化処理と精細化処理は、細かいグリッドの格子分割数と粗いグリッドの格子分割数の比が２になる場合が最も簡単であるため、格子分割数は２のべき乗になることが多い。また、特許文献２のように格子分割数が素因数分解可能な整数の積でもマルチグリッド法は適用可能である。この場合、Ｓ３１０では設定値以下の素数で格子分割数を割った数を粗いグリッドの格子分割数に設定し、Ｓ３２０、Ｓ３３０、Ｓ３５０で使用する粗視化処理と精細化処理は、当該素数に対応したものを使用する。
【先行技術文献】
【特許文献】
【０００７】
【特許文献１】特開２００６−２９３４８８号公報
【特許文献２】特開２００７−２８６７３８号公報
【非特許文献】
【０００８】
【非特許文献１】W.H.Press, S.A.Teukolsky, W.T.Vetterling, B.P.Flannery："NUMERICAL RECIPES in Fortran 77 Second Edition", Cambridge University Press, 1992, p.862
【非特許文献２】W.L.Briggs, V.E.Henson, S.F.McCormick："A Multigrid Tutorial Second Edition", SIAM, 2000.
【発明の概要】
【発明が解決しようとする課題】
【０００９】
しかし、上述した従来のマルチグリッド法の処理は、元の第１グリッドの格子分割数が２のべき乗または一般には素因数分解可能な整数の積という条件を満たす問題にしか適用できない。
解くべき問題の格子分割数がこの条件を満たさなければ、連立一次方程式の求解法として従来のマルチグリッド法が適用できないため、計算時間の長くかかる別の求解法を選択する必要がある。或いは従来のマルチグリッド法を適用するために前記条件を満たすよう格子分割数を変更する必要がある。そのため、余分な格子を追加することによって計算時間が増大したり、格子間隔が必要以上に狭まって物体表面と格子境界が揃わなかったりする場合がある。
また、解析空間の格子分割数が素数、例えば９７の場合、９７より大きい２のべき乗である１２８を用いて、元のグリッド上の物理量を格子分割数１２８のグリッド上の物理量に変換してからマルチグリッド法を用いて求解する。その後、格子分割数１２８のグリッド上から元の格子分割数９７のグリッド上に物理量を補間して戻すという処理を行って求解する方法も可能だが、補間に伴って誤差が生じる場合がある。
【００１０】
本発明はこのような背景のもとになされたものであり、任意の格子分割数の問題に対して適用可能なマルチグリッド法を用いた連立一次方程式の数値計算方法及び装置を提供することを目的とする。
【課題を解決するための手段】
【００１１】
本発明は、数値計算装置におけるマルチグリッド法を用いた連立一次方程式の数値計算方法であって、
前記数値計算装置が、
格子分割数Ｎｆの細かいグリッド上のｊ番目の格子ａｊにおける値ｘｊから、格子分割数Ｎｃ（＜Ｎｆ）の粗いグリッド上のｉ番目の格子Ａｉにおける値Ｘｉへの変換を
【数１】
2011154510.tif 000003
により行う粗視化工程と、
格子分割数Ｎｃの粗いグリッド上のｉ番目の格子Ｂｉにおける値Ｙｉから、格子分割数Ｎｆ（＞Ｎｃ）の細かいグリッド上のｊ番目の格子ｂｊにおける値ｙｊへの変換を
【数２】
2011154510.tif 000004
により行う精細化工程と、
を実行し、
前記粗視化工程は、
前記数値計算装置が、前記格子Ａｉを被覆する細かいグリッド上の格子ａｊの集合を格子グループ｛ａ｝ｉとして計算するする工程と、
前記数値計算装置が、前記格子グループ｛ａ｝ｉに属する格子ａｊ毎に、該格子ａｊと格子Ａｉとの共通部分の長さ、面積又は体積と、該格子ａｊの長さ、面積又は体積と、の比率を重なり分率ｐｉｊとして計算する工程と、
前記数値計算装置が、前記格子グループ｛ａ｝ｉに属する格子ａｊ毎に、該格子ａｊにおける値ｘｊの前記格子Ａｉにおける値Ｘｉへの寄与を重み付けｗｉｊとして所定の値に設定する工程と、
前記数値計算装置が、前記格子グループ｛ａ｝ｉに属する格子ａｊ毎に、該格子ａｊにおける値ｘｊ、重なり分率ｐｉｊ及び重み付けｗｉｊの積を計算する工程と、
前記数値計算装置が、前記積の前記格子グループ｛ａ｝ｉに属する全ての格子に亘る総和を格子Ａｉにおける値Ｘｉとして計算する工程と、
を有し、
前記精細化工程は、
前記数値計算装置が、前記格子ｂｊを被覆する粗いグリッド上の格子Ｂｉの集合を格子グループ｛Ｂ｝ｊとして計算する工程と、
前記数値計算装置が、前記格子グループ｛Ｂ｝ｊに属する格子Ｂｉ毎に、該格子Ｂｉと格子ｂｊとの共通部分の長さ、面積又は体積と、該格子Ｂｉの長さ、面積又は体積と、の比率を重なり分率ｑｊｉとして計算する工程と、
前記数値計算装置が、前記格子グループ｛Ｂ｝ｊに属する格子Ｂｉ毎に、該格子Ｂｉにおける値Ｙｉの前記格子ｂｊにおける値ｙｊへの寄与を重み付けｖｊｉとして所定の値に設定する工程と、
前記数値計算装置が、前記格子グループ｛Ｂ｝ｊに属する格子Ｂｉ毎に、該格子Ｂｉにおける値Ｙｉ、重なり分率ｑｊｉ及び重み付けｖｊｉの積を計算する工程と、
前記数値計算装置が、前記積の前記格子グループ｛Ｂ｝ｊに属する全ての格子に亘る総和を格子ｂｊにおける値ｙｊとして計算する工程と、
を有する
ことを特徴とするマルチグリッド法を用いた連立一次方程式の数値計算方法である。
【発明の効果】
【００１２】
本発明によれば、任意の格子分割数の問題に対して適用可能なマルチグリッド法を用いた連立一次方程式の数値計算方法及び装置を提供することができる。
【図面の簡単な説明】
【００１３】
【図１】実施例に係る数値計算手法のフローチャート
【図２】マルチグリッド法における粗視化及び精細化を説明する図
【図３】従来のマルチグリッド法の計算処理例を示すフローチャート
【図４】格子分割数比が２の場合の粗視化及び精細化作用を示す概念図
【図５】格子分割数比が７／４の場合の粗視化及び精細化作用を示す概念図
【図６】実施例の数値計算装置の基本構成を示す図
【発明を実施するための形態】
【００１４】
図６は、本実施形態に係る数値計算装置として機能するコンピュータの基本構成を示すブロック図である。該コンピュータは、大略、バス２０７を介して接続されたＣＰＵ２０１、表示装置２０２，入力装置２０３、ＲＡＭ２０４、ハードディスク装置２０５、ＮＩＣ２０６等から構成される。
ＣＰＵ（中央演算装置）２０１で、ＲＡＭ（主記憶装置）２０４に記憶されているプログラムやデータを用いて本装置全体の制御を行うと共に、後述する各処理を実行する。表示装置２０２は、液晶画面などにより構成されており、ＣＰＵ２０１による処理結果を文字又は画像などで表示することができる。入力装置２０３は、キーボードやマウスなどの操作系の装置により構成されており、各種の指示をＣＰＵ２０１に対して入力することができる。ＲＡＭ２０４は、ハードディスク装置２０５からロードされたプログラムやデータを一時的に記憶する為のエリアを備えると共に、ＣＰＵ２０１が各種の処理を実行する際に必要なワークエリアを備える。ハードディスク装置２０５には、ＯＳ（オペレーティングシステム）や、コンピュータが行うべき後述の各処理をＣＰＵ２０１に実行させるためのプログラムやデータが保存されている。これらのプログラムはＲＡＭ２０４にロードされて、ＣＰＵ２０１によって実行されることで、本コンピュータが連立一次方程式を数値的に解くための数値計算装置として機能する。ＮＩＣ（ネットワークインターフェース）２０６は、外部装置とのデータ通信を行う際のインターフェースとして機能するものであり、本コンピュータはこのＮＩＣ２０６を介して外部装置とプログラムやデータの送受信を行う。なお、ハードディスク装置２０５が保存しているプログラムやデータの代わりに、ＮＩＣ２０６を介して外部装置から受信したプログラムやデータをＣＰＵ２０１による処理対象としても良い。
【００１５】
次に、上記構成を備えるコンピュータが行う、連立一次方程式の解を求める為の数値解法について説明する。図１は、本実施形態に係るコンピュータが行う、連立一次方程式の解を求めるための数値計算処理のフローチャートである。なお、ここでは係数行列Ａｆ、解ベクトルｘ、定数項ベクトルｂとして、Ａｆ×ｘ＝ｂで表される連立一次方程式の求解処理に、本発明のマルチグリッド法を適用する。
ここでは、１次元空間での物理現象シミュレーションを行うために基礎方程式を格子中心（セル中心）に値が定義された等間隔のグリッド上で３点差分により離散化し、係数行列が３重対角となっている連立一次方程式を解くことを想定している。
まず、ステップＳ２１０ではグリッドの格子分割数Ｎｆについて偶数か奇数かを判定し、偶数の場合は２で割った値を、奇数の場合は所定の方法で算出した値を１段粗いグリッドの格子分割数Ｎｃとして粗いグリッドを生成する。
【００１６】
なお、３次元空間の場合は、ｘ軸，ｙ軸，ｚ軸方向の格子分割数Ｎｘ，Ｎｙ，Ｎｚについてそれぞれ偶奇を判定して粗いグリッドの各座標軸方向の格子分割数を計算して粗いグリッド生成すればよい。
Ｎｆが奇数の場合のＮｃの算出方法としては、Ｎｆに１を加算してから２で割った値をＮｃとするＮｃ＝（Ｎｆ＋１）÷２で表される方法や、逆にＮｆから１を減算してから２で割った値をＮｃとするＮｃ＝（Ｎｆ−１）÷２で表される方法がある。さらに半開区間［（Ｎｆ−１）÷２，Ｎｆ）に含まれる２のべき乗、２のべき乗と３の積、２のべき乗と５の積又は８の倍数のいずれかを選択することによってＮｃを決定する方法がある。
例えばＮｆが奇数の場合にＮｃ＝（Ｎｆ＋１）÷２とする方法では、Ｎｆ＝７ならばＮｃ＝（７＋１）÷２＝４、Ｎｆ＝１７ならばＮｃ＝（１７＋１）÷２＝９となる。
ステップＳ２２０では、ＤＣＡまたはＧＣＡによって解くべき連立一次方程式の係数行列Ａｆを粗視化して、粗いグリッド上の係数行列Ａｃを生成する。この際使用する粗視化作用素と精細化作用素は、ステップＳ２３０またはステップＳ２３１とステップＳ２５０またはステップＳ２５１のものと同等である。
【００１７】
ステップＳ２３０及びステップＳ２３１では、解くべき連立一次方程式Ａｆ×ｘ＝ｂの解ベクトルｘと定数項ベクトルｂに粗視化作用素を作用させて粗いグリッド上の連立一次方程式を準備する。細かいグリッドの格子分割数Ｎｆと粗いグリッドの格子分割数Ｎｃの比Ｎｆ÷Ｎｃが整数の場合、ステップＳ２３０の粗視化作用素を作用させる。細かいグリッドの格子分割数Ｎｆと粗いグリッドの格子分割数Ｎｃの比Ｎｆ÷Ｎｃが非整数有理数の場合、ステップＳ２３１の粗視化作用素を作用させる。
本実施例では、細かいグリッドの格子分割数Ｎｆが偶数の場合のみ格子分割数比Ｎｆ÷Ｎｃが整数（２）となり、ステップＳ２３０が実行される。細かいグリッドの格子分割数Ｎｆが奇数の場合に格子分割数比Ｎｆ÷Ｎｃが非整数有理数となり、ステップＳ２３１が実行される。
【００１８】
図４（Ａ）は、ステップＳ２３０で使用する格子分割数比Ｎｆ÷Ｎｃが整数の場合の粗視化作用素、特に、細かいグリッドの格子分割数Ｎｆが偶数の場合の粗視化作用素の粗視化作用を示す概念図である。図中の矢印に付随する数字は第１グリッド（細かいグリッド）上の値を第２グリッド（粗いグリッド）上の値に変換する際の分配比率を示している。
ステップＳ２３１で使用する格子分割数比Ｎｆ÷Ｎｃが非整数有理数の場合の粗視化作用素は、格子間の被覆状態を表す長さ分率ｐｉｊを用いて式（１）で定義する。
【数３】
2011154510.tif 000005
ここで、Ｘｉは粗いグリッド上のｉ番目の格子Ａｉにおける値、ｘｊは細かいグリッド上のｊ番目の格子ａｊでの値、｛ａ｝ｉは格子Ａｉを被覆する細かいグリッド上の格子の集合である。
【００１９】
ｐｉｊは格子Ａｉを被覆する格子グループ｛ａ｝ｉに属する各格子ａｊと格子Ａｉとの重なり具合を表す重なり分率である。重なり分率ｐｉｊは、格子グループ｛ａ｝ｉに属する格子ａｊ毎に、該格子ａｊと格子Ａｉとの共通部分の長さ、面積又は体積と、該格子ａ
ｊの長さ、面積又は体積と、の比率を計算して求められる。重なり分率ｐｉｊは、細かいグリッドと粗いグリッドの幾何学的関係により各格子について一意に定まる。一次元の場合は、格子ａｊの長さと、該格子ａｊのうち格子Ａｉに含まれる部分の長さと、の比率が、格子Ａｉにおける格子ａｊの重なり分率ｐｉｊであり、これを長さ分率と称する。同様に、二次元の場合は、細かいグリッドの格子ａｊの面積と、該格子ａｊのうち粗いグリッドの格子Ａｉに含まれる部分の面積と、の比率が重なり分率であり、これを面積分率と称する。三次元の場合は、細かいグリッドの格子ａｊの体積と、該格子ａｊのうち粗いグリッドの格子Ａｉに含まれる部分の体積と、の比率が重なり分率であり、これを体積分率と称する。
【００２０】
ｗｉｊは前記格子グループ｛ａ｝ｉに属する各格子ａｊにおける値ｘｊの、格子Ａｉの値Ｘｉへの寄与を表す重み付けであり、任意の値に設定できる。重み付けｗｉｊの設定方法としては、格子の位置に依存しない一定値Ｎｃ／Ｎｆとする方法や、格子Ａｉと格子ａｊの中心間距離の関数とする方法がある。また、格子Ａｉと格子ａｊの中心間距離が最短の格子以外は全てゼロに設定する方法、つまり格子Ａｉに最も近い格子ａｊ０における値のみを代表値として採用する方法がある。格子分割数比が非整数有理数の場合、上記の重なり分率ｐｉｊ、重み付けｗｉｊ及び格子グループに属する各格子ａｊにおける値ｘｊの積の、格子グループ｛ａ｝ｉに属する全ての格子に亘る総和を、粗いグリッドの格子Ａｉにおける値Ｘｉとする。
【００２１】
図５（Ａ）は、ステップＳ２３１で使用する格子分割数比Ｎｆ÷Ｎｃが非整数有理数（７／４）の場合の粗視化作用素の粗視化作用を示す概念図である。図中の矢印に付随する数字は長さ分率に基づいて第１グリッド（細かいグリッド）上の値を第２グリッド（粗いグリッド）上の値に変換する際の分配比率を示している。式１に示すように、細かいグリッドの格子ａｊにおける値ｘｊが粗いグリッドの格子Ａｉにおける値Ｘｉに分配される比率は、格子Ａｉにおける格子ａｊの長さ分率ｐｉｊと重み付けｗｉｊとの積である。図５（Ａ）に示した分配比率は、重み付けｗｉｊを一定値Ｎｃ／Ｎｆ＝４／７に設定した場合の値である。
【００２２】
例えば、粗視化作用により粗い第２グリッド上の３番目の格子Ａ２の値Ｘ２を求める方法は次のようになる。格子Ａ２を被覆する細かい第１グリッド上の格子はａ３，ａ４，ａ５の３個である。すなわち、｛ａ｝２＝｛ａ３，ａ４，ａ５｝である。格子ａ３はその１／２の部分が格子Ａ２に含まれるので、格子Ａ２における格子ａ３の長さ分率ｐ２３は１／２である。格子ａ４はその全部が格子Ａ２に含まれるので、格子Ａ２における格子ａ４の長さ分率ｐ２４は１である。格子ａ５はその１／４の部分が格子Ａ２に含まれるので、格子Ａ２における格子ａ５の長さ分率ｐ２５は１／４である。ここでは、重み付けｗｉｊは格子分割数比４／７で一定値とする。従って第２グリッド上の格子Ａ２における値Ｘ２は、Ｘ２＝ｘ３ｐ２３ｗ２３＋ｘ４ｐ２４ｗ２４＋ｘ５ｐ２５ｗ２５＝ｘ３×（１／２）×（４／７）＋ｘ４×１×（４／７）＋ｘ５×（１／４）×（４／７）＝ｘ３×２／７＋ｘ４×４／７＋ｘ５×１／７となる。
【００２３】
ここで、一般に粗いグリッド上のｉ番目の格子Ａｉを被覆する細かいグリッド上の格子のグループ｛ａ｝ｉに属する格子の格子番号を求めるには、格子位置を表す実数型の座標値から整数型に変換する処理または有理数型の計算処理が必要である。
しかし、Ｎｃ＝（Ｎｆ＋１）÷２の場合は、前記格子グループ｛ａ｝ｉに属する格子の格子番号は、ｉ＝０の場合はｉ×２とｉ×２＋１、すなわち｛ａ｝０＝｛ａ０，ａ１｝となる。また、ｉ＝Ｎｃ−１の場合はｉ×２−１とｉ×２、すなわち｛ａ｝Ｎｃ−１＝｛ａ２Ｎｃ−３，ａ２Ｎｃ−２｝となる。また、他の場合はｉ×２−１、ｉ×２及びｉ×２＋１、すなわち｛ａ｝ｉ＝｛ａ２ｉ−１，ａ２ｉ，ａ２ｉ＋１｝となり、簡単な整数型の計算のみで算出できる。また、Ｎｃ＝（Ｎｆ−１）÷２の場合は、前記格子グループ｛ａ｝
ｉに属する格子の格子番号はｉ×２、ｉ×２＋１及びｉ×２＋２、すなわち｛ａ｝ｉ＝｛ａ２ｉ，ａ２ｉ＋１，ａ２ｉ＋２｝となり、簡単な整数型の計算のみで算出できる。したがって、Ｎｃ＝（Ｎｆ±１）÷２の場合には高速に格子グループ｛ａ｝ｉに属する格子の格子番号を求めることが可能である。
【００２４】
また、ステップＳ２３１で使用する粗視化作用素の行列形式Ｒの１行当たりの非零要素数を数えると、Ｎｃ≧（Ｎｆ−１）÷２では３個以下だが、Ｎｃ＜（Ｎｆ−１）÷２では４個以上になる場合が生じてしまう。つまりＮｃ＝（Ｎｆ±１）÷２の場合は、粗視化作用素を作用させる時の演算量の増大を抑制できる。
ステップＳ２４０では、粗いグリッド上の連立一次方程式を解くことによって、誤差の長波長成分を減衰させ易くする。ここで用いる求解法はガウス−ザイデル法やＳＯＲ法等、適当に選んで良い。
【００２５】
ステップＳ２５０およびステップＳ２５１では、粗いグリッド上で解いた連立一次方程式の解ベクトルに精細化作用素を作用させて元の第１グリッド上の変数に戻す。細かいグリッドの格子分割数Ｎｆと粗いグリッドの格子分割数Ｎｃの比Ｎｆ÷Ｎｃが整数の場合、ステップＳ２５０の精細化作用素を作用させる。細かいグリッドの格子分割数Ｎｆと粗いグリッドの格子分割数Ｎｃの比Ｎｆ÷Ｎｃが非整数有理数の場合、ステップＳ２５１の精細化作用素を作用させる。本実施例では、細かいグリッドの格子分割数Ｎｆが偶数の場合のみ格子分割数比Ｎｆ÷Ｎｃが整数（２）となり、ステップＳ２５０が実行される。細かいグリッドの格子分割数Ｎｆが奇数の場合に格子分割数比Ｎｆ÷Ｎｃが非整数有理数となり、ステップＳ２５１が実行される。
【００２６】
図４（Ｂ）は、ステップＳ２５０で使用する格子分割数比Ｎｆ÷Ｎｃが整数の場合の精細化作用素、特に、細かいグリッドの格子分割数Ｎｆが偶数の場合の精細化作用素の精細化作用を示す概念図である。図中の矢印に付随する数字は第２グリッド（粗いグリッド）上の値を第１グリッド（細かいグリッド）上の値に変換する際の分配比率を示している。
Ｓ２５１で使用する格子分割数比Ｎｆ÷Ｎｃが非整数有理数の場合の精細化作用素は、格子間の被覆状態を表す長さ分率ｑｊｉを用いて式（２）で定義する。
【数４】
2011154510.tif 000006
【００２７】
ここで、ｙｊは細かいグリッド上のｊ番目の格子ｂｊにおける値、Ｙｉは粗いグリッド上のｉ番目の格子Ｂｉでの値、｛Ｂ｝ｊは格子ｂｊを被覆する粗いグリッド上の格子の集合である。
【００２８】
ｑｊｉは格子ｂｊを被覆する格子グループ｛Ｂ｝ｊに属する各格子Ｂｉと格子ｂｊとの重なり具合を表す重なり分率である。重なり分率ｑｊｉは、格子グループ｛Ｂ｝ｊに属する格子Ｂｉ毎に、該格子Ｂｉと格子ｂｊとの共通部分の長さ、面積又は体積と、該格子Ｂｉの長さ、面積又は体積と、の比率を計算して求められる。重なり分率ｑｊｉは、粗いグリッドと細かいグリッドとの幾何学的関係により各格子について一意に定まる。一次元の場合は、格子Ｂｉの長さと、該格子Ｂｉのうち格子ｂｊに含まれる部分の長さと、の比率が、格子ｂｊにおける格子Ｂｉの重なり分率ｑｊｉであり、これを長さ分率と称する。同様に、二次元の場合は、粗いグリッドの格子Ｂｉの面積と、該格子Ｂｉのうち細かいグリッドの格子ｂｊに含まれる部分の面積と、の比率が重なり分率であり、これを面積分率と称する。三次元の場合は、粗いグリッドの格子Ｂｉの体積と、該格子Ｂｉのうち細かいグリッドの格子ｂｊに含まれる部分の体積と、の比率が重なり分率であり、これを体積分率
と称する。
【００２９】
ｖｊｉは前記格子グループ｛Ｂ｝ｊに属する各格子Ｂｉにおける値Ｙｉの、格子ｂｊの値ｙｊへの寄与を表す重み付けであり、任意の値に設定できる。重み付けｖｊｉの設定方法としては、格子の位置に依存しない一定値Ｎｆ／Ｎｃとする方法や、格子ｂｊと格子Ｂｉの中心間距離の関数とする方法がある。また、格子ｂｊと格子Ｂｉの中心間距離が最短の格子以外は全てゼロに設定する方法、つまり格子ｂｊに最も近い格子Ｂｉ０における値のみを代表値として採用する方法がある。格子分割数比が非整数有理数の場合、上記の重なり分率ｑｊｉ、重み付けｖｊｉ及び格子グループに属する各格子Ｂｉにおける値Ｙｉの積の、格子グループ｛Ｂ｝ｊに属する全ての格子に亘る総和を、粗いグリッドの格子ｂｊにおける値ｙｊとする。
【００３０】
図５（Ｂ）は、ステップＳ２５１で使用する格子分割数比Ｎｆ÷Ｎｃが非整数有理数（７／４）の場合の精細化作用素の精細化作用を示す概念図である。図中の矢印に付随する数字は長さ分率に基づいて第２グリッド（粗いグリッド）上の値を第１グリッド（細かいグリッド）上の値に変換する際の分配比率を示している。式２に示すように、粗いグリッドの格子Ｂｉにおける値Ｙｉが細かいグリッドの格子ｂｊにおける値ｙｊに分配される比率は、格子ｂｊにおける格子Ｂｉの長さ分率ｑｊｉと重み付けｖｊｉとの積である。図５（Ｂ）に示した分配比率は、重み付けｖｊｉを一定値Ｎｆ／Ｎｃ＝７／４に設定した場合の値である。
例えば、精細化作用により細かい第１グリッド上の６番目の格子ｂ５の値ｙ５を求める方法は次のようになる。格子ｂ５を被覆する粗い第２グリッド上の格子はＢ２，Ｂ３の２個である。すなわち｛Ｂ｝５＝｛Ｂ２、Ｂ３｝である。格子Ｂ２はその１／７の部分が格子ｂ５に含まれるので、格子ｂ５における格子Ｂ２の長さ分率ｑ５２＝１／７である。格子Ｂ３はその３／７の部分が格子ｂ５に含まれるので、格子ｂ５における格子Ｂ３の長さ分率ｑ５３＝３／７である。ここでは、重み付けｖｊｉは格子分割数比７／４で一定値とする。従って第１グリッド上の格子ｂ５における値ｙ５は、ｙ５＝Ｙ２ｑ５２ｖ５２＋Ｙ３ｑ５３ｖ５３＝Ｙ２×（１／７）×（７／４）＋Ｙ３×（３／７）×（７／４）＝Ｙ２×１／４＋Ｙ３×３／４となる。
【００３１】
ここで、一般に細かいグリッド上のｊ番目の格子ｂｊを被覆する粗いグリッド上の格子のグループ｛Ｂ｝ｊに属する格子の格子番号を求めるには、格子位置を表す実数型の座標値から整数型に変換する処理または有理数型の計算処理が必要である。
しかし、Ｎｃ＝（Ｎｆ＋１）÷２の場合は、前記格子グループ｛Ｂ｝ｊに属する格子の格子番号は、ｊが偶数の場合はｊ÷２、すなわち｛Ｂ｝ｊ＝｛Ｂｊ／２｝となる。また、ｊが奇数の場合は（ｊ−１）÷２及び（ｊ−１）÷２＋１、すなわち｛Ｂ｝ｊ＝｛Ｂ（ｊ−１）／２、Ｂ（ｊ−１）／２＋１｝となり、簡単な整数型の計算のみで算出できる。また、Ｎｃ＝（Ｎｆ−１）÷２の場合は、前記格子グループ｛Ｂ｝ｊに属する格子の格子番号は、ｊ＝０の場合はｊ÷２、すなわち｛Ｂ｝０＝｛Ｂ０｝となる。また、ｊ＝Ｎｆ−１の場合はｊ÷２−１、すなわち｛Ｂ｝Ｎｆ−１＝｛Ｂ（Ｎｆ−１）／２−１｝となる。また、ｊが０又はＮｆ−１以外の偶数の場合はｊ÷２とｊ÷２−１、すなわち｛Ｂ｝ｊ＝｛Ｂｊ／２、Ｂｊ／２−１｝となる。また、ｊが奇数の場合は（ｊ−１）÷２、すなわち｛Ｂ｝ｊ＝｛Ｂ（ｊ−１）／２｝となり、簡単な整数型の計算のみで算出できる。したがって、Ｎｃ＝（Ｎｆ±１）÷２の場合には高速に格子グループ｛Ｂ｝ｊに属する格子の格子番号を求めることが可能である。
【００３２】
また、ステップＳ２５１で使用する精細化作用素の行列形式Ｐの１列当たりの非零要素数を数えると、Ｎｃ≧（Ｎｆ−１）÷２では最大３個以下だが、Ｎｃ＜（Ｎｆ−１）÷２では４個以上になる場合が生じてしまう。さらに粗いグリッド上の係数行列ＡｃをＡｃ＝Ｒ×Ａｆ×Ｐとして計算するＧＣＡの場合、Ｎｃ＞（Ｎｆ＋１）÷２では細かいグリッド
上の係数行列Ａｆが３重対角行列でもＡｃの１行あたりの非零要素数が４個以上になる場合が生じてしまう。つまりＮｃ＝（Ｎｆ±１）÷２の場合は、ＧＣＡにおいて係数行列の３重対角性が保たれるので、計算の際の演算量の増大を抑制できる。
ステップＳ２６０では、初めに与えられた連立一次方程式を解くことによって、誤差の短波長成分を減衰させ易くする。このときステップＳ２５０またはステップＳ２５１で精細化された解ベクトルを初期値として用いることによって、誤差の長波長成分と短波長成分を同時に減衰させ易くすることができる。
図１のフローでは第１グリッドと第２グリッドのみ計算しているが、実用上はＳ２３０またはＳ２３１からＳ２６０の処理を再帰的に実行してより粗いグリッド上でも計算することによって、誤差の様々な波長成分を同時に減衰させ易くしている。
【００３３】
以上説明したように、本実施例に係るマルチグリッド法による連立一次方程式の数値解法によれば、格子分割数が２のべき乗ではない任意の格子に対して、マルチグリッド法により連立一次方程式を解くことが可能になる。
図１で示した数値計算方法は、共役勾配法、共役残差法を含む非定常反復解法の前処理として利用することも可能である。また、上述した連立一次方程式の数値計算方法は、格子分割数が２のべき乗でない２次元問題や３次元問題にも適用可能である。
【００３４】
本発明は、上述した連立一次方程式の数値計算方法を実現するソフトウェアのプログラムコードを記憶した記憶媒体であっても良い。該記憶媒体をシステム或いは装置に供給し、そのシステム或いは装置のコンピュータ（またはＣＰＵやＭＰＵ等）が記憶媒体に格納されたプログラムコードを読み出し実行することによっても本発明の目的は達成される。
この場合、記憶媒体から読み出されたプログラムコード自体が前述した実施の形態の機能を実現することになり、そのプログラムコード及び該プログラムコードを記憶した記憶媒体は本発明を構成することになる。
また、プログラムコードを供給するための記憶媒体としては、例えば、フロッピー（登録商標）ディスク、ハードディスク、光磁気ディスク、ＣＤやＤＶＤ等の光ディスク、磁気テープ、不揮発性のメモリカード、ＲＯＭ等を用いることができる。または、プログラムコードをネットワークを介してダウンロードしてもよい。
【００３５】
また、コンピュータが読み出したプログラムコードの指示に基づき、コンピュータ上で稼動しているＯＳ（オペレーティングシステム）等が実際の処理の一部または全部を行い、その処理によって上述の実施例の機能が実現される実施形態でも良い。
さらに、記憶媒体から読み出されたプログラムコードが、コンピュータに挿入された機能拡張ボードやコンピュータに接続された機能拡張ユニットに備わるメモリに書き込まれるようにしても良い。その後、そのプログラムコードの指示に基づき、その拡張機能を拡張ボードや拡張ユニットに備わるＣＰＵ等が実際の処理の一部または全部を行い、その処理によって上述の実施例の機能が実現される実施形態でも良い。
【符号の説明】
【００３６】
２０１…ＣＰＵ，２０２…表示装置，２０３…入力装置，２０４…ＲＡＭ
【図１】
2011154510.tif 000011
【図２】
2011154510.tif 000012
【図３】
2011154510.tif 000013
【図４】
2011154510.tif 000014
【図５】
2011154510.tif 000015
【図６】
2011154510.tif 000016