(57)【要約】【課題】  ＡＭＧ法を内部ソルバに採用した残差切除法を用いた数値計算において、収束性を高める。【解決手段】  変数分離型ＡＭＧ法を利用する。ＡＭＧサイクルとして、ν−サイクルを用いる。ν−サイクルでは、緩和の反復回数ｎが、グリッドレベルを降りる場合は、レベルが浅いときは相対的に小さく、レベルが深いときは相対的に大きく設定される。一方、グリッドレベルを上がる場合は、緩和の反復回数ｎは所定の値ｎuに設定されている。これにより、少ない計算量で、Ｖ−サイクルに近い結果を得ることができる。【選択図】  図４
2005267418.tif 000002
【特許請求の範囲】
【請求項１】
物理量Ｕが満たすべき偏微分方程式を離散化したときの係数行列（Ｎ行Ｎ列、Ｎは正の整数）をＡ、非斉次項（ソース項）をｆとしたとき、コンピュータを用いて、
Ａ・Ｕ＝ｆ
を解き、物理量Ｕの数値計算を行う方法であって、
物理量Ｕの初期値Ｕ0を設定し、
繰り返し回数ｍに初期値として０を、摂動量φの初期値として０を与え、残差ｒの初期値ｒ0として（ｆ−Ａ・Ｕ0）を設定し、
繰り返し回数ｍをインクリメントしながら、
Ａ・φ＝ｒmの予測近似値Ψmを、内部ソルバを有する第１演算部によって、反復計算によって求める第１のステップと、
残差ｒm+1のＬ２ノルムを最小とする修正近似値φmを、第２演算部によって、最適化ルーチンによって、予測近似値Ψmから求め、近似解Ｕm+1として（Ｕm＋φm）を与え、残差ｒm+1として（ｒm−Ａ・φm）を与える第２のステップとを、近似解Ｕmが収束するまで、繰り返し実行するものであり、
前記第１のステップは、
前記内部ソルバにＡＭＧ法を用いるものであり、かつ、
ＡＭＧサイクルとして、緩和の反復回数ｎが、グリッドレベルを降りる場合は、レベルが相対的に浅いときは相対的に小さく、レベルが相対的に深いときは相対的に大きく設定される一方、グリッドレベルを上がる場合は、所定の値に設定されたν−サイクルを用いる
ことを特徴とする数値計算方法。
【請求項２】
請求項１において、
前記ν−サイクルは、緩和法の反復回数ｎが、
１）レベルｋからレベルｋ＋１に「制限」する前に行う緩和については、
ｎ＝ｍａｘ（ｋ−ｎb ，０）
２）レベルｋ＋１からレベルｋに「補間」する前に行う緩和については、
ｎ＝ｎu
と設定されているものである（ｎb ，ｎu は自然数）
ことを特徴とする数値計算方法。
【請求項３】
物理量Ｕが満たすべき偏微分方程式を離散化したときの係数行列（Ｎ行Ｎ列、Ｎは正の整数）をＡ、非斉次項（ソース項）をｆとしたとき、コンピュータを用いて、
Ａ・Ｕ＝ｆ
を解き、物理量Ｕの数値計算を行う方法であって、
物理量Ｕの初期値Ｕ0を設定し、
繰り返し回数ｍに初期値として０を、摂動量φの初期値として０を与え、残差ｒの初期値ｒ0として（ｆ−Ａ・Ｕ0）を設定し、
繰り返し回数ｍをインクリメントしながら、
Ａ・φ＝ｒmの予測近似値Ψmを、内部ソルバを有する第１演算部によって、反復計算によって求める第１のステップと、
残差ｒm+1のＬ２ノルムを最小とする修正近似値φmを、第２演算部によって、最適化ルーチンによって、予測近似値Ψmから求め、近似解Ｕm+1として（Ｕm＋φm）を与え、残差ｒm+1として（ｒm−Ａ・φm）を与える第２のステップとを、近似解Ｕmが収束するまで、繰り返し実行するものであり、
残差切除率が所定値よりも小さくなったとき、前記修正近似解φm を破棄し、前記第１のステップにおいて用いる内部ソルバのアルゴリズムの変更を伴うリセット処理を行う
ことを特徴とする数値計算方法。
【請求項４】
請求項３において、
前記第１のステップは、前記内部ソルバとして、ＡＭＧ法を用いるものであり、
前記リセット処理において、前記内部ソルバのアルゴリズムの変更として、ＡＭＧサイクルにおける緩和の反復回数ｎの設定の変更を、行うものとする
ことを特徴とする数値計算方法。
【請求項５】
物理量Ｕが満たすべき偏微分方程式を離散化したときの係数行列（Ｎ行Ｎ列、Ｎは正の整数）をＡ、非斉次項（ソース項）をｆとしたとき、コンピュータを用いて、
Ａ・Ｕ＝ｆ
を解き、物理量Ｕの数値計算を行う装置であって、
物理量Ｕの初期値Ｕ0を設定し、
繰り返し回数ｍに初期値として０を、摂動量φの初期値として０を与え、残差ｒの初期値ｒ0として（ｆ−Ａ・Ｕ0）を設定し、
繰り返し回数ｍをインクリメントしながら、
Ａ・φ＝ｒmの予測近似値Ψmを、内部ソルバを有する第１演算部によって、反復計算によって求める第１のステップと、
残差ｒm+1のＬ２ノルムを最小とする修正近似値φmを、第２演算部によって、最適化ルーチンによって、予測近似値Ψmから求め、近似解Ｕm+1として（Ｕm＋φm）を与え、残差ｒm+1として（ｒm−Ａ・φm）を与える第２のステップとを、近似解Ｕmが収束するまで、繰り返し実行するものであり、
前記第１のステップは、
前記内部ソルバにＡＭＧ法を用いるものであり、かつ、
ＡＭＧサイクルとして、緩和の反復回数ｎが、グリッドレベルを降りる場合は、レベルが相対的に浅いときは相対的に小さく、レベルが相対的に深いときは相対的に大きく設定される一方、グリッドレベルを上がる場合は、所定の値に設定されたν−サイクルを用いる
ことを特徴とする数値計算装置。
【請求項６】
物理量Ｕが満たすべき偏微分方程式を離散化したときの係数行列（Ｎ行Ｎ列、Ｎは正の整数）をＡ、非斉次項（ソース項）をｆとしたとき、コンピュータを用いて、
Ａ・Ｕ＝ｆ
を解き、物理量Ｕの数値計算を行う装置であって、
物理量Ｕの初期値Ｕ0を設定し、
繰り返し回数ｍに初期値として０を、摂動量φの初期値として０を与え、残差ｒの初期値ｒ0として（ｆ−Ａ・Ｕ0）を設定し、
繰り返し回数ｍをインクリメントしながら、
Ａ・φ＝ｒmの予測近似値Ψmを、内部ソルバを有する第１演算部によって、反復計算によって求める第１のステップと、
残差ｒm+1のＬ２ノルムを最小とする修正近似値φmを、第２演算部によって、最適化ルーチンによって、予測近似値Ψmから求め、近似解Ｕm+1として（Ｕm＋φm）を与え、残差ｒm+1として（ｒm−Ａ・φm）を与える第２のステップとを、近似解Ｕmが収束するまで、繰り返し実行するものであり、
残差切除率が所定値よりも小さくなったとき、前記修正近似解φm を破棄し、前記第１のステップにおいて用いる内部ソルバのアルゴリズムの変更を伴うリセット処理を行う
ことを特徴とする数値計算装置。
【請求項７】
コンピュータに、物理量Ｕが満たすべき偏微分方程式を離散化したときの係数行列（Ｎ行Ｎ列、Ｎは正の整数）をＡ、非斉次項（ソース項）をｆとしたとき、
Ａ・Ｕ＝ｆ
を解かせ、物理量Ｕの数値計算を行う数値計算用プログラムであって、
物理量Ｕの初期値Ｕ0を設定し、
繰り返し回数ｍに初期値として０を、摂動量φの初期値として０を与え、残差ｒの初期値ｒ0として（ｆ−Ａ・Ｕ0）を設定し、
繰り返し回数ｍをインクリメントしながら、
Ａ・φ＝ｒmの予測近似値Ψmを、内部ソルバを有する第１演算部によって、反復計算によって求める第１のステップと、
残差ｒm+1のＬ２ノルムを最小とする修正近似値φmを、第２演算部によって、最適化ルーチンによって、予測近似値Ψmから求め、近似解Ｕm+1として（Ｕm＋φm）を与え、残差ｒm+1として（ｒm−Ａ・φm）を与える第２のステップとを、近似解Ｕmが収束するまで、繰り返し実行するものであり、
前記第１のステップは、
前記内部ソルバにＡＭＧ法を用いるものであり、かつ、
ＡＭＧサイクルとして、緩和の反復回数ｎが、グリッドレベルを降りる場合は、レベルが相対的に浅いときは相対的に小さく、レベルが相対的に深いときは相対的に大きく設定される一方、グリッドレベルを上がる場合は、所定の値に設定されたν−サイクルを用いる
処理をコンピュータに実行させる数値計算用プログラム。
【請求項８】
コンピュータに、物理量Ｕが満たすべき偏微分方程式を離散化したときの係数行列（Ｎ行Ｎ列、Ｎは正の整数）をＡ、非斉次項（ソース項）をｆとしたとき、
Ａ・Ｕ＝ｆ
を解かせ、物理量Ｕの数値計算を行う数値計算用プログラムであって、
物理量Ｕの初期値Ｕ0を設定し、
繰り返し回数ｍに初期値として０を、摂動量φの初期値として０を与え、残差ｒの初期値ｒ0として（ｆ−Ａ・Ｕ0）を設定し、
繰り返し回数ｍをインクリメントしながら、
Ａ・φ＝ｒmの予測近似値Ψmを、内部ソルバを有する第１演算部によって、反復計算によって求める第１のステップと、
残差ｒm+1のＬ２ノルムを最小とする修正近似値φmを、第２演算部によって、最適化ルーチンによって、予測近似値Ψmから求め、近似解Ｕm+1として（Ｕm＋φm）を与え、残差ｒm+1として（ｒm−Ａ・φm）を与える第２のステップとを、近似解Ｕmが収束するまで、繰り返し実行するものであり、
残差切除率が所定値よりも小さくなったとき、前記修正近似解φm を破棄し、前記第１のステップにおいて用いる内部ソルバのアルゴリズムの変更を伴うリセット処理を行う
処理をコンピュータに実行させる数値計算用プログラム。
【発明の詳細な説明】
【技術分野】
【０００１】
本発明は、コンピュータを用いた数値計算技術に関するものである。
【背景技術】
【０００２】
従来から、圧力、温度などの物理量を対象とした制御や数値計算において、逐次近似解法アルゴリズムの利用がなされている。逐次近似解法アルゴリズムとは、様々な最適制御理論に基づいて、逐次近似を繰り返し行いつつ、数値計算により最適解を求めるものである。
【０００３】
そして、対象となる各パラメータの解を、高速かつ高精度に求めることができる収束性の高い逐次近似解法アルゴリズムとして、残差切除法が知られている。そして、特許文献１では、残差切除法において収束性の更なる向上を実現する技術が提案されている。
【０００４】
また、緩和法（ガウスザイデル法、ＳＯＲ法など）の加速手法の１つとして、マルチグリッド（Multi Grid：ＭＧ）法がある。一般に、緩和法は、残差の高周波成分を速く除去することができるが、低周波成分を減衰させるためには、非常に多くの反復回数を必要とする。ＭＧ法では、格子間隔（粗さ）の異なる複数の格子を用意し、各格子上で緩和法を用いることによって、各格子における残差の高周波成分を除去し、解を高速で求める。ある格子における残差の低周波成分は、より粗い格子においては高周波成分となる。
【０００５】
そして、ＭＧ法の派生手法として、代数的マルチグリッド（Algebraic Multi Grid：ＡＭＧ）法がある。上述のＭＧ法（ＡＭＧ法と区別するために幾何的マルチグリッド法と呼ばれることもある）を用いるためには、格子系を予め人手で準備しておく必要がある。これに対して、ＡＭＧ法を用いる場合は、格子系を与える必要はなく、連立一次方程式を与えるだけでよい。ＡＭＧ法では、与えられた連立一次方程式の係数行列から仮想の格子を作成し、疎グリッドの作成もまた係数行列の情報のみから行う。
【特許文献１】特開２００１−１３４３０４号公報
【非特許文献１】Andrew J. Cleary, et al. "ROBUSTNESS AND SCALABILITY OF ALGEBRAIC MULTIGRID", SIAM J. SCI. COMPUT.,2000,Vol.21,No.5, pp.1886-1908
【発明の開示】
【発明が解決しようとする課題】
【０００６】
本願発明者らは、上述のＡＭＧ法を、残差切除法の内部ソルバに利用することを検討している。
【０００７】
一方、構造解析や流体解析から得られる連立一次方程式では、未知変数に複数の物理量を含む場合がある。このような場合には、格子間の関係に物理量間の関係が加わるため、解くべき連立一次方程式の係数行列には、単一物理量の場合に比べて非零要素が多く含まれる。このため、問題が解きにくい（スティッフである）ことが多い。
【０００８】
このようなスティッフな問題に対して、上述したような、内部ソルバにＡＭＧ法を採用した残差切除法を単に適用すると、収束解が得られなかったり、収束解を得るまでに多大な計算時間を要したりする場合がある。
【０００９】
前記の問題に鑑み、本発明は、コンピュータによる逐次近似解法アルゴリズムを用いた数値計算、特に、ＡＭＧ法を内部ソルバに採用した残差切除法を用いた数値計算において、収束性を高めることを課題とする。
【課題を解決するための手段】
【００１０】
前記の課題を解決するために、第１の発明が講じた解決手段は、物理量Ｕが満たすべき偏微分方程式を離散化したときの係数行列（Ｎ行Ｎ列、Ｎは正の整数）をＡ、非斉次項（ソース項）をｆとしたとき、コンピュータを用いて、Ａ・Ｕ＝ｆを解き、物理量Ｕの数値計算を行うための方法、装置、およびプログラムとして、物理量Ｕの初期値Ｕ0を設定し、繰り返し回数ｍに初期値として０を、摂動量φの初期値として０を与え、残差ｒの初期値ｒ0として（ｆ−Ａ・Ｕ0）を設定し、繰り返し回数ｍをインクリメントしながら、Ａ・φ＝ｒmの予測近似値Ψmを、内部ソルバを有する第１演算部によって反復計算によって求める第１のステップと、残差ｒm+1のＬ２ノルムを最小とする修正近似値φmを、第２演算部によって最適化ルーチンによって予測近似値Ψmから求め、近似解Ｕm+1として（Ｕm＋φm）を与え、残差ｒm+1として（ｒm−Ａ・φm）を与える第２のステップとを、近似解Ｕmが収束するまで繰り返し実行するものであり、前記第１のステップは、前記内部ソルバにＡＭＧ法を用いるものであり、かつ、ＡＭＧサイクルとして、緩和の反復回数ｎが、グリッドレベルを降りる場合は、レベルが相対的に浅いときは相対的に小さく、レベルが相対的に深いときは相対的に大きく設定される一方、グリッドレベルを上がる場合は、所定の値に設定されたν−サイクルを用いるものである。
【００１１】
そして、前記第１の発明におけるν−サイクルは、緩和法の反復回数ｎが、１）レベルｋからレベルｋ＋１に「制限」する前に行う緩和については、ｎ＝ｍａｘ（ｋ−ｎb ，０）、２）レベルｋ＋１からレベルｋに「補間」する前に行う緩和については、ｎ＝ｎu と設定されている（ｎb ，ｎu は自然数）のが好ましい。
【００１２】
また、第２の発明が講じた解決手段は、物理量Ｕが満たすべき偏微分方程式を離散化したときの係数行列（Ｎ行Ｎ列、Ｎは正の整数）をＡ、非斉次項（ソース項）をｆとしたとき、コンピュータを用いて、Ａ・Ｕ＝ｆを解き、物理量Ｕの数値計算を行うための方法、装置、およびプログラムであって、物理量Ｕの初期値Ｕ0を設定し、繰り返し回数ｍに初期値として０を、摂動量φの初期値として０を与え、残差ｒの初期値ｒ0として（ｆ−Ａ・Ｕ0）を設定し、繰り返し回数ｍをインクリメントしながら、Ａ・φ＝ｒmの予測近似値Ψmを、内部ソルバを有する第１演算部によって反復計算によって求める第１のステップと、残差ｒm+1のＬ２ノルムを最小とする修正近似値φmを、第２演算部によって最適化ルーチンによって予測近似値Ψmから求め、近似解Ｕm+1として（Ｕm＋φm）を与え、残差ｒm+1として（ｒm−Ａ・φm）を与える第２のステップとを、近似解Ｕmが収束するまで繰り返し実行するものであり、残差切除率が所定値よりも小さくなったとき、前記修正近似解φm を破棄し、前記第１のステップにおいて用いる内部ソルバのアルゴリズムの変更を伴うリセット処理を行うものである。
【００１３】
そして、前記第２の発明において、前記第１のステップは、前記内部ソルバとしてＡＭＧ法を用いるものであり、前記リセット処理において、前記内部ソルバのアルゴリズムの変更として、ＡＭＧサイクルにおける緩和の反復回数ｎの設定の変更を行うのが好ましい。
【発明の効果】
【００１４】
第１の発明によると、残差切除法の内部ソルバにおけるＡＭＧサイクルとして、ν−サイクルを用いる。ν−サイクルでは、グリッドレベルを降りる際には、「制限」前に行う緩和の影響が小さいと考えられる浅いレベルでは、緩和の反復回数が少なく設定され、深いレベルでは多く設定されている。一方、グリッドレベルを昇る際には、「補間」後の緩和の反復回数は、Ｖ−サイクルと同様に、所定値に設定されている。したがって、より少ない計算量によって、Ｖ−サイクルとほぼ同等の性能を得ることができる。
【００１５】
第２の発明によると、残差切除率が所定値よりも小さくなったとき、リセット処理が行われ、内部ソルバのアルゴリズムが変更される。これにより、解の修正方向を見出す能力がなくなった状態から抜けだし、それまでとは異なる解修正方向を見出す可能性を高めることができる。したがって、収束性をより向上させることができる。
【発明を実施するための最良の形態】
【００１６】
以下、本発明の実施形態を図面に基づいて詳細に説明する。
【００１７】
本実施形態に係る数値計算方法は、内部ソルバにＡＭＧ法を利用した残差切除法によるものである。
【００１８】
＜ＡＭＧ法＞
まず初めに、ＡＭＧ法の概要について、説明する。
【００１９】
ＡＭＧ法の処理は、ＡＭＧセットアップと、ＡＭＧサイクルとの２つに大別される。ＡＭＧセットアップでは、疎グリッドの生成、疎密グリッド間の「制限」「補間」関数の生成、および疎グリッドにおける係数行列の生成が、行われる。これらの詳細については、上述した非特許文献１に詳述されているので、ここではその説明を省略する。
【００２０】
（ＡＭＧサイクル）
ＡＭＧサイクルは、ＡＭＧセットアップにおいて生成した疎グリッド、「制限」「補間」関数、および係数行列を用いて、連立一次方程式を解くフェーズである。最も簡単なサイクルは、グリッドが２段のみの場合であり、次の手順Ａ〜Ｅで１サイクルとなる。
【００２１】
Ａ．Ａ0 Ｕ0 ＝ｆ0 に対して、ｎ回緩和法（ガウスザイデル法など）を適用する（制限前緩和）
Ｂ．制限関数Ｉ01を用いて、ｆ1 ＝Ｉ01（ｆ0 −Ａ0 Ｕ0 ）を計算する。
【００２２】
Ｃ．Ａ1 Ｕ1 ＝ｆ1 を解いて、Ｕ1 を求める。
【００２３】
Ｄ．補間関数Ｉ10を用いて、近似解の補正Ｕ0 ←Ｕ0 ＋Ｉ10Ｕ1 を行う。
【００２４】
Ｅ．Ａ0 Ｕ0 ＝ｆ0 に対してｎ回緩和法を適用する（補間後緩和）。
【００２５】
ここでの上付添字は、グリッドレベルを示している。添字が０の式、すなわちグリッドレベル０における式が、上述の添字なしの式  Ａ・Ｕ＝ｆ  に対応している。なお、紛らわしいが、ＡＭＧサイクル以外での上付添字は、反復回数を表している。
【００２６】
上の手順Ａ〜Ｅのサイクルは、図１のように表現される。図１において、縦方向はグリッドレベル、横方向は処理順序である。手順Ａ，Ｅは同じグリッドレベル０における処理であり、手順Ｃはグリッドレベル１における処理である。また、手順Ｂは「制限」の関数によってグリッドレベルをレベル０からレベル１に遷移する処理、手順Ｄは「補間」の関数によってグリッドレベルをレベル１からレベル０に遷移する処理である。図１に示すように、このような処理は、点と線とで表される形状が文字「Ｖ」に似ていることから、「Ｖ−サイクル」と呼ばれる。
【００２７】
また、グリッドが３段以上の場合のＶ−サイクルも、２段の場合と同様である。すなわち、手順Ａ，Ｂを繰り返すことによってグリッドレベルを下がって行き、最深グリッドで解を求め（手順Ｃ）、手順Ｄ，Ｅを繰り返すことによってグリッドレベルを昇って行き、１サイクルを終える。
【００２８】
また、Ｖ−サイクル以外にも、Ｗ−サイクルや鋸歯−サイクルなどが提案されている。図２はグリッドが３段の場合のＷ−サイクル（ａ）と鋸歯−サイクル（ｂ）を示す概念図である。
【００２９】
（変数分離型ＡＭＧ法）
上述したように、ＡＭＧ法では、疎グリッドの生成（格子の間引き）を、係数行列という情報に基づいて行う。解くべき連立一次方程式の元になる微分方程式に未知量が１種類しかない場合、すなわち、１つのグリッド上に１つの未知量しか存在しない場合には、係数行列からグリッド情報を類推できる場合が多い。しかし、多変量の問題を扱う場合には、係数行列のみからグリッド情報を読みとることは一般には困難である。この困難を回避するために考案された手法が、変数分離型ＡＭＧ法である。
【００３０】
変数分離型ＡＭＧ法では、疎グリッドの生成を行うために、係数行列以外に、未知量の種類に関する情報を用いる。具体的な手順は次のとおりである。
【００３１】
Ａ．密グリッドにおいて、未知量を種類毎に分離し、係数行列において種類の異なる未知量間の相互作用を表す項を無視することによって、各種類の未知量毎に連立一次方程式を構築する。
【００３２】
Ｂ．一般的なＡＭＧセットアップのアルゴリズムを用いて、各種類の未知量毎に、疎グリッド、および疎密グリッド間の「制限」「補間」関数を生成する。
【００３３】
Ｃ．手順Ｂにおいて生成した疎グリッド、および「制限」「補間」関数、並びに未知量間の相互作用を表す項を無視しない元の係数行列を用いて、疎グリッド上における連立一次方程式の係数行列を構築する。
【００３４】
図３は本実施形態に係る数値計算方法を示すフローチャートである。また図４は本実施形態に係る数値計算装置１の構成を示すブロック図である。本実施形態では、解修正ベクトルΨを求める残差切除法の内部ソルバに、上述の変数分離型ＡＭＧ法を利用する。この数値計算方法は、例えばＣＤ−ＲＯＭ２のような記録媒体に記録されている、当該方法を実現するためのプログラムを、コンピュータに実行させることによって、実施することができる。また、当該方法を実現するためのプログラムを実行するコンピュータを備えた装置によって、容易に実現することができる。
【００３５】
本実施形態に係る数値計算方法には、従来にはない特徴が、いくつか含まれている。その１つが、「ν−ＡＭＧサイクル」であり、もう１つが、「リセット発生時のＡＭＧサイクルの変更」である。
【００３６】
＜ν−ＡＭＧサイクル＞
本実施形態では、内部ソルバにおける１回の処理において、次のようなサイクル１サイクルのみを用いるものとする。
【００３７】
当初未知変数が属しているグリッドのレベルをレベル０とし、グリッドが粗くなるにつれてレベル１，レベル２，…とグリッド番号が増加するものとする。ｎb ，ｎu を自然数とし、緩和法の反復回数をｎとしたとき、
１）レベルｋからレベルｋ＋１に「制限」する前に行う緩和については、
ｎ＝ｍａｘ（ｋ−ｎb ，０）
２）レベルｋ＋１からレベルｋに「補間」する前に行う緩和については、
ｎ＝ｎu
とする。
【００３８】
図５は本実施形態において用いるＡＭＧサイクルの一例である。図５では、６段のグリッドを用い、ｎb ＝１としている。レベル０からレベル１への制限、およびレベル１からレベル２への制限においては、緩和法の反復回数ｎ＝０となっている。本願発明者は、このようなＡＭＧサイクルを、「ν−サイクル」と名付けた。
【００３９】
本願発明者が考案したν−サイクルの利点について説明する。
【００４０】
残差切除法における内部ソルバの役割は、反復過去の解修正ベクトル列とは独立な方向を見つけることであり、内部ソルバによって残差方程式が正確に解かれる必要はない。極端な表現をすれば、残差切除法の内部ソルバは、単体では連立一次方程式を解く能力がなくてもよい。
【００４１】
内部ソルバにＡＭＧ法を用いた場合、数値実験の結果から次のようなことが分かった。Ｖ−サイクルよりも手順の多いサイクル（Ｗ−サイクルなど）を用いたときは、いずれのサイクルを用いても残差切除法の反復１回あたりの収束率に大きな変化はない。一方、鋸歯−サイクルを用いたときは、収束率が顕著に低下する。Ｖ−サイクルでは、反復過去の解修正ベクトル列とは独立な方向を見つけることができるが、鋸歯−サイクルでは見つけることが困難である場合が多いと考えられる。ν−サイクルは、Ｖ−サイクルよりも鋸歯−サイクルに近いが、鋸歯−サイクルよりも計算量は少なくて済む。これは、鋸歯−サイクルでは最密グリッドから一段粗いグリッドに「制限」する前に行われる緩和が、ν−サイクルでは行われないからである。
【００４２】
そして、次のような理由によって、鋸歯−サイクルよりも計算量が少ないν−サイクルを用いた場合に、良好な結果が得られると考えられる。すなわち、残差切除法の内部ソルバによって計算される問題は初期値に零ベクトルを用いることが多いため、グリッドレベルが浅いときは、「制限」前に行う緩和の影響は小さいと考えられる。ν−サイクルでは、緩和の影響が小さいと考えられる浅いグリッドレベルでは、「制限」前に緩和を行わないか、または、緩和の反復回数を少なくしている。一方、深いグリッドレベルでは、「制限」前に行う緩和の反復回数を多くしている。また、グリッドレベルを昇る際の「補間」後の緩和は、Ｖ−サイクルと同様に、反復回数を一定回数にしている。したがって、ν−サイクルは、Ｖ−サイクルに近い結果を生むと考えられる。
【００４３】
＜リセット時のＡＭＧサイクル変更＞
ｍ反復目における相対残差εm （＝‖ｒm ‖／‖ｂ‖）に対して、残差切除率δm ＝（εm −εm+1）／εm が所定の値Δよりも小さくなったとき、次の処理を行う。
【００４４】
処理Ａ  ｍ反復目までの反復過去の解修正ベクトルの列φl （ｌ＝ｍ−Ｌmax＋１，…，ｍ）を破棄する。
【００４５】
処理Ｂ  内部ソルバのアルゴリズムを変更する。
【００４６】
処理Ａは、一般に反復法のリセットとして知られる手法である。反復法が局所解に陥った場合に対する対処法として知られており、所定の反復回数目にリセットを行う手法や、所定の反復回数毎にリセットを行う手法が報告されている。
【００４７】
残差切除法では、内部ソルバによって、解を修正すべき方向（ベクトル）を検知し、近似解が解修正ベクトルと反復過去の解修正ベクトルとの線形結合で表されるとの仮定の基に、次の反復における相対残差が最小になるように、線形結合の係数を決定している。ここで、上述の残差切除率が極めて小さくなるということは、利用している内部ソルバに、それ以上解を修正すべき方向を見い出す能力がない、と考えられる。このため、内部ソルバのアルゴリズムを変更することによって、それまでとは異なる解修正方向を見いだす可能性を高める。
【００４８】
そして、本実施形態では、アルゴリズムの変更の一例として、ＡＭＧサイクルにおける緩和の反復回数ｎの設定を変更するものとする。
【００４９】
具体的には、例えば次のように変更を行う。すなわち、計算開始から第１回目のリセットが発生するまでは、図５のようなν−サイクルを用いたＡＭＧアルゴリズムを用いる。第１回目のリセットの後は、「制限」前および「補間」後に行われる緩和の反復回数を、全てｎ＝ｎu 回とする。第２回目以降のリセットの後は、それまでの反復回数ｎに自然数ｎs を加え、これを新たな反復回数ｎとする。ｎが所定値Ｎに達したときは、未収束として、全ての処理を終了する。
【００５０】
本願発明者が考案した、リセット時にＡＭＧサイクルを変更する手法の利点は、次のように考えられる。上述の特許文献１でも論じられているように、残差切除法において内部ソルバにおける緩和回数を決定することは、困難かつ重要な問題である。緩和回数を小さくすると、反復１回当たりの計算量は小さくなるが、収束率が悪化するため、収束解に辿り着けない可能性が高くなる。逆に、緩和回数を大きくすると、反復１回当たりの計算量は大きくなるが、収束率が高くなり、収束解を安定的に得ることができる。
【００５１】
上述したＡＭＧサイクルの変更では、当初は計算量が少ないν−サイクルを用いて計算し、収束率（残差切除率）が小さくなったとき、計算量は多いが解修正方向を見出す能力が高いと考えられるＶ−サイクルに変更し、さらに収束率が小さくなったときは、Ｖ−サイクルの各段の緩和回数を増やしていく。これによって、できるだけ少ない計算量で、かつ、安定的に収束解を得ることができる。
【００５２】
図３に示すフローチャートに従って、本実施形態に係る数値計算方法について説明する。
【００５３】
いま求めたい物理量をＵとし、その物理量Ｕが満たすべき偏微分方程式を離散化したときの係数行列をＡ（Ｎ行Ｎ列）、非斉次項（ソース項）をｆとすると、解くべき式は、
Ａ・Ｕ＝ｆ  …（１）
で表される。この式は、一般的には、ＳＯＲ法、ＡＤＩ法等の逐次近似法によって解くことができる。
【００５４】
式（１）の解Ｕ∞を、近似解Ｕと摂動量（真の解との差）φとによって、次のように表す。
Ｕ∞＝Ｕ＋φ  …（２）
そして、摂動量φを求めて近似解Ｕを修正していくことによって、式（１）の真の解Ｕ∞を求める。
【００５５】
ここで、近似解Ｕに対する残差ｒを次のように定義する。
ｒ＝ｆ−Ａ・Ｕ  …（３）
式（１）〜（３）から、
Ａ・（Ｕ＋φ）＝ｆ
∴  Ａ・φ＝ｆ−Ａ・Ｕ
＝ｒ
したがって、摂動量φを求めるためには、次の式を解けばよい。
Ａ・φ＝ｒ  …（４）
【００５６】
本実施形態では、式（４）の収束解を求めるのではなく、最も収束勾配の急な最小単位の反復で予測近似値Ψを求め、求めた予測近似値Ψを、最適化制御ルーチンに供給する。そして、最適化制御ルーチンの実行結果が所定の条件を満たすまで、予測近似値Ψの算出を繰り返し実行する。なお、処理の詳細については、特許文献１に記載されている。
【００５７】
まず、ステップＳ１において、対象物理量Ｕの初期値Ｕ0を設定する。そして、ステップＳ２において、ＡＭＧセットアップを行う。すなわち、疎グリッドの生成、疎密グリッド間の「制限」「補間」関数の生成、および疎グリッドにおける係数行列の生成が、行われる。
【００５８】
次に、ステップＳ３において、繰り返し回数ｍに初期値として０を、摂動量φの初期値として０を与え、残差ｒの初期値ｒ0として（ｆ−Ａ・Ｕ0）を設定する。また、相対残差εの初期値ε0として、‖ｒ‖／‖ｆ‖を与える。さらに、ＡＭＧサイクルにおける緩和の反復回数に関して、ｎ＝０とし、またｎu，ｎsを所定値に設定する。また、リセットを判断するための所定値Δを設定する。
【００５９】
以下、近似解Ｕmが収束するまで、繰り返し回数ｍをインクリメントしながら、ステップＳ４〜Ｓ８を繰り返し実行する（Ｓ９，Ｓ１１）。また、残差切除率δmが所定値Δよりも小さくなったとき、修正近似解φm を破棄し、内部ソルバのアルゴリズムを変更する。ここでは、ＡＭＧサイクルにおける緩和の反復回数ｎの設定の変更を行うものとする。
【００６０】
ステップＳ４〜Ｓ７（第１のステップ）では、ＡＭＧ法を採用した内部ソルバ１１を有する第１演算部１０によって、
Ａ・φ＝ｒm
の予測近似値Ψm を求める。計算当初から第１回目のリセットがかかるまでは、ｎ＝０に設定されているので（ステップＳ４でＮｏ）、ステップＳ５において、上述したν−サイクルをＡＭＧサイクルとして実行する。すなわち、緩和法の反復回数ｎが、
１）レベルｋからレベルｋ＋１に「制限」する前に行う緩和については、
ｎ＝ｍａｘ（ｋ−ｎb ，０）
２）レベルｋ＋１からレベルｋに「補間」する前に行う緩和については、
ｎ＝ｎu
と設定されている（ｎbは自然数）。
【００６１】
ステップＳ４〜Ｓ７において予測近似値Ψm が求められると、ステップＳ８において、この予測近似値Ψm を用いて、第２演算部２０における最適化ルーチンによって、次回の残差ｒm+1のＬ２ノルムが最小となるように、最適化された修正近似値φmを求める。残差最小化係数αl（ｌ＝１，…，Ｌ）の計算など最適化ルーチンにおける処理については、特許文献１で詳述されており、ここでは省略する。
【００６２】
そして、この最適化された修正近似値φmによって、第（ｍ＋１）回目の繰り返しにおける近似値Ｕm+1および残差ｒm+1は、次式で求められる。
【００６３】
Ｕm+1＝Ｕm＋φm
ｒm+1＝ｒm−Ａ・φm
また、相対残差εm+1と、残差切除率δmもそれぞれ求められる。
【００６４】
そして、近似解Ｕmがまだ収束しないとき（ステップＳ９でＮｏ）、残差切除率δmを所定値Δと比較する（ステップＳ１０）。そして、残差切除率δmが所定値Δを上回っているときは（Ｓ１０でＹｅｓ）、繰り返し回数ｍをインクリメントし（Ｓ１１）、ステップＳ４に戻る。
【００６５】
一方、残差切除率δmが所定値Δを下回ったときは（Ｓ１０でＮｏ）、リセット処理を行い、内部ソルバのアルゴリズムを変更する。
【００６６】
第１回目のリセットのときは、ｎ＝０であるので（ステップＳ１２でＮｏ）、ｎ＝ｎuとする。また、反復過去の解修正ベクトルφl（ｌ＝ｍ−Ｌmax＋１，…，ｍ）を零とする（ステップＳ１３）。これにより、以降の繰り返しでは、ステップＳ４の分岐においてステップＳ６にすすみ、内部ソルバにおいて、Ｖ−サイクルがＡＭＧサイクルとして実行される。Ｖ−サイクルの各段における緩和の反復回数は、ｎ＝ｎuである。
【００６７】
第２回目以降のリセットのときは、ｎ＞０であるので（ステップＳ１２でＹｅｓ）、ｎuにｎsが加えられる（ステップＳ１４）。そして、ｎ＝ｎuとする。また、反復過去の解修正ベクトルφl（ｌ＝ｍ−Ｌmax＋１，…，ｍ）を零とする（ステップＳ１３）。これにより、以降の繰り返しでは、内部ソルバにおいて、緩和の反復回数がｎsだけ増えた，Ｖ−サイクルによるＡＭＧ法が実行される。
【００６８】
なお、ここでは、リセット時の内部ソルバアルゴリズムの変更を、ＡＭＧサイクルにおける緩和の反復回数の設定変更によって行うものとしたが、本発明はこれに限られるものではなく、例えば、あるアルゴリズムから全く別のアルゴリズムに変更することも考えられる。また、緩和の反復回数の変更の仕方も、ここで示したものに限られるものではなく、例えば、当初からＶ−サイクルを採用してリセット毎に緩和回数を増やすようにしてもよい。
【００６９】
また、リセット処理を行わない場合であっても、ν−サイクルを採用することによって、少ない計算量によって、Ｖ−サイクルとほぼ同等の性能を得ることができる。
【００７０】
＜数値実験＞
本実施形態で示した数値計算方法を用いて、コンピュータによって、流体解析によって得られた連立一次方程式を実際に解いた。未知変数は、速度３成分と圧力の計４種類であり、その個数は合計１４３，３８６元である。
【００７１】
図６は数値実験の結果得られた収束曲線を示すグラフである。収束判定条件は、相対残差が１．０ｅ−５以下となることとした。緩和法にはＳＯＲ法（緩和係数ω＝０．２）を用い、ｎb＝２，ｎu＝２，ｎs＝２とした。また比較例として、ＡＭＧサイクル内における緩和法の反復回数を４回に固定した場合を、破線で示している。
【００７２】
図６から分かるように、本実施形態によると、比較例よりも速い段階から収束を始め、収束が一旦止まった後も、３回目のリセットの後に、再び収束を始めている。これに対して比較例では、本実施形態よりも収束速度は遅く、また、一旦は本実施形態よりも精度が上回ったものの、やがて収束が止まってしまっている。したがって、本実施形態の方が、短時間での収束性が高く、また最終的に得られる精度も比較例を大きく上回っている。
【産業上の利用可能性】
【００７３】
本発明では、ＡＭＧ法を内部ソルバに採用した残差切除法を用いた数値計算において、収束性を高めることができるので、例えば、構造解析や流体解析などの数値計算の精度向上や速度向上に、有効である。
【図面の簡単な説明】
【００７４】
【図１】ＡＭＧサイクルの一例としてのＶ−サイクルを概念的に示す図である。
【図２】ＡＭＧサイクルの一例を概念的に示す図であり、（ａ）はＷ−サイクル、（ｂ）は鋸歯−サイクルである。
【図３】本発明の一実施形態に係る数値計算方法を示すフローチャートである。
【図４】本発明の一実施形態に係る数値計算装置の構成を概念的に示すブロック図である。
【図５】本実施形態において用いるＡＭＧサイクルの一例であるν−サイクルを概念的に示す図である。
【図６】数値実験の結果得られた収束曲線を示すグラフである。
【符号の説明】
【００７５】
１  数値計算装置
１０  第１演算部
１１  内部ソルバ
２０  第２演算部
【図１】
2005267418.tif 000003
【図２】
2005267418.tif 000004
【図３】
2005267418.tif 000005
【図４】
2005267418.tif 000006
【図５】
2005267418.tif 000007
【図６】
2005267418.tif 000008