明細書 :数値計算方法および数値計算装置、並びに数値計算用プログラムを記録した記録媒体

発行国	日本国特許庁(JP)
公報種別	特許公報(B2)
特許番号	特許第3394989号 (P3394989)
公開番号	特開2001-134304 (P2001-134304A)
登録日	平成15年2月7日(2003.2.7)
発行日	平成15年4月7日(2003.4.7)
公開日	平成13年5月18日(2001.5.18)
発明の名称または考案の名称	数値計算方法および数値計算装置、並びに数値計算用プログラムを記録した記録媒体
国際特許分類	G05B 13/02       G06F 17/00       G06F 17/13       G08G  1/16
ＦＩ	G05B 13/02 K G08G 1/16 G06F 15/20 G06F 15/328
請求項の数または発明の数	6
全頁数	14
出願番号	特願平11-310770 (P1999-310770)
出願日	平成11年11月1日(1999.11.1)
審査請求日	平成12年10月5日(2000.10.5)
特許権者または実用新案権者	【識別番号】501137577 【氏名又は名称】独立行政法人航空宇宙技術研究所 【識別番号】599154087 【氏名又は名称】株式会社ヴァイナス
発明者または考案者	【氏名】菊地 一雄 【氏名】高橋 匡康 【氏名】田村 敦宏 【氏名】谷口 幸二
個別代理人の代理人	【識別番号】100077931、【弁理士】、【氏名又は名称】前田 弘 （外４名）
審査官	【審査官】森林 克郎
参考文献・文献	特開 平６－９５７０７（ＪＰ，Ａ） 特開 平８－５４５０（ＪＰ，Ａ）
調査した分野	G05B 11/00 - 13/04 G06F 15/20,15/328
特許請求の範囲	【請求項1】 物理量Ｕが満たすべき偏微分方程式を離散化したときの係数行列をＡ、非斉次項（ソース項）をｆとしたとき、コンピュータを用いて、 Ａ・Ｕ＝ｆ を解き、物理量Ｕの数値計算を行う方法であって、 物理量Ｕの初期値Ｕ0を設定し、 繰り返し回数ｍに初期値として０を、摂動量φの初期値として０を与え、残差ｒの初期値ｒ0として（ｆ－Ａ・Ｕ0）を設定し、 繰り返し回数ｍをインクリメントしながら、 Ａ・φ＝ｒmの予測近似値Ψmを、内部ソルバを有する第１演算部によって、反復計算によって求めるステップと、 残差ｒmのＬ２ノルムを最小とする最適化された近似値φmを、第２演算部によって、最適化ルーチンによって、予測近似値Ψmから求め、近似解Ｕm+1として（Ｕm＋φm）を与え、残差ｒm+1として（ｒm－Ａ・φm）を与えるステップとを、近似解Ｕmが収束するまで、繰り返し実行するものであり、 前記予測近似値Ψmを求めるステップは、 繰り返し回数ｍ毎の残差切除率を表す変数κmが所定値Ｋを超えたとき、反復計算を終えるものであることを特徴とする数値計算方法。 【請求項2】 請求項１記載の数値計算方法において、 前記変数κmは、 【数１】 の式で与えられることを特徴とする数値計算方法。 【請求項3】 物理量Ｕが満たすべき偏微分方程式を離散化したときの係数行列をＡ、非斉次項（ソース項）をｆとしたとき、 Ａ・Ｕ＝ｆ を解き、物理量Ｕの数値計算を行う装置であって、 物理量Ｕの初期値Ｕ0を設定し、 繰り返し回数ｍに初期値として０を、摂動量φの初期値として０を与え、残差ｒの初期値ｒ0として（ｆ－Ａ・Ｕ0）を設定し、 繰り返し回数ｍをインクリメントしながら、 Ａ・φ＝ｒmの予測近似値Ψmを、内部ソルバを有する第１演算部によって、反復計算によって求めるステップと、 残差ｒmのＬ２ノルムを最小とする最適化された近似値φmを、第２演算部によって、最適化ルーチンによって、予測近似値Ψmから求め、近似解Ｕm+1として（Ｕm＋φm）を与え、残差ｒm+1として（ｒm－Ａ・φm）を与えるステップとを、近似解Ｕmが収束するまで、繰り返し実行するものであり、 前記予測近似値Ψmを求めるステップは、 繰り返し回数ｍ毎の残差切除率を表す変数κmが所定値Ｋを超えたとき、反復計算を終えるものであることを特徴とする数値計算装置。 【請求項4】 請求項３記載の数値計算装置において、 前記変数κmは、 【数１】 の式で与えられることを特徴とする数値計算装置。 【請求項5】 コンピュータに、物理量Ｕが満たすべき偏微分方程式を離散化したときの係数行列をＡ、非斉次項（ソース項）をｆとしたとき、 Ａ・Ｕ＝ｆ を解かせ、物理量Ｕの数値計算を行う数値計算用プログラムを記録した記録媒体であって、 物理量Ｕの初期値Ｕ0を設定し、 繰り返し回数ｍに初期値として０を、摂動量φの初期値として０を与え、残差ｒの初期値ｒ0として（ｆ－Ａ・Ｕ0）を設定し、 繰り返し回数ｍをインクリメントしながら、 Ａ・φ＝ｒmの予測近似値Ψmを、内部ソルバを有する第１演算部によって、反復計算によって求めるステップと、 残差ｒmのＬ２ノルムを最小とする最適化された近似値φmを、第２演算部によって、最適化ルーチンによって、予測近似値Ψmから求め、近似解Ｕm+1として（Ｕm＋φm）を与え、残差ｒm+1として（ｒm－Ａ・φm）を与えるステップとを、近似解Ｕmが収束するまで、繰り返し実行するものであり、 前記予測近似値Ψmを求めるステップは、 繰り返し回数ｍ毎の残差切除率を表す変数κmが所定値Ｋを超えたとき、反復計算を終えるものである処理をコンピュータに実行させる数値計算用プログラムを記録した記録媒体。 【請求項6】 請求項５記載の数値計算用プログラムを記録した記録媒体において、 前記変数κmは、 【数１】 の式で与えられることを特徴とする数値計算用プログラムを記録した記録媒体。
発明の詳細な説明	【発明の詳細な説明】 【0001】 【発明の属する技術分野】本発明は、圧力、温度などの物理量を対象とした数値計算技術に関するものであり、特に、逐次近似解法アルゴリズムを用いた数値計算技術に属する。 【0002】 【発明が解決しようとする課題】従来から、圧力、温度などの物理量を対象とした制御技術において、逐次近似解法アルゴリズムの利用がなされている。逐次近似解法アルゴリズムとは、様々な最適制御理論に基づいて、逐次近似を繰り返し行いつつ、数値計算により最適解を求めるものである。 【0003】 ところが、逐次近似解法アルゴリズムを用いた場合には、制御対象の各パラメータの解を求める際に、計算の反復回数が多くなるために、計算時間の点から、高速かつ高精度の制御が困難になっている。このため、様々な分野の制御技術に応用できるようにするために、逐次近似解法アルゴリズムの収束性の向上が、望まれている。 【0004】 前記の問題に鑑み、本発明は、従来よりも収束性の高い逐次近似解法アルゴリズムを用いた数値計算方法、数値計算装置および数値計算用プログラムを記録した記録媒体を提供することを課題とする。 【0005】 【課題を解決するための手段】前記の課題を解決するために、請求項１の発明が講じた解決手段は、物理量Ｕが満たすべき偏微分方程式を離散化したときの係数行列をＡ、非斉次項（ソース項）をｆとしたとき、コンピュータを用いて、Ａ・Ｕ＝ｆを解き、物理量Ｕの数値計算を行う方法として、物理量Ｕの初期値Ｕ0を設定し、繰り返し回数ｍに初期値として０を、摂動量φの初期値として０を与え、残差ｒの初期値ｒ0として（ｆ－Ａ・Ｕ0）を設定し、繰り返し回数ｍをインクリメントしながら、Ａ・φ＝ｒmの予測近似値Ψmを内部ソルバを有する第１の演算部によって反復計算によって求めるステップと、残差ｒmのＬ２ノルムを最小とする最適化された近似値φmを第２演算部によって最適化ルーチンによって予測近似値Ψmから求め、近似解Ｕm+1として（Ｕm＋φm）を与え、残差ｒm+1として（ｒm－Ａ・φm）を与えるステップとを近似解Ｕmが収束するまで繰り返し実行するものであり、前記予測近似値Ψmを求めるステップは、繰り返し回数ｍ毎の残差切除率を表す変数κmが所定値Ｋを超えたとき、反復計算を終えるものである。 【0006】 そして、請求項２の発明では、前記請求項１の数値計算方法における変数κmは、 【数１】 の式で与えられるものとする。 【0007】 また、請求項３の発明が講じた解決手段は、物理量Ｕが満たすべき偏微分方程式を離散化したときの係数行列をＡ、非斉次項（ソース項）をｆとしたとき、Ａ・Ｕ＝ｆを解き、物理量Ｕの数値計算を行う装置として、物理量Ｕの初期値Ｕ0を設定し、繰り返し回数ｍに初期値として０を、摂動量φの初期値として０を与え、残差ｒの初期値ｒ0として（ｆ－Ａ・Ｕ0）を設定し、繰り返し回数ｍをインクリメントしながら、Ａ・φ＝ｒmの予測近似値Ψmを内部ソルバを有する第１演算部によって反復計算によって求めるステップと、残差ｒmのＬ２ノルムを最小とする最適化された近似値φmを第２演算部によって最適化ルーチンによって予測近似値Ψmから求め、近似解Ｕm+1として（Ｕm＋φm）を与え、残差ｒm+1として（ｒm－Ａ・φm）を与えるステップとを近似解Ｕmが収束するまで、繰り返し実行するものであり、前記予測近似値Ψmを求めるステップは、繰り返し回数ｍ毎の残差切除率を表す変数κmが所定値Ｋを超えたとき反復計算を終えるものである。 【0008】 そして、請求項４の発明では、前記請求項３の数値計算装置における変数κmは、 【数１】 の式で与えられるものとする。 【0009】 また、請求項５の発明が講じた解決手段は、コンピュータに、物理量Ｕが満たすべき偏微分方程式を離散化したときの係数行列をＡ、非斉次項（ソース項）をｆとしたとき、Ａ・Ｕ＝ｆを解かせ、物理量Ｕの数値計算を行う数値計算用プログラムを記録した記録媒体として、物理量Ｕの初期値Ｕ0を設定し、繰り返し回数ｍに初期値として０を、摂動量φの初期値として０を与え、残差ｒの初期値ｒ0として（ｆ－Ａ・Ｕ0）を設定し、繰り返し回数ｍをインクリメントしながら、Ａ・φ＝ｒmの予測近似値Ψmを内部ソルバを有する第１演算部によって反復計算によって求めるステップと、残差ｒmのＬ２ノルムを最小とする最適化された近似値φmを第２演算部によって最適化ルーチンによって予測近似値Ψmから求め、近似解Ｕm+1として（Ｕm＋φm）を与え、残差ｒm+1として（ｒm－Ａ・φm）を与えるステップとを近似解Ｕmが収束するまで繰り返し実行するものであり、前記予測近似値Ψmを求めるステップは、繰り返し回数ｍ毎の残差切除率を表す変数κmが所定値Ｋを超えたとき反復計算を終えるものである処理をコンピュータに実行させる数値計算用プログラムを記録したものである。 【0010】 請求項６の発明では、前記請求項５の数値計算用プログラムを記録した記録媒体において、前記数値計算用プログラムで用いる変数κmは、 【数１】 の式で与えられるものとする。 【0011】 【発明の実施の形態】以下、本発明の一実施形態について、図面を参照して説明する。 【0012】 図１は本発明の一実施形態に係る制御装置の構成を示すブロック図であり、本発明に係る数値計算技術を車両の自動走行制御に応用した例を示している。図１において、車両は、目的値・経路設定部１１によって設定される軌道上を許容誤差範囲内で目的地まで自動走行する回路を備えているものとする。走行中は、位置情報検出器１２によって、位置・方位情報や速度等の走行情報が検出される。表示器１８では設定値や走行状態などの情報がリアルタイムでモニター表示される。 【0013】 車両の位置と軌道とのずれφは情報管理装置１３によって常に監視されている。そして、このずれφが許容値δmax を越えたとき、そのときの走行情報を出発値として、方向角・角加速度・速度・加速度などの最適化パラメータ制御、すなわち軌道復帰への最適化手段が、本発明に係るアルゴリズムを用いてシミュレーションされる。 【0014】 図２は本発明の一実施形態に係る数値計算方法を示すフローチャートである。 【0015】 （残差方程式の定式化） いま求めたい物理量をＵとし、その物理量が満たすべき偏微分方程式を離散化したときの係数行列をＡ、非斉次項（ソース項）をｆとすると、解くべき式は、 Ａ・Ｕ＝ｆ …（１） で表される。この式は、一般的には、ＳＯＲ法、ＡＤＩ法等の逐次近似法によって解くことができる。 【0016】 本発明では、以下のような解法を用いる。 【0017】 式（１）の解Ｕ∞を、近似解Ｕと摂動量（真の解との差）φとによって、次のように表す。 Ｕ∞＝Ｕ＋φ …（２） 本発明に係る解法では、摂動量φを求めて近似解Ｕを修正していくことによって、式（１）の真の解Ｕ∞を求める。 【0018】 ここで、近似解Ｕに対する残差ｒを次のように定義する。 【0019】 ｒ＝ｆ－Ａ・Ｕ …（３） 式（１）～（３）から、 Ａ・（Ｕ＋φ）＝ｆ ∴ Ａ・φ＝ｆ－Ａ・Ｕ＝ｒ したがって、摂動量φを求めるためには、次の式を解けばよい。 Ａ・φ＝ｒ …（４） 【0020】 （残差切除法のアルゴリズム） 本発明に係るアルゴリズムでは、式（４）の収束解を求めるのではなく、ＡＤＩ法などによって最も収束勾配の急な最小単位の反復で予測近似値Ψを求め、求めた予測近似値Ψを、最適化制御ルーチンに供給する。そして、最適化制御ルーチンの実行結果が所定の条件を満たすまで、予測近似値Ψの算出を繰り返し実行する。 【0021】 ＜最適化制御ルーチン＞ いま、繰り返し回数をｍとしたとき、残差のＬ２ノルムを最小とする合成摂動量φmと新しい近似解Ｕm+1を次のように定義する。 【数２】 Ｕm+1＝Ｕm＋φm …（６） αl（ｌ＝１，２，３，…，Ｌ）は残差最小化係数であり、後述する計算方法によって求まる定数である。残差のＬ２ノルムがより小さくなるように、次のように残差の最小化を行う。 【0022】 新しい近似解Ｕm+1が式（６）で表されるとき、これに対する残差ｒm+1は、式（３）から、次のように表される。 【0023】 ｒm+1＝ｆ－Ａ・Ｕm+1 …（７） 式（７）に式（６）を代入し、式（３），（５）を用いると、 【数３】 【００２４】近似解Ｕm+1に対する残差ｒm+1のＬ２ノルム∥ｒm+1∥は、三次元の場合、（ｒm+1）2の内点全ての総和の平方根として、式（９）で与えられる。 【数４】 【００２５】この残差ｒm+1のＬ２ノルム∥ｒm+1∥を最小にするように、式（９）の残差最小化係数αl（ｌ＝１，２，３，…，Ｌ）を最小二乗法で定める。すなわち、 【数５】 式（１０）はαのＬ元連立方程式となり、これを数値的に解くことによってαlが定まる。αlが定まると、式（５）からφmが求まり、式（６）からＵm+1が求まる。残差最小化係数αlの計算方法については後述する。 【0026】 ｍをインクリメントしながら、同様のルーチンを繰り返すことによって、式（９）の残差のＬ２ノルムを零または最小にする解Ｕ∞に収束する。 【0027】 図２に示すフローチャートに従って、本実施形態に係る数値計算方法について説明する。 【0028】 まず、ステップＳ１において、対象物理量Ｕの初期値Ｕ0を設定する。そして、ステップＳ２において、繰り返し回数ｍに初期値として０を、摂動量φの初期値として０を与え、残差ｒの初期値ｒ0として（ｆ－Ａ・Ｕ0）を設定する。 【0029】 以下、近似解Ｕmが収束するまで、繰り返し回数ｍをインクリメントしながら、ステップＳ３～Ｓ９を繰り返し実行する。 【0030】 ステップＳ３～Ｓ８では、ＡＤＩ法、ＳＯＲ法、ＣＧ法などの逐次近似法を実行する内部ソルバによって、 Ａ・φ＝ｒm の近似解（予測近似値）Ψm を求める。ところがこの式は、既存の逐次近似法で解こうとすると、実用問題に対しては収束までに多くの繰り返しを必要とし、必然的に計算機使用時間も膨大となる。 【0031】 そこで本実施形態では、逐次計算の反復回数の上限Ｎを設定する（ステップＳ７，Ｓ８）とともに、残差の減少の割合すなわち残差切除率を表す変数κmを導入し（ステップＳ５）、この変数κmが所定値Ｋを超えたとき、反復計算を終えるものとしている（ステップＳ６）。 【0032】 逐次計算の反復回数を固定した場合には、収束判定までのループの繰り返しｍ毎に、残差の減少の割合（残差切除率）が大きく変動する。この変動は、各繰り返しｍにおける予測近似値Ψm の近似度合に起因するものである。このために、解が収束するまでのループの繰り返し回数が必要以上に増加してしまうおそれがある。 【0033】 これに対して、本実施形態では、予測近似値Ψm を求めるステップにおいて、残差切除率を表す変数κm が所定値Ｋを超えたときに反復計算を終えるものとしているので、ループの繰り返しｍ毎の残差切除率を一定値以上に保つことができ、したがって、解が収束するまでに要するループの繰り返し回数を減少させることができる。 【0034】 ここでは、残差切除率を表す変数κmを次のように定義する。 【数６】 【００３５】この変数κmは、反復計算を繰り返すにつれて、予測近似値Ψmの精度が高くなるために、１．０に限りなく近づく。したがって、この変数κmの値が任意に設定した所定値Ｋを超えるまで反復計算を行うことによって、ループの繰り返し回数ｍ毎の残差切除率はほぼ一定になる。 【0036】 ステップＳ３～Ｓ８において予測近似値Ψ0 が求められると、ステップＳ９において、この予測近似値Ψ0 を用いて、最適化ルーチンによって、次回の残差ｒ1が最小となるように、最適化された近似値φ0を求める。この最適化された近似値φ0によって、第１次の物理量の近似値Ｕ1および第１次の残差ｒ1は、式（２），（３）から、次式で求められる。 Ｕ1＝Ｕ0＋φ0 …（１２） ｒ1＝ｆ－Ａ・Ｕ1 …（１３） 式（１３）に式（１２）を代入して、 ｒ1＝ｆ－Ａ・（Ｕ0＋φ0） ＝ｆ－Ａ・Ｕ0－Ａ・φ0 ＝ｒ0－Ａ・φ0 …（１４） よって、式（１４）により、第１次の残差ｒ1が求められる。 【0037】 ステップＳ１１でｍをインクリメントしてステップＳ３に戻り、同様の計算を行う。このような処理を繰り返すことによって、 （Ｕ2 ，ｒ2 ），（Ｕ3 ，ｒ3 ）…（Ｕm ，ｒm ）… が順次求められ、残差ノルム∥ｒm∥が減少していく。 【0038】 繰り返し回数がｍのとき、ステップＳ４において予測近似値Ψmを求める方程式は、 Ａ・φ＝ｒm となり、ステップＳ９では、最適化制御ルーチンによって最適化された近似値φmによって、 Ｕm+1＝Ｕm＋φm ｒm+1＝ｒm－Ａ・φm と求められる。 【0039】 なお、∥ｒm∥＝０となったときが解であるが、境界条件によっては、最終残差ｒfinal が残る場合がある。すなわち、境界条件が全て物理量の微分で与えられる場合（ノイマン問題）、数学的にソース項ｆと境界条件との間で満たすべき関係（ノイマンの拘束条件）があり、これが満たされないために、残差ノルムが ∥ｒm∥＝∥ｒm-1∥＝∥ｒm-2∥… となってしまう。 【0040】 他の解法では、このようなノイマン問題の場合、そのままでは収束しないため、最終残差ｒfinal によって、ソース項ｆあるいは境界条件を修正して収束解を得ようとする。これは、問題を変更していることになる。ところが、本実施形態に係る解法では、このような修正を行わないでも自動的に最終残差ｒfinal が確定し、収束解が得られる特性を持っている。 【0041】 （残差最小化係数αlの求め方） いま、表記を簡単にするために、次のようにδを導入する。 δm ＝Ψm δm-l+1＝φm-l+1 （ｌ＝２，３，…，Ｌ） とおくと、最適化された近似値φmは、 【数７】 と表される。よって、 【数８】 Σ記号を外して表記すると、 ｒm+1＝ｒm－（α1Ａ・δm＋α2Ａ・δm-1＋α3Ａ・δm-2＋…＋αLＡ・δm-L+1） となる。 【0042】 残差最小化係数αlは、例えば最小自乗法を用いて、（ｍ＋１）次の残差ノルム∥ｒm+1∥（自乗和の平方根）を最小にするという条件を与えることによって、求められる。 【数９】 これをαlで微分すると、 【数１０】 【００４３】この式をαlについて解くことによって、残差最小化係数αlを求めることができる。微分を実行すれば、 【数１１】 となる。ここで、ｌ'＝１，２，…，Ｌである。この連立一次方程式を解くことによって、残差最小化係数αlを求めることができる。 【0044】 例えば、Ｌ＝３のときは、上式は、 【数１２】 のような３元連立一次方程式になる。そして、 φm ＝α1δm＋α2δm-1＋α3δm-2 ｒm+1＝ｒm－（α1Ａ・δm＋α2Ａ・δm-1＋α3Ａ・δm-2） ＝ｒm－Ａ・φm Ｕm+1＝Ｕm＋φm となる。 【0045】 図１に示す車両の走行制御においては、まず予測修正量演算部１４によって、現在の走行情報を基にして予測修正量Ψm が反復計算される。ここでは、繰り返し回数毎の残差切除率κm が所定値Ｋを越えるまで反復計算される。次に、最適加算修正量演算部１５によってｍ次の加算修正量φm が演算される。予測修正量演算部１４および最適加算修正量演算部１５によって求められた解は、最終的に修正量収束判定部１６によって収束が判定され、所定の軌道に復帰するための最適制御方法が決定される。決定された制御情報は方向制御部１７に送られ、車両の走行経路が修正される。 【0046】 なお、ここでは車両の走行制御を例にとって本発明について説明したが、他の分野の様々な制御にも本発明は利用可能であることはいうまでもない。 【0047】 また、上記の数値計算方法は、当該方法を実現するためのプログラムを実行するコンピュータを備えた装置によって容易に実現することができる。また、当該方法を実現するためのプログラムをコンピュータ読み取り可能な記録媒体に記録して、この記録媒体に記録したプログラムをコンピュータに実行させることによって実現することができる。 【0048】 ここで、従来技術との比較について補足説明する。従来、機器の制御に用いられてきた技術には、比例制御やファジィ制御等がある。位置制御を例にとると、比例制御の場合には、基本的には現在位置と設定位置との比較を２値化して判断し、そのずれを限りなく零に近づけようとするため、微調整動作が頻繁に行われる。このため、動作が安定しない、また外乱に対して柔軟に対処できないという欠点があった。 【0049】 このような欠点を解消する技術として、ファジィ理論を採用した制御技術が開発された。これは、その状態量を「大きくずれている」「少しずれている」「微妙にずれている」などと「あいまいに」判断し、これらの状態に対して「大きく修正する」「少し修正する」「そのまま」のような「あいまいな」ルールに基づいて制御を行うものである。この技術によって、ＡＩ制御に近い「滑らかな制御」が実現されたが、情報量と制御パラメータに比例してルール数が増加するため、その設定とチューニング作業に手間がかかる、または主観的な設定が可能なために制御結果にばらつきが生じるという欠点もあった。 【0050】 本発明は、制御対象を目的の状態に制御するための最適手段を決定するために用いられる技術であり、その基となるのは制御パラメータに対応した非線形方程式である。このため、対象によっては目的の状態が数値的に明確でないケースへも適応可能であり、また最適な状態を自動的にチューニングする機能も備えている。 【0051】 （他の応用例その１） 図３は本発明に係るアルゴリズムの他の応用例を示す図である。図３では、数値計算への応用例として、２次元または３次元楕円型方程式法の格子生成手法への適用について示すフローである。 【0052】 図３で用いた符号の意味は下記のとおりである。 Ｕ ：物理面における格子点の座標値（ｘ，ｙ，ｚ） Ｐ ：格子生成の制御関数 Ａ ：線型化された係数行列 Ψ(m) ：ｍ次ステップでの予測修正値 φ(m) ：ｍ次ステップでの最適化された修正値 αｌ ：未定係数行列 ｒ(m) ：ｍ次ステップでの残差 Ｌ ：未定係数行列の個数 【0053】 ステップＳ２１では、ｍ次ステップにおける値に基づいて、離散化された楕円型方程式系の係数行列を線型化させる。ステップＳ２２では、κ技術によって繰り返し回数を制御して、（ｍ＋１）次ステップにおける予測修正量Ψ(m) を求める。ステップＳ２３では、（ｍ＋１）次ステップの残差を最小（ゼロ）にするための未定係数行列αｉを求める。ステップＳ２４では、予測修正値Ψ(m) 、（ｍ－１）次以前の最適化された修正値（φ(m-1) ，…，φ(m-L+1) ）および未定係数行列に基づいて（ｍ＋１）次ステップでの最適化された修正値φ(m+1) を求める。ステップＳ２５では、最適化された修正値φ(m+1) に基づいて、物理面における格子点の座標値を修正する。 【0054】 図４は本応用例によって生成された格子、図５は比較例としてのＳＯＲ法によって生成された格子である。図４および図５では、気流中におかれた楕円柱周りの２次元格子生成を行っている。格子点数は９０×９０である。本応用例では、、残差切除法の内部ソルバとしてＡＤＩ法を用いた。 【0055】 図６は収束速度を比較するためのグラフである。図６において、ａ１，ａ２は本応用例におけるＸ方向およびＹ方向の残差、ｂ１，ｂ２は比較例におけるＸ方向およびＹ方向の残差である。図６から分かるように、本発明に係るアルゴリズムを利用した場合には、ＳＯＲ法と比較して６倍以上高速に収束している。 【0056】 （他の応用例その２） 図７は本発明に係るアルゴリズムの他の応用例を示す図である。図７では、数値計算への応用例として、流体解析のＭＡＣ法の圧力計算部に残差切除法を適用した例を示すフローである。 【0057】 図７で用いた符号の意味は下記の通りである。 ｔ ：時刻 Ｕ ：流速 Ｐ ：圧力 Ψm ：ｍ次ステップでの予測修正量 φm ：ｍ次ステップでの最適化された修正量 α ：最適化係数 ｒm ：ｍ次ステップでの残差 【0058】 Ｓ３１は初期残差ｒ0 の算出を含む初期化部分、Ｓ３２は予測近似値Ψm の算出部、Ｓ３３はκm の算出部、Ｓ３４はＳ３３で算出されたκm によって繰り返し数を制御されたＳ３２およびＳ３３を含むループ、Ｓ３５は最適化係数α、最適化された修正量φm および次ステップでの残差ｒm+1 の算出部である。 【0059】 本応用例による効果を評価するために、ＭＡＣ法によって、一様流中に置かれた円柱周りの二次元解析を行った。比較例としてはＳＯＲ法を用いた。格子点数は６０×６０、主流速度は５ｃｍ／ｓである。本応用例では、残差切除法の内部ソルバとしてＳＯＲ法を用いた。 【0060】 図８はタイムステップ数とＣＰＵ時間との関係を示すグラフである。図８において、ａは本応用例による結果、ｂは比較例としてのＳＯＲ法による結果を示す。また図９はステップ数と相対残差／相対誤差との関係を示すグラフである。図９において、ａ１，ａ２は本応用例による結果、ｂ１，ｂ２は比較例としてのＳＯＲ法による結果を示す。図８および図９から、本応用例によると、ＳＯＲ法に比べて５倍に高速化されていることが分かる。 【0061】 【発明の効果】以上のように本発明によると、予測近似値Ψmを求めるステップにおいて、繰り返し回数ｍ毎の残差切除率を表す変数κmが所定値Ｋを超えたとき、反復計算を終えるので、ループの繰り返しｍ毎の残差切除率を一定値以上に保つことができ、これにより、収束に必要なｍの値を減少させることができる。したがって、逐次近似解法アルゴリズムの収束性が向上するので、例えば制御に用いた場合に、高速かつ高精度の制御が容易になる。
図面	【図1】 【図8】 【図4】 【図5】 【図2】 【図3】 【図6】 【図9】 【図7】