JP2011154510A - 連立一次方程式の数値計算方法及び装置 - Google Patents

Abstract

【課題】任意の格子分割数の問題に対して適用可能なマルチグリッド法を用いた連立一次方程式の数値計算方法及び装置を提供する。
【解決手段】細かいグリッド上の値を粗いグリッド上の値へ変換する粗視化処理と、粗いグリッド上の値を細かいグリッド上の値へ変換する精細化処理と、において、変換先のグリッド上の格子を被覆する変換元のグリッド上の格子のグループを求め、該格子グループに属する各格子の大きさと、該格子と変換先の格子との共通部分の大きさと、の比率を重なり分率として求め、該格子グループに属する各格子における値の、変換先の格子における値への寄与を表す重み付けを設定し、該格子グループに属する全ての格子に亘って、各格子における値、重なり分率及び重み付けをの積の総和を求め、これを変換先の格子における値とする。細かいグリッドと粗いグリッドとの格子分割数比が非整数有理数の場合も粗視化及び精細化の処理が可能になる。
【選択図】図１

Description

本発明は、連立一次方程式の数値計算方法及び装置に関し、特にマルチグリッド法を用いた連立一次方程式の数値解法に関する。

一般に、ある物理系に対して物理現象の問題を解く場合は、系を離散格子からなるグリッド上に離散化する。そして、物理系を表す基礎方程式は、グリッド上の差分方程式になり、大規模な連立一次方程式Ａ×ｘ＝ｂを解く問題に帰着する。Ａは係数行列、ｂは定数項ベクトル、ｘは解ベクトルである。
このような大規模な連立一時方程式を解く数値計算手法としては、直接法であるガウス消去法や、反復法であるヤコビ法、ガウス−ザイデル法、ＳＯＲ法、共役勾配法などが知られている。特に連立一次方程式の係数行列が零成分の多い疎行列の場合には反復解法が有効である。
しかし、通常用いられるＳＯＲ法や共役勾配法は、主に隣接格子を介して誤差情報が伝播するので、誤差の短波長成分の減衰は速いが、誤差の長波長成分の減衰は遅い。そのため、結果的に十分満足できる解へ収束するまでの反復回数が多くなり、長い計算時間を要する。

そこで、誤差の長波長成分の減衰を速くして反復回数を少なくする方法として、マルチグルッド法（多重格子法、ＭＧ法）が提案されている。この方法は、与えられたグリッド上の連立一次方程式と、それをより粗いグリッド上に変数変換したものをそれぞれに解くことによって、誤差の短波長成分と長波長成分の両方を効率よく減衰させることができる（非特許文献１、非特許文献２を参照）。
図２はマルチグリッド法において、粗さの異なるグリッド間で変数を変換する処理について説明するための図である。解くべき連立一次方程式Ａ×ｘ＝ｂが与えられている最も細かいグリッドを第１グリッドとし、第１グリッドよりも粗いグリッドを順番に第２グリッド、第３グリッドと階層化する。細かいグリッド上の変数を粗いグリッド上に変換する操作が粗視化処理であり、リストリクション（ｒｅｓｔｒｉｃｔｉｏｎ）とも呼ばれている。粗いグリッド上の変数を細かいグリッド上に変換する操作が精細化処理であり、プロロンゲーション（ｐｒｏｌｏｎｇａｔｉｏｎ）とも呼ばれている。

マルチグリッド法の典型的な処理例を図３を用いて説明する。
ステップＳ３１０では、与えられた第１グリッド（細かいグリッド）の格子分割数を２で割った値が格子分割数となる第２グリッド（粗いグリッド）を生成する。格子分割数が２で割れない場合は処理を終了する。
ステップＳ３２０では、解くべき連立一次方程式の係数行列（細かいグリッド上の係数行列）を粗視化して、粗いグリッド上の係数行列を生成する。ここでは、細かいグリッドの格子分割数と粗いグリッドの格子分割数の比が２になる場合の粗視化作用素と精細化作用素を用いる。
細かいグリッド上の係数行列Ａ_ｆから粗いグリッド上の係数行列Ａ_ｃを生成する方法としては、ＤＣＡとＧＣＡがある。ＤＣＡ（Ｄｉｒｅｃｔ Ｃｏａｒｓｅ ｇｒｉｄ Ａｐｐｒｏｘｉｍａｔｉｏｎ）は、物理系を表す基礎方程式を粗いグリッド上で直接離散化して連立一次方程式を導出する方法である。ＧＣＡ（Ｇａｌｅｒｋｉｎ Ｃｏａｒｓｅ ｇｒｉｄ Ａｐｐｒｏｘｉｍａｔｉｏｎ）は、細かいグリッド上の係数行列Ａ_ｆを粗視化作用素の行列形式Ｒと精細化作用素の行列形式Ｐで挟んで行列積Ｒ×Ａ_ｆ×Ｐを計算して導出する方法である。細かいグリッドの分割数をＮ_ｆ、粗いグリッドの分割数をＮ_ｃとすると、Ａ_ｆはＮ_ｆ行Ｎ_ｆ列の正方行列、ＲはＮ_ｃ行Ｎ_ｆ列の行列、ＰはＮ_ｆ行Ｎ_ｃ列の行列
であり、Ｒ×Ａ_ｆ×Ｐ＝Ａ_ｃはＮ_ｃ行Ｎ_ｃ列の正方行列になる。

ステップＳ３３０では、解くべき連立一次方程式Ａ×ｘ＝ｂの解ベクトルｘと定数項ベクトルｂについて、細かいグリッドと粗いグリッドとの格子分割数比が２になる場合の粗視化作用素を作用させて粗いグリッド上の連立一次方程式を準備する。
ステップＳ３４０では、粗いグリッド上の連立一次方程式を解くことによって、誤差の長波長成分を減衰させ易くする。ここで用いる求解法はガウス−ザイデル法やＳＯＲ法等、適当に選んで良い。
ステップＳ３５０では、粗いグリッド上で解いた連立一次方程式の解ベクトルについて、細かいグリッドの格子分割数と粗いグリッドの格子分割数の比が２になる場合の精細化作用素を作用させて、元の第１グリッド上の変数に戻す。
ステップＳ３６０では、初めに与えられた連立一次方程式を解くことによって、誤差の短波長成分を減衰させ易くする。このときステップＳ３５０で精細化された解ベクトルを初期値として用いることによって、誤差の長波長成分と短波長成分を同時に減衰させ易くすることができる。
図３の処理例は第１グリッドと第２グリッドのみ計算しているが、実用上はＳ３３０からＳ３６０の処理を再帰的に実行してより粗いグリッド上でも計算することによって、誤差の様々な波長成分を同時に減衰させ易くしている。

Ｓ３２０、Ｓ３３０、Ｓ３５０で使用する粗視化処理と精細化処理は、細かいグリッドの格子分割数と粗いグリッドの格子分割数の比が２になる場合が最も簡単であるため、格子分割数は２のべき乗になることが多い。また、特許文献２のように格子分割数が素因数分解可能な整数の積でもマルチグリッド法は適用可能である。この場合、Ｓ３１０では設定値以下の素数で格子分割数を割った数を粗いグリッドの格子分割数に設定し、Ｓ３２０、Ｓ３３０、Ｓ３５０で使用する粗視化処理と精細化処理は、当該素数に対応したものを使用する。

特開２００６−２９３４８８号公報 特開２００７−２８６７３８号公報

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.

しかし、上述した従来のマルチグリッド法の処理は、元の第１グリッドの格子分割数が２のべき乗または一般には素因数分解可能な整数の積という条件を満たす問題にしか適用できない。
解くべき問題の格子分割数がこの条件を満たさなければ、連立一次方程式の求解法として従来のマルチグリッド法が適用できないため、計算時間の長くかかる別の求解法を選択する必要がある。或いは従来のマルチグリッド法を適用するために前記条件を満たすよう格子分割数を変更する必要がある。そのため、余分な格子を追加することによって計算時間が増大したり、格子間隔が必要以上に狭まって物体表面と格子境界が揃わなかったりする場合がある。
また、解析空間の格子分割数が素数、例えば９７の場合、９７より大きい２のべき乗である１２８を用いて、元のグリッド上の物理量を格子分割数１２８のグリッド上の物理量に変換してからマルチグリッド法を用いて求解する。その後、格子分割数１２８のグリッド上から元の格子分割数９７のグリッド上に物理量を補間して戻すという処理を行って求解する方法も可能だが、補間に伴って誤差が生じる場合がある。

本発明はこのような背景のもとになされたものであり、任意の格子分割数の問題に対して適用可能なマルチグリッド法を用いた連立一次方程式の数値計算方法及び装置を提供することを目的とする。

本発明は、数値計算装置におけるマルチグリッド法を用いた連立一次方程式の数値計算方法であって、
前記数値計算装置が、
格子分割数Ｎ_ｆの細かいグリッド上のｊ番目の格子ａ_ｊにおける値ｘ_ｊから、格子分割数Ｎ_ｃ（＜Ｎ_ｆ）の粗いグリッド上のｉ番目の格子Ａ_ｉにおける値Ｘ_ｉへの変換を

により行う粗視化工程と、
格子分割数Ｎ_ｃの粗いグリッド上のｉ番目の格子Ｂ_ｉにおける値Ｙ_ｉから、格子分割数Ｎ_ｆ（＞Ｎ_ｃ）の細かいグリッド上のｊ番目の格子ｂ_ｊにおける値ｙ_ｊへの変換を

により行う精細化工程と、
を実行し、
前記粗視化工程は、
前記数値計算装置が、前記格子Ａ_ｉを被覆する細かいグリッド上の格子ａ_ｊの集合を格子グループ｛ａ｝_ｉとして計算するする工程と、
前記数値計算装置が、前記格子グループ｛ａ｝_ｉに属する格子ａ_ｊ毎に、該格子ａ_ｊと格子Ａ_ｉとの共通部分の長さ、面積又は体積と、該格子ａ_ｊの長さ、面積又は体積と、の比率を重なり分率ｐ_ｉｊとして計算する工程と、
前記数値計算装置が、前記格子グループ｛ａ｝_ｉに属する格子ａ_ｊ毎に、該格子ａ_ｊにおける値ｘ_ｊの前記格子Ａ_ｉにおける値Ｘ_ｉへの寄与を重み付けｗ_ｉｊとして所定の値に設定する工程と、
前記数値計算装置が、前記格子グループ｛ａ｝_ｉに属する格子ａ_ｊ毎に、該格子ａ_ｊにおける値ｘ_ｊ、重なり分率ｐ_ｉｊ及び重み付けｗ_ｉｊの積を計算する工程と、
前記数値計算装置が、前記積の前記格子グループ｛ａ｝_ｉに属する全ての格子に亘る総和を格子Ａ_ｉにおける値Ｘ_ｉとして計算する工程と、
を有し、
前記精細化工程は、
前記数値計算装置が、前記格子ｂ_ｊを被覆する粗いグリッド上の格子Ｂ_ｉの集合を格子グループ｛Ｂ｝_ｊとして計算する工程と、
前記数値計算装置が、前記格子グループ｛Ｂ｝_ｊに属する格子Ｂ_ｉ毎に、該格子Ｂ_ｉと格子ｂ_ｊとの共通部分の長さ、面積又は体積と、該格子Ｂ_ｉの長さ、面積又は体積と、の比率を重なり分率ｑ_ｊｉとして計算する工程と、
前記数値計算装置が、前記格子グループ｛Ｂ｝_ｊに属する格子Ｂ_ｉ毎に、該格子Ｂ_ｉにおける値Ｙ_ｉの前記格子ｂ_ｊにおける値ｙ_ｊへの寄与を重み付けｖ_ｊｉとして所定の値に設定する工程と、
前記数値計算装置が、前記格子グループ｛Ｂ｝_ｊに属する格子Ｂ_ｉ毎に、該格子Ｂ_ｉにおける値Ｙ_ｉ、重なり分率ｑ_ｊｉ及び重み付けｖ_ｊｉの積を計算する工程と、
前記数値計算装置が、前記積の前記格子グループ｛Ｂ｝_ｊに属する全ての格子に亘る総和を格子ｂ_ｊにおける値ｙ_ｊとして計算する工程と、
を有する
ことを特徴とするマルチグリッド法を用いた連立一次方程式の数値計算方法である。

本発明によれば、任意の格子分割数の問題に対して適用可能なマルチグリッド法を用いた連立一次方程式の数値計算方法及び装置を提供することができる。

実施例に係る数値計算手法のフローチャート マルチグリッド法における粗視化及び精細化を説明する図 従来のマルチグリッド法の計算処理例を示すフローチャート 格子分割数比が２の場合の粗視化及び精細化作用を示す概念図 格子分割数比が７／４の場合の粗視化及び精細化作用を示す概念図 実施例の数値計算装置の基本構成を示す図

図６は、本実施形態に係る数値計算装置として機能するコンピュータの基本構成を示すブロック図である。該コンピュータは、大略、バス２０７を介して接続されたＣＰＵ２０１、表示装置２０２，入力装置２０３、ＲＡＭ２０４、ハードディスク装置２０５、ＮＩＣ２０６等から構成される。
ＣＰＵ（中央演算装置）２０１で、ＲＡＭ（主記憶装置）２０４に記憶されているプログラムやデータを用いて本装置全体の制御を行うと共に、後述する各処理を実行する。表示装置２０２は、液晶画面などにより構成されており、ＣＰＵ２０１による処理結果を文字又は画像などで表示することができる。入力装置２０３は、キーボードやマウスなどの操作系の装置により構成されており、各種の指示をＣＰＵ２０１に対して入力することができる。ＲＡＭ２０４は、ハードディスク装置２０５からロードされたプログラムやデータを一時的に記憶する為のエリアを備えると共に、ＣＰＵ２０１が各種の処理を実行する際に必要なワークエリアを備える。ハードディスク装置２０５には、ＯＳ（オペレーティングシステム）や、コンピュータが行うべき後述の各処理をＣＰＵ２０１に実行させるためのプログラムやデータが保存されている。これらのプログラムはＲＡＭ２０４にロードされて、ＣＰＵ２０１によって実行されることで、本コンピュータが連立一次方程式を数値的に解くための数値計算装置として機能する。ＮＩＣ（ネットワークインターフェース）２０６は、外部装置とのデータ通信を行う際のインターフェースとして機能するものであり、本コンピュータはこのＮＩＣ２０６を介して外部装置とプログラムやデータの送受信を行う。なお、ハードディスク装置２０５が保存しているプログラムやデータの代わりに、ＮＩＣ２０６を介して外部装置から受信したプログラムやデータをＣＰＵ２０１による処理対象としても良い。

次に、上記構成を備えるコンピュータが行う、連立一次方程式の解を求める為の数値解法について説明する。図１は、本実施形態に係るコンピュータが行う、連立一次方程式の解を求めるための数値計算処理のフローチャートである。なお、ここでは係数行列Ａ_ｆ、解ベクトルｘ、定数項ベクトルｂとして、Ａ_ｆ×ｘ＝ｂで表される連立一次方程式の求解処理に、本発明のマルチグリッド法を適用する。
ここでは、１次元空間での物理現象シミュレーションを行うために基礎方程式を格子中心（セル中心）に値が定義された等間隔のグリッド上で３点差分により離散化し、係数行列が３重対角となっている連立一次方程式を解くことを想定している。
まず、ステップＳ２１０ではグリッドの格子分割数Ｎ_ｆについて偶数か奇数かを判定し、偶数の場合は２で割った値を、奇数の場合は所定の方法で算出した値を１段粗いグリッドの格子分割数Ｎ_ｃとして粗いグリッドを生成する。

なお、３次元空間の場合は、ｘ軸，ｙ軸，ｚ軸方向の格子分割数Ｎ_ｘ，Ｎ_ｙ，Ｎ_ｚについてそれぞれ偶奇を判定して粗いグリッドの各座標軸方向の格子分割数を計算して粗いグリッド生成すればよい。
Ｎ_ｆが奇数の場合のＮ_ｃの算出方法としては、Ｎ_ｆに１を加算してから２で割った値をＮ_ｃとするＮ_ｃ＝（Ｎ_ｆ＋１）÷２で表される方法や、逆にＮ_ｆから１を減算してから２で割った値をＮ_ｃとするＮ_ｃ＝（Ｎ_ｆ−１）÷２で表される方法がある。さらに半開区間［（Ｎ_ｆ−１）÷２，Ｎ_ｆ）に含まれる２のべき乗、２のべき乗と３の積、２のべき乗と５の積又は８の倍数のいずれかを選択することによってＮ_ｃを決定する方法がある。
例えばＮ_ｆが奇数の場合にＮ_ｃ＝（Ｎ_ｆ＋１）÷２とする方法では、Ｎ_ｆ＝７ならばＮ_ｃ＝（７＋１）÷２＝４、Ｎ_ｆ＝１７ならばＮ_ｃ＝（１７＋１）÷２＝９となる。
ステップＳ２２０では、ＤＣＡまたはＧＣＡによって解くべき連立一次方程式の係数行列Ａ_ｆを粗視化して、粗いグリッド上の係数行列Ａ_ｃを生成する。この際使用する粗視化作用素と精細化作用素は、ステップＳ２３０またはステップＳ２３１とステップＳ２５０またはステップＳ２５１のものと同等である。

ステップＳ２３０及びステップＳ２３１では、解くべき連立一次方程式Ａ_ｆ×ｘ＝ｂの解ベクトルｘと定数項ベクトルｂに粗視化作用素を作用させて粗いグリッド上の連立一次方程式を準備する。細かいグリッドの格子分割数Ｎ_ｆと粗いグリッドの格子分割数Ｎ_ｃの比Ｎ_ｆ÷Ｎ_ｃが整数の場合、ステップＳ２３０の粗視化作用素を作用させる。細かいグリッドの格子分割数Ｎ_ｆと粗いグリッドの格子分割数Ｎ_ｃの比Ｎ_ｆ÷Ｎ_ｃが非整数有理数の場合、ステップＳ２３１の粗視化作用素を作用させる。
本実施例では、細かいグリッドの格子分割数Ｎ_ｆが偶数の場合のみ格子分割数比Ｎ_ｆ÷Ｎ_ｃが整数（２）となり、ステップＳ２３０が実行される。細かいグリッドの格子分割数Ｎ_ｆが奇数の場合に格子分割数比Ｎ_ｆ÷Ｎ_ｃが非整数有理数となり、ステップＳ２３１が実行される。

図４（Ａ）は、ステップＳ２３０で使用する格子分割数比Ｎ_ｆ÷Ｎ_ｃが整数の場合の粗視化作用素、特に、細かいグリッドの格子分割数Ｎ_ｆが偶数の場合の粗視化作用素の粗視化作用を示す概念図である。図中の矢印に付随する数字は第１グリッド（細かいグリッド）上の値を第２グリッド（粗いグリッド）上の値に変換する際の分配比率を示している。
ステップＳ２３１で使用する格子分割数比Ｎ_ｆ÷Ｎ_ｃが非整数有理数の場合の粗視化作用素は、格子間の被覆状態を表す長さ分率ｐ_ｉｊを用いて式（１）で定義する。

ここで、Ｘ_ｉは粗いグリッド上のｉ番目の格子Ａ_ｉにおける値、ｘ_ｊは細かいグリッド上のｊ番目の格子ａ_ｊでの値、｛ａ｝_ｉは格子Ａ_ｉを被覆する細かいグリッド上の格子の集合である。

ｐ_ｉｊは格子Ａ_ｉを被覆する格子グループ｛ａ｝_ｉに属する各格子ａ_ｊと格子Ａ_ｉとの重なり具合を表す重なり分率である。重なり分率ｐ_ｉｊは、格子グループ｛ａ｝_ｉに属する格子ａ_ｊ毎に、該格子ａ_ｊと格子Ａ_ｉとの共通部分の長さ、面積又は体積と、該格子ａ
_ｊの長さ、面積又は体積と、の比率を計算して求められる。重なり分率ｐ_ｉｊは、細かいグリッドと粗いグリッドの幾何学的関係により各格子について一意に定まる。一次元の場合は、格子ａ_ｊの長さと、該格子ａ_ｊのうち格子Ａ_ｉに含まれる部分の長さと、の比率が、格子Ａ_ｉにおける格子ａ_ｊの重なり分率ｐ_ｉｊであり、これを長さ分率と称する。同様に、二次元の場合は、細かいグリッドの格子ａ_ｊの面積と、該格子ａ_ｊのうち粗いグリッドの格子Ａ_ｉに含まれる部分の面積と、の比率が重なり分率であり、これを面積分率と称する。三次元の場合は、細かいグリッドの格子ａ_ｊの体積と、該格子ａ_ｊのうち粗いグリッドの格子Ａ_ｉに含まれる部分の体積と、の比率が重なり分率であり、これを体積分率と称する。

ｗ_ｉｊは前記格子グループ｛ａ｝_ｉに属する各格子ａ_ｊにおける値ｘ_ｊの、格子Ａ_ｉの値Ｘ_ｉへの寄与を表す重み付けであり、任意の値に設定できる。重み付けｗ_ｉｊの設定方法としては、格子の位置に依存しない一定値Ｎ_ｃ／Ｎ_ｆとする方法や、格子Ａ_ｉと格子ａ_ｊの中心間距離の関数とする方法がある。また、格子Ａ_ｉと格子ａ_ｊの中心間距離が最短の格子以外は全てゼロに設定する方法、つまり格子Ａ_ｉに最も近い格子ａ_ｊ０における値のみを代表値として採用する方法がある。格子分割数比が非整数有理数の場合、上記の重なり分率ｐ_ｉｊ、重み付けｗ_ｉｊ及び格子グループに属する各格子ａ_ｊにおける値ｘ_ｊの積の、格子グループ｛ａ｝_ｉに属する全ての格子に亘る総和を、粗いグリッドの格子Ａ_ｉにおける値Ｘ_ｉとする。

図５（Ａ）は、ステップＳ２３１で使用する格子分割数比Ｎ_ｆ÷Ｎ_ｃが非整数有理数（７／４）の場合の粗視化作用素の粗視化作用を示す概念図である。図中の矢印に付随する数字は長さ分率に基づいて第１グリッド（細かいグリッド）上の値を第２グリッド（粗いグリッド）上の値に変換する際の分配比率を示している。式１に示すように、細かいグリッドの格子ａ_ｊにおける値ｘ_ｊが粗いグリッドの格子Ａ_ｉにおける値Ｘ_ｉに分配される比率は、格子Ａ_ｉにおける格子ａ_ｊの長さ分率ｐ_ｉｊと重み付けｗ_ｉｊとの積である。図５（Ａ）に示した分配比率は、重み付けｗ_ｉｊを一定値Ｎ_ｃ／Ｎ_ｆ＝４／７に設定した場合の値である。

例えば、粗視化作用により粗い第２グリッド上の３番目の格子Ａ_２の値Ｘ_２を求める方法は次のようになる。格子Ａ_２を被覆する細かい第１グリッド上の格子はａ_３，ａ_４，ａ_５の３個である。すなわち、｛ａ｝_２＝｛ａ_３，ａ_４，ａ_５｝である。格子ａ_３はその１／２の部分が格子Ａ_２に含まれるので、格子Ａ_２における格子ａ_３の長さ分率ｐ_２３は１／２である。格子ａ_４はその全部が格子Ａ_２に含まれるので、格子Ａ_２における格子ａ_４の長さ分率ｐ_２４は１である。格子ａ_５はその１／４の部分が格子Ａ_２に含まれるので、格子Ａ_２における格子ａ_５の長さ分率ｐ_２５は１／４である。ここでは、重み付けｗ_ｉｊは格子分割数比４／７で一定値とする。従って第２グリッド上の格子Ａ_２における値Ｘ_２は、Ｘ_２＝ｘ_３ｐ_２３ｗ_２３＋ｘ_４ｐ_２４ｗ_２４＋ｘ_５ｐ_２５ｗ_２５＝ｘ_３×（１／２）×（４／７）＋ｘ_４×１×（４／７）＋ｘ_５×（１／４）×（４／７）＝ｘ_３×２／７＋ｘ_４×４／７＋ｘ_５×１／７となる。

ここで、一般に粗いグリッド上のｉ番目の格子Ａ_ｉを被覆する細かいグリッド上の格子のグループ｛ａ｝_ｉに属する格子の格子番号を求めるには、格子位置を表す実数型の座標値から整数型に変換する処理または有理数型の計算処理が必要である。
しかし、Ｎ_ｃ＝（Ｎ_ｆ＋１）÷２の場合は、前記格子グループ｛ａ｝_ｉに属する格子の格子番号は、ｉ＝０の場合はｉ×２とｉ×２＋１、すなわち｛ａ｝_０＝｛ａ_０，ａ_１｝となる。また、ｉ＝Ｎ_ｃ−１の場合はｉ×２−１とｉ×２、すなわち｛ａ｝_Ｎｃ−１＝｛ａ_{２Ｎｃ−３}，ａ_{２Ｎｃ−２}｝となる。また、他の場合はｉ×２−１、ｉ×２及びｉ×２＋１、すなわち｛ａ｝_ｉ＝｛ａ_２ｉ−１，ａ_２ｉ，ａ_２ｉ＋１｝となり、簡単な整数型の計算のみで算出できる。また、Ｎ_ｃ＝（Ｎ_ｆ−１）÷２の場合は、前記格子グループ｛ａ｝
_ｉに属する格子の格子番号はｉ×２、ｉ×２＋１及びｉ×２＋２、すなわち｛ａ｝_ｉ＝｛ａ_２ｉ，ａ_２ｉ＋１，ａ_２ｉ＋２｝となり、簡単な整数型の計算のみで算出できる。したがって、Ｎ_ｃ＝（Ｎ_ｆ±１）÷２の場合には高速に格子グループ｛ａ｝_ｉに属する格子の格子番号を求めることが可能である。

また、ステップＳ２３１で使用する粗視化作用素の行列形式Ｒの１行当たりの非零要素数を数えると、Ｎ_ｃ≧（Ｎ_ｆ−１）÷２では３個以下だが、Ｎ_ｃ＜（Ｎ_ｆ−１）÷２では４個以上になる場合が生じてしまう。つまりＮ_ｃ＝（Ｎ_ｆ±１）÷２の場合は、粗視化作用素を作用させる時の演算量の増大を抑制できる。
ステップＳ２４０では、粗いグリッド上の連立一次方程式を解くことによって、誤差の長波長成分を減衰させ易くする。ここで用いる求解法はガウス−ザイデル法やＳＯＲ法等、適当に選んで良い。

ステップＳ２５０およびステップＳ２５１では、粗いグリッド上で解いた連立一次方程式の解ベクトルに精細化作用素を作用させて元の第１グリッド上の変数に戻す。細かいグリッドの格子分割数Ｎ_ｆと粗いグリッドの格子分割数Ｎ_ｃの比Ｎ_ｆ÷Ｎ_ｃが整数の場合、ステップＳ２５０の精細化作用素を作用させる。細かいグリッドの格子分割数Ｎ_ｆと粗いグリッドの格子分割数Ｎ_ｃの比Ｎ_ｆ÷Ｎ_ｃが非整数有理数の場合、ステップＳ２５１の精細化作用素を作用させる。本実施例では、細かいグリッドの格子分割数Ｎ_ｆが偶数の場合のみ格子分割数比Ｎ_ｆ÷Ｎ_ｃが整数（２）となり、ステップＳ２５０が実行される。細かいグリッドの格子分割数Ｎ_ｆが奇数の場合に格子分割数比Ｎ_ｆ÷Ｎ_ｃが非整数有理数となり、ステップＳ２５１が実行される。

図４（Ｂ）は、ステップＳ２５０で使用する格子分割数比Ｎ_ｆ÷Ｎ_ｃが整数の場合の精細化作用素、特に、細かいグリッドの格子分割数Ｎ_ｆが偶数の場合の精細化作用素の精細化作用を示す概念図である。図中の矢印に付随する数字は第２グリッド（粗いグリッド）上の値を第１グリッド（細かいグリッド）上の値に変換する際の分配比率を示している。
Ｓ２５１で使用する格子分割数比Ｎ_ｆ÷Ｎ_ｃが非整数有理数の場合の精細化作用素は、格子間の被覆状態を表す長さ分率ｑ_ｊｉを用いて式（２）で定義する。

ここで、ｙ_ｊは細かいグリッド上のｊ番目の格子ｂ_ｊにおける値、Ｙ_ｉは粗いグリッド上のｉ番目の格子Ｂ_ｉでの値、｛Ｂ｝_ｊは格子ｂ_ｊを被覆する粗いグリッド上の格子の集合である。

ｑ_ｊｉは格子ｂ_ｊを被覆する格子グループ｛Ｂ｝_ｊに属する各格子Ｂ_ｉと格子ｂ_ｊとの重なり具合を表す重なり分率である。重なり分率ｑ_ｊｉは、格子グループ｛Ｂ｝_ｊに属する格子Ｂ_ｉ毎に、該格子Ｂ_ｉと格子ｂ_ｊとの共通部分の長さ、面積又は体積と、該格子Ｂ_ｉの長さ、面積又は体積と、の比率を計算して求められる。重なり分率ｑ_ｊｉは、粗いグリッドと細かいグリッドとの幾何学的関係により各格子について一意に定まる。一次元の場合は、格子Ｂ_ｉの長さと、該格子Ｂ_ｉのうち格子ｂ_ｊに含まれる部分の長さと、の比率が、格子ｂ_ｊにおける格子Ｂ_ｉの重なり分率ｑ_ｊｉであり、これを長さ分率と称する。同様に、二次元の場合は、粗いグリッドの格子Ｂ_ｉの面積と、該格子Ｂ_ｉのうち細かいグリッドの格子ｂ_ｊに含まれる部分の面積と、の比率が重なり分率であり、これを面積分率と称する。三次元の場合は、粗いグリッドの格子Ｂ_ｉの体積と、該格子Ｂ_ｉのうち細かいグリッドの格子ｂ_ｊに含まれる部分の体積と、の比率が重なり分率であり、これを体積分率
と称する。

ｖ_ｊｉは前記格子グループ｛Ｂ｝_ｊに属する各格子Ｂ_ｉにおける値Ｙ_ｉの、格子ｂ_ｊの値ｙ_ｊへの寄与を表す重み付けであり、任意の値に設定できる。重み付けｖ_ｊｉの設定方法としては、格子の位置に依存しない一定値Ｎ_ｆ／Ｎ_ｃとする方法や、格子ｂ_ｊと格子Ｂ_ｉの中心間距離の関数とする方法がある。また、格子ｂ_ｊと格子Ｂ_ｉの中心間距離が最短の格子以外は全てゼロに設定する方法、つまり格子ｂ_ｊに最も近い格子Ｂ_ｉ０における値のみを代表値として採用する方法がある。格子分割数比が非整数有理数の場合、上記の重なり分率ｑ_ｊｉ、重み付けｖ_ｊｉ及び格子グループに属する各格子Ｂ_ｉにおける値Ｙ_ｉの積の、格子グループ｛Ｂ｝_ｊに属する全ての格子に亘る総和を、粗いグリッドの格子ｂ_ｊにおける値ｙ_ｊとする。

図５（Ｂ）は、ステップＳ２５１で使用する格子分割数比Ｎ_ｆ÷Ｎ_ｃが非整数有理数（７／４）の場合の精細化作用素の精細化作用を示す概念図である。図中の矢印に付随する数字は長さ分率に基づいて第２グリッド（粗いグリッド）上の値を第１グリッド（細かいグリッド）上の値に変換する際の分配比率を示している。式２に示すように、粗いグリッドの格子Ｂ_ｉにおける値Ｙ_ｉが細かいグリッドの格子ｂ_ｊにおける値ｙ_ｊに分配される比率は、格子ｂ_ｊにおける格子Ｂ_ｉの長さ分率ｑ_ｊｉと重み付けｖ_ｊｉとの積である。図５（Ｂ）に示した分配比率は、重み付けｖ_ｊｉを一定値Ｎ_ｆ／Ｎ_ｃ＝７／４に設定した場合の値である。
例えば、精細化作用により細かい第１グリッド上の６番目の格子ｂ_５の値ｙ_５を求める方法は次のようになる。格子ｂ_５を被覆する粗い第２グリッド上の格子はＢ_２，Ｂ_３の２個である。すなわち｛Ｂ｝_５＝｛Ｂ_２、Ｂ_３｝である。格子Ｂ_２はその１／７の部分が格子ｂ_５に含まれるので、格子ｂ_５における格子Ｂ_２の長さ分率ｑ_５２＝１／７である。格子Ｂ_３はその３／７の部分が格子ｂ_５に含まれるので、格子ｂ_５における格子Ｂ_３の長さ分率ｑ_５３＝３／７である。ここでは、重み付けｖ_ｊｉは格子分割数比７／４で一定値とする。従って第１グリッド上の格子ｂ_５における値ｙ_５は、ｙ_５＝Ｙ_２ｑ_５２ｖ_５２＋Ｙ_３ｑ_５３ｖ_５３＝Ｙ_２×（１／７）×（７／４）＋Ｙ_３×（３／７）×（７／４）＝Ｙ_２×１／４＋Ｙ_３×３／４となる。

ここで、一般に細かいグリッド上のｊ番目の格子ｂ_ｊを被覆する粗いグリッド上の格子のグループ｛Ｂ｝_ｊに属する格子の格子番号を求めるには、格子位置を表す実数型の座標値から整数型に変換する処理または有理数型の計算処理が必要である。
しかし、Ｎ_ｃ＝（Ｎ_ｆ＋１）÷２の場合は、前記格子グループ｛Ｂ｝_ｊに属する格子の格子番号は、ｊが偶数の場合はｊ÷２、すなわち｛Ｂ｝_ｊ＝｛Ｂ_ｊ／２｝となる。また、ｊが奇数の場合は（ｊ−１）÷２及び（ｊ−１）÷２＋１、すなわち｛Ｂ｝_ｊ＝｛Ｂ_{（ｊ−１）／２}、Ｂ_{（ｊ−１）／２＋１}｝となり、簡単な整数型の計算のみで算出できる。また、Ｎ_ｃ＝（Ｎ_ｆ−１）÷２の場合は、前記格子グループ｛Ｂ｝_ｊに属する格子の格子番号は、ｊ＝０の場合はｊ÷２、すなわち｛Ｂ｝_０＝｛Ｂ_０｝となる。また、ｊ＝Ｎ_ｆ−１の場合はｊ÷２−１、すなわち｛Ｂ｝_Ｎｆ−１＝｛Ｂ_{（Ｎｆ−１）／２−１}｝となる。また、ｊが０又はＮ_ｆ−１以外の偶数の場合はｊ÷２とｊ÷２−１、すなわち｛Ｂ｝_ｊ＝｛Ｂ_ｊ／２、Ｂ_{ｊ／２−１}｝となる。また、ｊが奇数の場合は（ｊ−１）÷２、すなわち｛Ｂ｝_ｊ＝｛Ｂ_{（ｊ−１）／２}｝となり、簡単な整数型の計算のみで算出できる。したがって、Ｎ_ｃ＝（Ｎ_ｆ±１）÷２の場合には高速に格子グループ｛Ｂ｝_ｊに属する格子の格子番号を求めることが可能である。

また、ステップＳ２５１で使用する精細化作用素の行列形式Ｐの１列当たりの非零要素数を数えると、Ｎ_ｃ≧（Ｎ_ｆ−１）÷２では最大３個以下だが、Ｎ_ｃ＜（Ｎ_ｆ−１）÷２では４個以上になる場合が生じてしまう。さらに粗いグリッド上の係数行列Ａ_ｃをＡ_ｃ＝Ｒ×Ａ_ｆ×Ｐとして計算するＧＣＡの場合、Ｎ_ｃ＞（Ｎ_ｆ＋１）÷２では細かいグリッド
上の係数行列Ａ_ｆが３重対角行列でもＡ_ｃの１行あたりの非零要素数が４個以上になる場合が生じてしまう。つまりＮ_ｃ＝（Ｎ_ｆ±１）÷２の場合は、ＧＣＡにおいて係数行列の３重対角性が保たれるので、計算の際の演算量の増大を抑制できる。
ステップＳ２６０では、初めに与えられた連立一次方程式を解くことによって、誤差の短波長成分を減衰させ易くする。このときステップＳ２５０またはステップＳ２５１で精細化された解ベクトルを初期値として用いることによって、誤差の長波長成分と短波長成分を同時に減衰させ易くすることができる。
図１のフローでは第１グリッドと第２グリッドのみ計算しているが、実用上はＳ２３０またはＳ２３１からＳ２６０の処理を再帰的に実行してより粗いグリッド上でも計算することによって、誤差の様々な波長成分を同時に減衰させ易くしている。

以上説明したように、本実施例に係るマルチグリッド法による連立一次方程式の数値解法によれば、格子分割数が２のべき乗ではない任意の格子に対して、マルチグリッド法により連立一次方程式を解くことが可能になる。
図１で示した数値計算方法は、共役勾配法、共役残差法を含む非定常反復解法の前処理として利用することも可能である。また、上述した連立一次方程式の数値計算方法は、格子分割数が２のべき乗でない２次元問題や３次元問題にも適用可能である。

本発明は、上述した連立一次方程式の数値計算方法を実現するソフトウェアのプログラムコードを記憶した記憶媒体であっても良い。該記憶媒体をシステム或いは装置に供給し、そのシステム或いは装置のコンピュータ（またはＣＰＵやＭＰＵ等）が記憶媒体に格納されたプログラムコードを読み出し実行することによっても本発明の目的は達成される。
この場合、記憶媒体から読み出されたプログラムコード自体が前述した実施の形態の機能を実現することになり、そのプログラムコード及び該プログラムコードを記憶した記憶媒体は本発明を構成することになる。
また、プログラムコードを供給するための記憶媒体としては、例えば、フロッピー（登録商標）ディスク、ハードディスク、光磁気ディスク、ＣＤやＤＶＤ等の光ディスク、磁気テープ、不揮発性のメモリカード、ＲＯＭ等を用いることができる。または、プログラムコードをネットワークを介してダウンロードしてもよい。

また、コンピュータが読み出したプログラムコードの指示に基づき、コンピュータ上で稼動しているＯＳ（オペレーティングシステム）等が実際の処理の一部または全部を行い、その処理によって上述の実施例の機能が実現される実施形態でも良い。
さらに、記憶媒体から読み出されたプログラムコードが、コンピュータに挿入された機能拡張ボードやコンピュータに接続された機能拡張ユニットに備わるメモリに書き込まれるようにしても良い。その後、そのプログラムコードの指示に基づき、その拡張機能を拡張ボードや拡張ユニットに備わるＣＰＵ等が実際の処理の一部または全部を行い、その処理によって上述の実施例の機能が実現される実施形態でも良い。

２０１…ＣＰＵ，２０２…表示装置，２０３…入力装置，２０４…ＲＡＭ

Claims (7)

1. 数値計算装置におけるマルチグリッド法を用いた連立一次方程式の数値計算方法であって、
   前記数値計算装置が、
   格子分割数Ｎ_ｆの細かいグリッド上のｊ番目の格子ａ_ｊにおける値ｘ_ｊから、格子分割数Ｎ_ｃ（＜Ｎ_ｆ）の粗いグリッド上のｉ番目の格子Ａ_ｉにおける値Ｘ_ｉへの変換を
   

   

   により行う粗視化工程と、
   格子分割数Ｎ_ｃの粗いグリッド上のｉ番目の格子Ｂ_ｉにおける値Ｙ_ｉから、格子分割数Ｎ_ｆ（＞Ｎ_ｃ）の細かいグリッド上のｊ番目の格子ｂ_ｊにおける値ｙ_ｊへの変換を
   

   

   により行う精細化工程と、
   を実行し、
   前記粗視化工程は、
   前記数値計算装置が、前記格子Ａ_ｉを被覆する細かいグリッド上の格子ａ_ｊの集合を格子グループ｛ａ｝_ｉとして計算するする工程と、
   前記数値計算装置が、前記格子グループ｛ａ｝_ｉに属する格子ａ_ｊ毎に、該格子ａ_ｊと格子Ａ_ｉとの共通部分の長さ、面積又は体積と、該格子ａ_ｊの長さ、面積又は体積と、の比率を重なり分率ｐ_ｉｊとして計算する工程と、
   前記数値計算装置が、前記格子グループ｛ａ｝_ｉに属する格子ａ_ｊ毎に、該格子ａ_ｊにおける値ｘ_ｊの前記格子Ａ_ｉにおける値Ｘ_ｉへの寄与を重み付けｗ_ｉｊとして所定の値に設定する工程と、
   前記数値計算装置が、前記格子グループ｛ａ｝_ｉに属する格子ａ_ｊ毎に、該格子ａ_ｊにおける値ｘ_ｊ、重なり分率ｐ_ｉｊ及び重み付けｗ_ｉｊの積を計算する工程と、
   前記数値計算装置が、前記積の前記格子グループ｛ａ｝_ｉに属する全ての格子に亘る総和を格子Ａ_ｉにおける値Ｘ_ｉとして計算する工程と、
   を有し、
   前記精細化工程は、
   前記数値計算装置が、前記格子ｂ_ｊを被覆する粗いグリッド上の格子Ｂ_ｉの集合を格子グループ｛Ｂ｝_ｊとして計算する工程と、
   前記数値計算装置が、前記格子グループ｛Ｂ｝_ｊに属する格子Ｂ_ｉ毎に、該格子Ｂ_ｉと格子ｂ_ｊとの共通部分の長さ、面積又は体積と、該格子Ｂ_ｉの長さ、面積又は体積と、の比率を重なり分率ｑ_ｊｉとして計算する工程と、
   前記数値計算装置が、前記格子グループ｛Ｂ｝_ｊに属する格子Ｂ_ｉ毎に、該格子Ｂ_ｉにおける値Ｙ_ｉの前記格子ｂ_ｊにおける値ｙ_ｊへの寄与を重み付けｖ_ｊｉとして所定の値に設定する工程と、
   前記数値計算装置が、前記格子グループ｛Ｂ｝_ｊに属する格子Ｂ_ｉ毎に、該格子Ｂ_ｉにおける値Ｙ_ｉ、重なり分率ｑ_ｊｉ及び重み付けｖ_ｊｉの積を計算する工程と、
   前記数値計算装置が、前記積の前記格子グループ｛Ｂ｝_ｊに属する全ての格子に亘る総和を格子ｂ_ｊにおける値ｙ_ｊとして計算する工程と、
   を有する
   ことを特徴とするマルチグリッド法を用いた連立一次方程式の数値計算方法。
2. 前記数値計算装置が、前記粗視化工程において重み付けｗ_ｉｊを格子の位置に依存しない一定値に設定し、
   前記数値計算装置が、前記精細化工程において重み付けｖ_ｊｉを格子の位置に依存しない一定値に設定する
   ことを特徴とする請求項１に記載のマルチグリッド法を用いた連立一次方程式の数値計算方法。
3. 前記数値計算装置が、前記粗視化工程において重み付けｗ_ｉｊを前記格子グループ｛ａ｝_ｉに属する格子ａ_ｊの中で前記格子Ａ_ｉとの中心間距離が最短の格子以外は全てゼロに設定し、
   前記数値計算装置が、前記精細化工程において重み付けｖ_ｊｉを前記格子グループ｛Ｂ｝_ｊに属する格子Ｂ_ｉの中で前記格子ｂ_ｊとの中心間距離が最短の格子以外は全てゼロに設定する
   ことを特徴とする請求項１に記載のマルチグリッド法を用いた連立一次方程式の数値計算方法。
4. 前記数値計算装置は、細かいグリッドから粗いグリッドを生成するグリッド生成工程を実行するものであって、
   前記グリッド生成工程は、
   前記細かいグリッドの各座標軸方向の格子分割数Ｎ_ｆが偶数の場合に、前記数値計算装置が該格子分割数Ｎｆを２で割った値を前記粗いグリッドの各座標軸方向の格子分割数Ｎ_ｃとして決定し格子分割数比Ｎ_ｆ÷Ｎ_ｃが整数となる粗いグリッドを生成する工程と、
   前記細かいグリッドの各座標軸方向の格子分割数Ｎ_ｆが奇数の場合に、前記数値計算装置が該格子分割数Ｎ_ｆに１を加算または減算してから２で割った値を粗いグリッドの各座標軸方向の格子分割数Ｎ_ｃとして決定し格子分割数比Ｎ_ｆ÷Ｎ_ｃが整数ではない有理数となる粗いグリッドを生成する工程と、
   を有する
   ことを特徴とする請求項１から３の何れか１項に記載のマルチグリッド法を用いた連立一次方程式の数値計算方法。
5. 前記数値計算装置が、前記グリッド生成工程において粗いグリッドの格子分割数Ｎ_ｃを（Ｎ_ｆ＋１）÷２に決定した場合、前記数値計算装置は、
   前記粗視化工程において前記格子グループ｛ａ｝_ｉを、
   （Ａ）ｉ＝０の場合、細かいグリッド上の格子番号がｉ×２及びｉ×２＋１の格子の集合、
   （Ｂ）ｉ＝Ｎ_ｃ−１の場合、細かいグリッド上の格子番号がｉ×２−１及びｉ×２の格子の集合、
   （Ｃ）それ以外の場合、細かいグリッド上の格子番号がｉ×２−１、ｉ×２及びｉ×２＋１の格子の集合、
   として計算し、
   前記精細化工程において前記格子グループ｛Ｂ｝_ｊを、
   （Ａ）ｊが偶数の場合、粗いグリッド上の格子番号がｊ÷２の格子の集合、
   （Ｂ）ｊが奇数の場合、粗いグリッド上の格子番号が（ｊ−１）÷２及び（ｊ−１）÷２＋１の格子の集合、
   として計算し、
   前記数値計算装置が、前記グリッド生成工程において粗いグリッドの格子分割数Ｎ_ｃを（Ｎ_ｆ−１）÷２に決定した場合、前記数値計算装置は、
   前記粗視化工程において前記格子グループ｛ａ｝_ｉを、細かいグリッド上の格子番号がｉ×２、ｉ×２＋１及びｉ×２＋２の格子の集合として計算し、
   前記精細化工程において前記格子グループ｛Ｂ｝_ｊを、
   （Ａ）ｊ＝０の場合、粗いグリッド上の格子番号がｊ÷２の格子の集合、
   （Ｂ）ｊ＝Ｎ_ｆ−１の場合、粗いグリッド上の格子番号がｊ÷２−１の格子の集合、
   （Ｃ）ｊが０またはＮ_ｆ−１以外の偶数の場合、粗いグリッド上の格子番号がｊ÷２及びｊ÷２−１の格子の集合、
   （Ｄ）ｊが奇数の場合、粗いグリッド上の格子番号が（ｊ−１）÷２の格子の集合、として計算する
   ことを特徴とする請求項４に記載のマルチグリッド法を用いた連立一次方程式の数値計算方法。
6. 前記数値計算装置は、細かいグリッドから粗いグリッドを生成するグリッド生成工程を実行するものであって、
   前記グリッド生成工程は、
   前記細かいグリッドの各座標軸方向の格子分割数Ｎ_ｆが偶数の場合に、前記数値計算装置が該格子分割数Ｎ_ｆを２で割った値を前記粗いグリッドの各座標軸方向の格子分割数Ｎ_ｃとして決定し格子分割数比Ｎ_ｆ÷Ｎ_ｃが整数となる粗いグリッドを生成する工程と、
   前記細かいグリッドの各座標軸方向の格子分割数Ｎ_ｆが奇数の場合に、半開区間［（Ｎ_ｆ−１）÷２，Ｎ_ｆ）に含まれる２のべき乗、２のべき乗と３の積、２のべき乗と５の積または８の倍数のいずれかを粗いグリッドの各座標軸方向の格子分割数Ｎ_ｃとして決定し格子分割数比Ｎ_ｆ÷Ｎ_ｃが整数ではない有理数となる粗いグリッドを生成する工程と、を有する
   ことを特徴とする請求項１から３の何れか１項に記載のマルチグリッド法を用いた連立一次方程式の数値計算方法。
7. マルチグリッド法を用いて連立一次方程式を数値計算する数値計算装置であって、
   格子分割数Ｎ_ｆの細かいグリッド上のｊ番目の格子ａ_ｊにおける値ｘ_ｊから、格子分割数Ｎ_ｃ（＜Ｎ_ｆ）の粗いグリッド上のｉ番目の格子Ａ_ｉにおける値Ｘ_ｉへの変換を
   

   

   により行う粗視化手段と、
   格子分割数Ｎ_ｃの粗いグリッド上のｉ番目の格子Ｂ_ｉにおける値Ｙ_ｉから、格子分割数Ｎ_ｆ（＞Ｎ_ｃ）の細かいグリッド上のｊ番目の格子ｂ_ｊにおける値ｙ_ｊへの変換を
   

   

   により行う精細化手段と、
   を備え、
   前記粗視化手段は、
   前記格子Ａ_ｉを被覆する細かいグリッド上の格子ａ_ｊの集合を格子グループ｛ａ｝_ｉとして計算するする手段と、
   前記格子グループ｛ａ｝_ｉに属する格子ａ_ｊ毎に、該格子ａ_ｊと格子Ａ_ｉとの共通部分の長さ、面積又は体積と、該格子ａ_ｊの長さ、面積又は体積と、の比率を重なり分率ｐ_ｉｊとして計算する手段と、
   前記格子グループ｛ａ｝_ｉに属する格子ａ_ｊ毎に、該格子ａ_ｊにおける値ｘ_ｊの前記格子Ａ_ｉにおける値Ｘ_ｉへの寄与を重み付けｗ_ｉｊとして所定の値に設定する手段と、
   前記格子グループ｛ａ｝_ｉに属する格子ａ_ｊ毎に、該格子ａ_ｊにおける値ｘ_ｊ、重な
   り分率ｐ_ｉｊ及び重み付けｗ_ｉｊの積を計算する手段と、
   前記積の前記格子グループ｛ａ｝_ｉに属する全ての格子に亘る総和を格子Ａ_ｉにおける値Ｘ_ｉとして計算する手段と、
   を備え、
   前記精細化手段は、
   前記格子ｂ_ｊを被覆する粗いグリッド上の格子Ｂ_ｉの集合を格子グループ｛Ｂ｝_ｊとして計算する手段と、
   前記格子グループ｛Ｂ｝_ｊに属する格子Ｂ_ｉ毎に、該格子Ｂ_ｉと格子ｂ_ｊとの共通部分の長さ、面積又は体積と、該格子Ｂ_ｉの長さ、面積又は体積と、の比率を重なり分率ｑ_ｊｉとして計算する手段と、
   前記格子グループ｛Ｂ｝_ｊに属する格子Ｂ_ｉ毎に、該格子Ｂ_ｉにおける値Ｙ_ｉの前記格子ｂ_ｊにおける値ｙ_ｊへの寄与を重み付けｖ_ｊｉとして所定の値に設定する手段と、
   前記格子グループ｛Ｂ｝_ｊに属する格子Ｂ_ｉ毎に、該格子Ｂ_ｉにおける値Ｙ_ｉ、重なり分率ｑ_ｊｉ及び重み付けｖ_ｊｉの積を計算する手段と、
   前記積の前記格子グループ｛Ｂ｝_ｊに属する全ての格子に亘る総和を格子ｂ_ｊにおける値ｙ_ｊとして計算する手段と、
   を備えることを特徴とする数値計算装置。

Images