TOP > 国内特許検索 > 数値計算方法および数値計算装置、並びに数値計算用プログラムを記録した記録媒体 > 明細書

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

発行国 日本国特許庁(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      
FI 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、【弁理士】、【氏名又は名称】前田 弘 (外4名)
審査官 【審査官】森林 克郎
参考文献・文献 特開 平6-95707(JP,A)
特開 平8-5450(JP,A)
調査した分野 G05B 11/00 - 13/04
G06F 15/20,15/328
特許請求の範囲 【請求項1】
物理量Uが満たすべき偏微分方程式を離散化したときの係数行列をA、非斉次項(ソース項)をfとしたとき、コンピュータを用いて、
A・U=f
を解き、物理量Uの数値計算を行う方法であって、
物理量Uの初期値U0を設定し、
繰り返し回数mに初期値として0を、摂動量φの初期値として0を与え、残差rの初期値r0として(f-A・U0)を設定し、
繰り返し回数mをインクリメントしながら、
A・φ=rmの予測近似値Ψmを、内部ソルバを有する第1演算部によって、反復計算によって求めるステップと、
残差rmのL2ノルムを最小とする最適化された近似値φmを、第2演算部によって、最適化ルーチンによって、予測近似値Ψmから求め、近似解Um+1として(Um+φm)を与え、残差rm+1として(rm-A・φm)を与えるステップとを、近似解Umが収束するまで、繰り返し実行するものであり、
前記予測近似値Ψmを求めるステップは、
繰り返し回数m毎の残差切除率を表す変数κmが所定値Kを超えたとき、反復計算を終えるものであることを特徴とする数値計算方法。

【請求項2】
請求項1記載の数値計算方法において、
前記変数κmは、
【数1】
JP0003394989B2_000002t.gifの式で与えられることを特徴とする数値計算方法。

【請求項3】
物理量Uが満たすべき偏微分方程式を離散化したときの係数行列をA、非斉次項(ソース項)をfとしたとき、
A・U=f
を解き、物理量Uの数値計算を行う装置であって、
物理量Uの初期値U0を設定し、
繰り返し回数mに初期値として0を、摂動量φの初期値として0を与え、残差rの初期値r0として(f-A・U0)を設定し、
繰り返し回数mをインクリメントしながら、
A・φ=rmの予測近似値Ψmを、内部ソルバを有する第1演算部によって、反復計算によって求めるステップと、
残差rmのL2ノルムを最小とする最適化された近似値φmを、第2演算部によって、最適化ルーチンによって、予測近似値Ψmから求め、近似解Um+1として(Um+φm)を与え、残差rm+1として(rm-A・φm)を与えるステップとを、近似解Umが収束するまで、繰り返し実行するものであり、
前記予測近似値Ψmを求めるステップは、
繰り返し回数m毎の残差切除率を表す変数κmが所定値Kを超えたとき、反復計算を終えるものであることを特徴とする数値計算装置。

【請求項4】
請求項3記載の数値計算装置において、
前記変数κmは、
【数1】
JP0003394989B2_000003t.gifの式で与えられることを特徴とする数値計算装置。

【請求項5】
コンピュータに、物理量Uが満たすべき偏微分方程式を離散化したときの係数行列をA、非斉次項(ソース項)をfとしたとき、
A・U=f
を解かせ、物理量Uの数値計算を行う数値計算用プログラムを記録した記録媒体であって、
物理量Uの初期値U0を設定し、
繰り返し回数mに初期値として0を、摂動量φの初期値として0を与え、残差rの初期値r0として(f-A・U0)を設定し、
繰り返し回数mをインクリメントしながら、
A・φ=rmの予測近似値Ψmを、内部ソルバを有する第1演算部によって、反復計算によって求めるステップと、
残差rmのL2ノルムを最小とする最適化された近似値φmを、第2演算部によって、最適化ルーチンによって、予測近似値Ψmから求め、近似解Um+1として(Um+φm)を与え、残差rm+1として(rm-A・φm)を与えるステップとを、近似解Umが収束するまで、繰り返し実行するものであり、
前記予測近似値Ψmを求めるステップは、
繰り返し回数m毎の残差切除率を表す変数κmが所定値Kを超えたとき、反復計算を終えるものである処理をコンピュータに実行させる数値計算用プログラムを記録した記録媒体。

【請求項6】
請求項5記載の数値計算用プログラムを記録した記録媒体において、
前記変数κmは、
【数1】
JP0003394989B2_000004t.gifの式で与えられることを特徴とする数値計算用プログラムを記録した記録媒体。
発明の詳細な説明 【発明の詳細な説明】

【1】

【発明の属する技術分野】本発明は、圧力、温度などの物理量を対象とした数値計算技術に関するものであり、特に、逐次近似解法アルゴリズムを用いた数値計算技術に属する

【10】
請求項6の発明では、前記請求項5の数値計算用プログラムを記録した記録媒体において、前記数値計算用プログラムで用いる変数κmは、
【数1】
JP0003394989B2_000007t.gifの式で与えられるものとする。

【11】

【発明の実施の形態】以下、本発明の一実施形態について、図面を参照して説明する。

【12】
図1は本発明の一実施形態に係る制御装置の構成を示すブロック図であり、本発明に係る数値計算技術を車両の自動走行制御に応用した例を示している。図1において、車両は、目的値・経路設定部11によって設定される軌道上を許容誤差範囲内で目的地まで自動走行する回路を備えているものとする。走行中は、位置情報検出器12によって、位置・方位情報や速度等の走行情報が検出される。表示器18では設定値や走行状態などの情報がリアルタイムでモニター表示される。

【13】
車両の位置と軌道とのずれφは情報管理装置13によって常に監視されている。そして、このずれφが許容値δmax を越えたとき、そのときの走行情報を出発値として、方向角・角加速度・速度・加速度などの最適化パラメータ制御、すなわち軌道復帰への最適化手段が、本発明に係るアルゴリズムを用いてシミュレーションされる。

【14】
図2は本発明の一実施形態に係る数値計算方法を示すフローチャートである。

【15】
(残差方程式の定式化)
いま求めたい物理量をUとし、その物理量が満たすべき偏微分方程式を離散化したときの係数行列をA、非斉次項(ソース項)をfとすると、解くべき式は、
A・U=f …(1)
で表される。この式は、一般的には、SOR法、ADI法等の逐次近似法によって解くことができる。

【16】
本発明では、以下のような解法を用いる。

【17】
式(1)の解U∞を、近似解Uと摂動量(真の解との差)φとによって、次のように表す。
U∞=U+φ …(2)
本発明に係る解法では、摂動量φを求めて近似解Uを修正していくことによって、式(1)の真の解U∞を求める。

【18】
ここで、近似解Uに対する残差rを次のように定義する。

【19】
r=f-A・U …(3)
式(1)~(3)から、
A・(U+φ)=f
∴ Aφ=f-AU=r
したがって、摂動量φを求めるためには、次の式を解けばよい。
φ=r …(4)

【2】

【発明が解決しようとする課題】従来から、圧力、温度などの物理量を対象とした制御技術において、逐次近似解法アルゴリズムの利用がなされている。逐次近似解法アルゴリズムとは、様々な最適制御理論に基づいて、逐次近似を繰り返し行いつつ、数値計算により最適解を求めるものである。

【20】
(残差切除法のアルゴリズム)
本発明に係るアルゴリズムでは、式(4)の収束解を求めるのではなく、ADI法などによって最も収束勾配の急な最小単位の反復で予測近似値Ψを求め、求めた予測近似値Ψを、最適化制御ルーチンに供給する。そして、最適化制御ルーチンの実行結果が所定の条件を満たすまで、予測近似値Ψの算出を繰り返し実行する。

【21】
<最適化制御ルーチン>
いま、繰り返し回数をmとしたとき、残差のL2ノルムを最小とする合成摂動量φmと新しい近似解Um+1を次のように定義する。
【数2】
JP0003394989B2_000008t.gifm+1=Um+φm …(6)
αl(l=1,2,3,…,L)は残差最小化係数であり、後述する計算方法によって求まる定数である。残差のL2ノルムがより小さくなるように、次のように残差の最小化を行う。

【22】
新しい近似解Um+1が式(6)で表されるとき、これに対する残差rm+1は、式(3)から、次のように表される。

【23】
m+1=f-A・Um+1 …(7)
式(7)に式(6)を代入し、式(3),(5)を用いると、
【数3】
JP0003394989B2_000009t.gif【0024】近似解Um+1に対する残差rm+1のL2ノルム∥rm+1∥は、三次元の場合、(rm+12の内点全ての総和の平方根として、式(9)で与えられる。
【数4】
JP0003394989B2_000010t.gif【0025】この残差rm+1のL2ノルム∥rm+1∥を最小にするように、式(9)の残差最小化係数αll=1,2,3,…,L)を最小二乗法で定める。すなわち、
【数5】
JP0003394989B2_000011t.gif式(10)はαのL元連立方程式となり、これを数値的に解くことによってαlが定まる。αlが定まると、式(5)からφmが求まり、式(6)からUm+1が求まる。残差最小化係数αlの計算方法については後述する。

【26】
mをインクリメントしながら、同様のルーチンを繰り返すことによって、式(9)の残差のL2ノルムを零または最小にする解U∞に収束する。

【27】
図2に示すフローチャートに従って、本実施形態に係る数値計算方法について説明する。

【28】
まず、ステップS1において、対象物理量Uの初期値U0を設定する。そして、ステップS2において、繰り返し回数mに初期値として0を、摂動量φの初期値として0を与え、残差rの初期値r0として(f-A0)を設定する。

【29】
以下、近似解Umが収束するまで、繰り返し回数mをインクリメントしながら、ステップS3~S9を繰り返し実行する。

【3】
ところが、逐次近似解法アルゴリズムを用いた場合には、制御対象の各パラメータの解を求める際に、計算の反復回数が多くなるために、計算時間の点から、高速かつ高精度の制御が困難になっている。このため、様々な分野の制御技術に応用できるようにするために、逐次近似解法アルゴリズムの収束性の向上が、望まれている。

【30】
ステップS3~S8では、ADI法、SOR法、CG法などの逐次近似法を実行する内部ソルバによって、
A・φ=rm
の近似解(予測近似値)Ψm を求める。ところがこの式は、既存の逐次近似法で解こうとすると、実用問題に対しては収束までに多くの繰り返しを必要とし、必然的に計算機使用時間も膨大となる。

【31】
そこで本実施形態では、逐次計算の反復回数の上限Nを設定する(ステップS7,S8)とともに、残差の減少の割合すなわち残差切除率を表す変数κmを導入し(ステップS5)、この変数κmが所定値Kを超えたとき、反復計算を終えるものとしている(ステップS6)。

【32】
逐次計算の反復回数を固定した場合には、収束判定までのループの繰り返しm毎に、残差の減少の割合(残差切除率)が大きく変動する。この変動は、各繰り返しmにおける予測近似値Ψm の近似度合に起因するものである。このために、解が収束するまでのループの繰り返し回数が必要以上に増加してしまうおそれがある。

【33】
これに対して、本実施形態では、予測近似値Ψm を求めるステップにおいて、残差切除率を表す変数κm が所定値Kを超えたときに反復計算を終えるものとしているので、ループの繰り返しm毎の残差切除率を一定値以上に保つことができ、したがって、解が収束するまでに要するループの繰り返し回数を減少させることができる。

【34】
ここでは、残差切除率を表す変数κmを次のように定義する。
【数6】
JP0003394989B2_000012t.gif【0035】この変数κmは、反復計算を繰り返すにつれて、予測近似値Ψmの精度が高くなるために、1.0に限りなく近づく。したがって、この変数κmの値が任意に設定した所定値Kを超えるまで反復計算を行うことによって、ループの繰り返し回数m毎の残差切除率はほぼ一定になる。

【36】
ステップS3~S8において予測近似値Ψ0 が求められると、ステップS9において、この予測近似値Ψ0 を用いて、最適化ルーチンによって、次回の残差r1が最小となるように、最適化された近似値φ0を求める。この最適化された近似値φ0によって、第1次の物理量の近似値U1および第1次の残差r1は、式(2),(3)から、次式で求められる。
1=U0+φ0 …(12)
1=f-A1 …(13)
式(13)に式(12)を代入して、
1=f-A(U0+φ0
=f-A0-Aφ0
=r0-Aφ0 …(14)
よって、式(14)により、第1次の残差r1が求められる。

【37】
ステップS11でmをインクリメントしてステップS3に戻り、同様の計算を行う。このような処理を繰り返すことによって、
(U2 ,r2 ),(U3 ,r3 )…(Um ,rm )…
が順次求められ、残差ノルム∥rm∥が減少していく。

【38】
繰り返し回数がmのとき、ステップS4において予測近似値Ψmを求める方程式は、
φ=rm
となり、ステップS9では、最適化制御ルーチンによって最適化された近似値φmによって、
m+1=Um+φm
m+1=rm-Aφm
と求められる。

【39】
なお、∥rm∥=0となったときが解であるが、境界条件によっては、最終残差rfinal が残る場合がある。すなわち、境界条件が全て物理量の微分で与えられる場合(ノイマン問題)、数学的にソース項fと境界条件との間で満たすべき関係(ノイマンの拘束条件)があり、これが満たされないために、残差ノルムが
∥rm∥=∥rm-1∥=∥rm-2∥…
となってしまう。

【4】
前記の問題に鑑み、本発明は、従来よりも収束性の高い逐次近似解法アルゴリズムを用いた数値計方法、数値計算装置および数値計算用プログラムを記録した記録媒体を提供することを課題とする。

【40】
他の解法では、このようなノイマン問題の場合、そのままでは収束しないため、最終残差rfinal によって、ソース項fあるいは境界条件を修正して収束解を得ようとする。これは、問題を変更していることになる。ところが、本実施形態に係る解法では、このような修正を行わないでも自動的に最終残差rfinal が確定し、収束解が得られる特性を持っている。

【41】
(残差最小化係数αlの求め方)
いま、表記を簡単にするために、次のようにδを導入する。
δm =Ψm
δm-l+1=φm-l+1 (l=2,3,…,L)
とおくと、最適化された近似値φmは、
【数7】
JP0003394989B2_000013t.gifと表される。よって、
【数8】
JP0003394989B2_000014t.gifΣ記号を外して表記すると、
m+1=rm-(α1A・δm+α2A・δm-1+α3A・δm-2+…+αLA・δm-L+1
となる。

【42】
残差最小化係数αlは、例えば最小自乗法を用いて、(m+1)次の残差ノルム∥rm+1∥(自乗和の平方根)を最小にするという条件を与えることによって、求められる。
【数9】
JP0003394989B2_000015t.gifこれをαlで微分すると、
【数10】
JP0003394989B2_000016t.gif【0043】この式をαlについて解くことによって、残差最小化係数αlを求めることができる。微分を実行すれば、
【数11】
JP0003394989B2_000017t.gifとなる。ここで、l'=1,2,…,Lである。この連立一次方程式を解くことによって、残差最小化係数αlを求めることができる。

【44】
例えば、L=3のときは、上式は、
【数12】
JP0003394989B2_000018t.gifのような3元連立一次方程式になる。そして、
φm =α1δm+α2δm-1+α3δm-2
m+1=rm-(α1δm+α2δm-1+α3δm-2
=rm-Aφm
m+1=Um+φm
となる。

【45】
図1に示す車両の走行制御においては、まず予測修正量演算部14によって、現在の走行情報を基にして予測修正量Ψm が反復計算される。ここでは、繰り返し回数毎の残差切除率κm が所定値Kを越えるまで反復計算される。次に、最適加算修正量演算部15によってm次の加算修正量φm が演算される。予測修正量演算部14および最適加算修正量演算部15によって求められた解は、最終的に修正量収束判定部16によって収束が判定され、所定の軌道に復帰するための最適制御方法が決定される。決定された制御情報は方向制御部17に送られ、車両の走行経路が修正される。

【46】
なお、ここでは車両の走行制御を例にとって本発明について説明したが、他の分野の様々な制御にも本発明は利用可能であることはいうまでもない。

【47】
また、上記の数値計算方法は、当該方法を実現するためのプログラムを実行するコンピュータを備えた装置によって容易に実現することができる。また、当該方法を実現するためのプログラムをコンピュータ読み取り可能な記録媒体に記録して、この記録媒体に記録したプログラムをコンピュータに実行させることによって実現することができる。

【48】
ここで、従来技術との比較について補足説明する。従来、機器の制御に用いられてきた技術には、比例制御やファジィ制御等がある。位置制御を例にとると、比例制御の場合には、基本的には現在位置と設定位置との比較を2値化して判断し、そのずれを限りなく零に近づけようとするため、微調整動作が頻繁に行われる。このため、動作が安定しない、また外乱に対して柔軟に対処できないという欠点があった。

【49】
このような欠点を解消する技術として、ファジィ理論を採用した制御技術が開発された。これは、その状態量を「大きくずれている」「少しずれている」「微妙にずれている」などと「あいまいに」判断し、これらの状態に対して「大きく修正する」「少し修正する」「そのまま」のような「あいまいな」ルールに基づいて制御を行うものである。この技術によって、AI制御に近い「滑らかな制御」が実現されたが、情報量と制御パラメータに比例してルール数が増加するため、その設定とチューニング作業に手間がかかる、または主観的な設定が可能なために制御結果にばらつきが生じるという欠点もあった。

【5】

【課題を解決するための手段】前記の課題を解決するために、請求項1の発明が講じた解決手段は、物理量Uが満たすべき偏微分方程式を離散化したときの係数行列をA、非斉次項(ソース項)をfとしたとき、コンピュータを用いて、A・U=fを解き、物理量Uの数値計算を行う方法として、物理量Uの初期値U0を設定し、繰り返し回数mに初期値として0を、摂動量φの初期値として0を与え、残差rの初期値r0として(f-A・U0)を設定し、繰り返し回数mをインクリメントしながら、A・φ=rmの予測近似値Ψmを内部ソルバを有する第1の演算部によって反復計算によって求めるステップと、残差rmのL2ノルムを最小とする最適化された近似値φmを第2演算部によって最適化ルーチンによって予測近似値Ψmから求め、近似解Um+1として(Um+φm)を与え、残差rm+1として(rm-A・φm)を与えるステップとを近似解Umが収束するまで繰り返し実行するものであり、前記予測近似値Ψmを求めるステップは、繰り返し回数m毎の残差切除率を表す変数κmが所定値Kを超えたとき、反復計算を終えるものである。

【50】
本発明は、制御対象を目的の状態に制御するための最適手段を決定するために用いられる技術であり、その基となるのは制御パラメータに対応した非線形方程式である。このため、対象によっては目的の状態が数値的に明確でないケースへも適応可能であり、また最適な状態を自動的にチューニングする機能も備えている。

【51】
(他の応用例その1)
図3は本発明に係るアルゴリズムの他の応用例を示す図である。図3では、数値計算への応用例として、2次元または3次元楕円型方程式法の格子生成手法への適用について示すフローである。

【52】
図3で用いた符号の意味は下記のとおりである。
U :物理面における格子点の座標値(x,y,z)
P :格子生成の制御関数
A :線型化された係数行列
Ψ(m) :m次ステップでの予測修正値
φ(m) :m次ステップでの最適化された修正値
αl :未定係数行列
(m) :m次ステップでの残差
L :未定係数行列の個数

【53】
ステップS21では、m次ステップにおける値に基づいて、離散化された楕円型方程式系の係数行列を線型化させる。ステップS22では、κ技術によって繰り返し回数を制御して、(m+1)次ステップにおける予測修正量Ψ(m) を求める。ステップS23では、(m+1)次ステップの残差を最小(ゼロ)にするための未定係数行列αiを求める。ステップS24では、予測修正値Ψ(m) 、(m-1)次以前の最適化された修正(φ(m-1) ,…,φ(m-L+1) )および未定係数行列に基づいて(m+1)次ステップでの最適化された修正値φ(m+1) を求める。ステップS25では、最適化された修正値φ(m+1) に基づいて、物理面における格子点の座標値を修正する。

【54】
図4は本応用例によって生成された格子、図5は比較例としてのSOR法によって生成された格子である。図4および図5では、気流中におかれた楕円柱周りの2次元格子生成を行っている。格子点数は90×90である。本応用例では、、残差切除法の内部ソルバとしてADI法を用いた。

【55】
図6は収束速度を比較するためのグラフである。図6において、a1,a2は本応用例におけるX方向およびY方向の残差、b1,b2は比較例におけるX方向およびY方向の残差である。図6から分かるように、本発明に係るアルゴリズムを利用した場合には、SOR法と比較して6倍以上高速に収束している。

【56】
(他の応用例その2)
図7は本発明に係るアルゴリズムの他の応用例を示す図である。図7では、数値計算への応用例として、流体解析のMAC法の圧力計算部に残差切除法を適用した例を示すフローである。

【57】
図7で用いた符号の意味は下記の通りである。
t :時刻
U :流速
P :圧力
Ψm :m次ステップでの予測修正量
φm :m次ステップでの最適化された修正量
α :最適化係数
m :m次ステップでの残差

【58】
S31は初期残差r0 の算出を含む初期化部分、S32は予測近似値Ψm の算出部、S33はκm の算出部、S34はS33で算出されたκm によって繰り返し数を制御されたS32およびS33を含むループ、S35は最適化係数α、最適化された修正量φm および次ステップでの残差rm+1 の算出部である。

【59】
本応用例による効果を評価するために、MAC法によって、一様流中に置かれた円柱周りの二次元解析を行った。比較例としてはSOR法を用いた。格子点数は60×60、主流速度は5cm/sである。本応用例では、残差切除法の内部ソルバとしてSOR法を用いた。

【6】
そして、請求項2の発明では、前記請求項の数値計算方法における変数κmは、
【数1】
JP0003394989B2_000005t.gifの式で与えられるものとする。

【60】
図8はタイムステップ数とCPU時間との関係を示すグラフである。図8において、aは本応用例による結果、bは比較例としてのSOR法による結果を示す。また図9はステップ数と相対残差/相対誤差との関係を示すグラフである。図9において、a1,a2は本応用例による結果、b1,b2は比較例としてのSOR法による結果を示す。図8および図9から、本応用例によると、SOR法に比べて5倍に高速化されていることが分かる。

【61】

【発明の効果】以上のように本発明によると、予測近似値Ψmを求めるステップにおいて、繰り返し回数m毎の残差切除率を表す変数κmが所定値Kを超えたとき、反復計算を終えるので、ループの繰り返しm毎の残差切除率を一定値以上に保つことができ、これにより、収束に必要なmの値を減少させることができる。したがって、逐次近似解法アルゴリズムの収束性が向上するので、えば制御に用いた場合に、高速かつ高精度の制御が容易になる。

【7】
また、請求項3の発明が講じた解決手段は、物理量Uが満たすべき偏微分方程式を離散化したときの係数行列をA、非斉次項(ソース項)をfとしたとき、A・U=fを解き、物理量Uの数値計算を行う装置として、物理量Uの初期値U0を設定し、繰り返し回数mに初期値として0を、摂動量φの初期値として0を与え、残差rの初期値r0として(f-A・U0)を設定し、繰り返し回数mをインクリメントしながら、A・φ=rmの予測近似値Ψmを内部ソルバを有する第1演算部によって反復計算によって求めるステップと、残差rmのL2ノルムを最小とする最適化された近似値φmを第2演算部によって最適化ルーチンによって予測近似値Ψmから求め、近似解Um+1として(Um+φm)を与え、残差rm+1として(rm-A・φm)を与えるステップとを近似解Umが収束するまで、繰り返し実行するものであり、前記予測近似値Ψmを求めるステップは、繰り返し回数m毎の残差切除率を表す変数κmが所定値Kを超えたとき反復計算を終えるものである。

【8】
そして、請求項4の発明では、前記請求項の数値計算装置における変数κmは、
【数1】
JP0003394989B2_000006t.gifの式で与えられるものとする。

【9】
また、請求項5の発明が講じた解決手段は、コンピュータに、物理量Uが満たすべき偏微分方程式を離散化したときの係数行列をA、非斉次項(ソース項)をfとしたとき、A・U=fを解かせ、物理量Uの数値計算を行う数値計算用プログラムを記録した記録媒体として、物理量Uの初期値U0を設定し、繰り返し回数mに初期値として0を、摂動量φの初期値として0を与え、残差rの初期値r0として(f-A・U0)を設定し、繰り返し回数mをインクリメントしながら、A・φ=rmの予測近似値Ψmを内部ソルバを有する第1演算部によって反復計算によって求めるステップと、残差rmのL2ノルムを最小とする最適化された近似値φmを第2演算部によって最適化ルーチンによって予測近似値Ψmから求め、近似解Um+1として(Um+φm)を与え、残差rm+1として(rm-A・φm)を与えるステップとを近似解Umが収束するまで繰り返し実行するものであり、前記予測近似値Ψmを求めるステップは、繰り返し回数m毎の残差切除率を表す変数κmが所定値Kを超えたとき反復計算を終えるものである処理をコンピュータに実行させる数値計算用プログラムを記録したものである。
図面
【図1】
0
【図8】
1
【図4】
2
【図5】
3
【図2】
4
【図3】
5
【図6】
6
【図9】
7
【図7】
8