(57)【要約】【課題】逐次近似解法アルゴリズムを用いた数値計算において、計算精度を保ちつつ、計算に必要な記憶領域を抑えるとともに、計算時間を短縮する。【解決手段】ベクトル列Ａ・φｍ　の要素を、サンプリングしてメモリに保存する。このとき、空間的なサンプリング（△）と、物理現象を踏まえた局所的なサンプリング（◎○）とを組み合わせる（例Ｃ）。そして、修正近似値φｍ　を求めるために用いる残差最小化係数αｌｍ（ｌ＝１，…，Ｌ）を、メモリに保存したＡ・φｋ　（ｋ＝ｍ−Ｌ＋１，…，ｍ−１）の要素を用いて、近似的に求める。【選択図】　図２
2004118618.tif 000002
【特許請求の範囲】
【請求項１】
物理量Ｕが満たすべき偏微分方程式を離散化したときの係数行列（Ｎ行Ｎ列、Ｎは正の整数）をＡ、非斉次項（ソース項）をｆとしたとき、コンピュータを用いて、
Ａ・Ｕ＝ｆ
を解き、物理量Ｕの数値計算を行う方法であって、
物理量Ｕの初期値Ｕ０を設定し、
繰り返し回数ｍに初期値として０を、摂動量φの初期値として０を与え、残差ｒの初期値ｒ０として（ｆ−Ａ・Ｕ０）を設定し、
繰り返し回数ｍをインクリメントしながら、
Ａ・φ＝ｒｍの予測近似値Ψｍを、内部ソルバを有する第１演算部によって、反復計算によって求める第１のステップと、
残差ｒｍのＬ２ノルムを最小とする修正近似値φｍを、第２演算部によって、最適化ルーチンによって、予測近似値Ψｍから求め、近似解Ｕｍ＋１として（Ｕｍ＋φｍ）を与え、残差ｒｍ＋１として（ｒｍ−Ａ・φｍ）を与える第２のステップとを、近似解Ｕｍが収束するまで、繰り返し実行するものであり、
前記第２のステップは、
求めたベクトル列Ａ・φｍ　の要素を、所定のサンプリング方法によってサンプリングして、メモリに保存するものであり、かつ、
修正近似値φｍ　を求めるために用いる残差最小化係数αｌｍ（ｌ＝１，…，Ｌ）を、前記メモリに保存されているベクトル列Ａ・φｋ　（ｋ＝ｍ−Ｌ＋１，…，ｍ−１）の要素を用いて、近似的に求める
ことを特徴とする数値計算方法。
【請求項２】
請求項１記載の数値計算方法において、
前記第２のステップは、
ベクトル列Ａ・φｍ　の要素ｂ１　，ｂ２　，…，ｂＮ　のサンプリングにおいて、要素ｂｉ（ｉ∈Ω）を選択する（ただし、
【数１】
2004118618.tif 000018
ｌｇは整数、βは実数、ｆｉ　はソース項ｆの要素、ａｉｉは行列Ａのｉ行ｉ列の対角項）
ことを特徴とする数値計算方法。
【請求項３】
物理量Ｕが満たすべき偏微分方程式を離散化したときの係数行列（Ｎ行Ｎ列、Ｎは正の整数）をＡ、非斉次項（ソース項）をｆとしたとき、
Ａ・Ｕ＝ｆ
を解き、物理量Ｕの数値計算を行う装置であって、
物理量Ｕの初期値Ｕ０を設定し、
繰り返し回数ｍに初期値として０を、摂動量φの初期値として０を与え、残差ｒの初期値ｒ０として（ｆ−Ａ・Ｕ０）を設定し、
繰り返し回数ｍをインクリメントしながら、
Ａ・φ＝ｒｍの予測近似値Ψｍを、内部ソルバを有する第１演算部によって、反復計算によって求める第１のステップと、
残差ｒｍのＬ２ノルムを最小とする修正近似値φｍを、第２演算部によって、最適化ルーチンによって、予測近似値Ψｍから求め、近似解Ｕｍ＋１として（Ｕｍ＋φｍ）を与え、残差ｒｍ＋１として（ｒｍ−Ａ・φｍ）を与える第２のステップとを、近似解Ｕｍが収束するまで、繰り返し実行するものであり、
前記第２のステップは、
求めたベクトル列Ａ・φｍ　の要素を、所定のサンプリング方法によってサンプリングして、メモリに保存するものであり、かつ、
修正近似値φｍ　を求めるために用いる残差最小化係数αｌｍ（ｌ＝１，…，Ｌ）を、前記メモリに保存されているベクトル列Ａ・φｋ　（ｋ＝ｍ−Ｌ＋１，…，ｍ−１）の要素を用いて、近似的に求める
ことを特徴とする数値計算装置。
【請求項４】
請求項３記載の数値計算装置において、
前記第２のステップは、
ベクトル列Ａ・φｍ　の要素ｂ１　，ｂ２　，…，ｂＮ　のサンプリングにおいて、要素ｂｉ（ｉ∈Ω）を選択する（ただし、
【数１】
2004118618.tif 000019
ｌｇは整数、βは実数、ｆｉ　はソース項ｆの要素、ａｉｉは行列Ａのｉ行ｉ列の対角項）
ことを特徴とする数値計算装置。
【請求項５】
コンピュータに、物理量Ｕが満たすべき偏微分方程式を離散化したときの係数行列（Ｎ行Ｎ列、Ｎは正の整数）をＡ、非斉次項（ソース項）をｆとしたとき、
Ａ・Ｕ＝ｆ
を解かせ、物理量Ｕの数値計算を行う数値計算用プログラムであって、
物理量Ｕの初期値Ｕ０を設定し、
繰り返し回数ｍに初期値として０を、摂動量φの初期値として０を与え、残差ｒの初期値ｒ０として（ｆ−Ａ・Ｕ０）を設定し、
繰り返し回数ｍをインクリメントしながら、
Ａ・φ＝ｒｍの予測近似値Ψｍを、内部ソルバを有する第１演算部によって、反復計算によって求める第１のステップと、
残差ｒｍのＬ２ノルムを最小とする修正近似値φｍを、第２演算部によって、最適化ルーチンによって、予測近似値Ψｍから求め、近似解Ｕｍ＋１として（Ｕｍ＋φｍ）を与え、残差ｒｍ＋１として（ｒｍ−Ａ・φｍ）を与える第２のステップとを、近似解Ｕｍが収束するまで、繰り返し実行するものであり、
前記第２のステップは、
求めたベクトル列Ａ・φｍ　の要素を、所定のサンプリング方法によってサンプリングして、メモリに保存するものであり、かつ、
修正近似値φｍ　を求めるために用いる残差最小化係数αｌｍ（ｌ＝１，…，Ｌ）を、前記メモリに保存されているベクトル列Ａ・φｋ　（ｋ＝ｍ−Ｌ＋１，…，ｍ−１）の要素を用いて、近似的に求める
処理をコンピュータに実行させる数値計算用プログラム。
【請求項６】
請求項５記載の数値計算用プログラムにおいて、
前記第２のステップは、
ベクトル列Ａ・φｍ　の要素ｂ１　，ｂ２　，…，ｂＮ　のサンプリングにおいて、要素ｂｉ（ｉ∈Ω）を選択する（ただし、
【数１】
2004118618.tif 000020
ｌｇは整数、βは実数、ｆｉ　はソース項ｆの要素、ａｉｉは行列Ａのｉ行ｉ列の対角項）
ことを特徴とする数値計算用プログラム。
【発明の詳細な説明】
【０００１】
【発明の属する技術分野】
本発明は、コンピュータを用いた数値計算技術に関するものである。
【０００２】
【従来の技術】
従来から、圧力、温度などの物理量を対象とした制御技術において、逐次近似解法アルゴリズムの利用がなされている。そして、制御対象の各パラメータの解を、高速かつ高精度に求めることができる収束性の高い逐次近似解法アルゴリズムとして、残差切除法が知られている。さらに、特開２００１−１３４３０４号公報では、残差切除法において収束性の更なる向上を実現するような技術が、すでに提案されれている。
【０００３】
【特許文献１】
特開２００１−１３４３０４号公報
【０００４】
【発明が解決しようとする課題】
残差切除法のような逐次近似解法アルゴリズムでは、近似解に修正量を加えていくことによって解の精度を高めているが、修正量の計算には過去の修正量の列が用いられ、この列が長くなればなるほど修正量の精度が向上する。ところがその一方で、過去の修正量の列を長くとればとるほど、計算に必要な記憶領域は増大し、かつ、修正量を求めるための演算時間も長くなる、という問題が生じる。
【０００５】
前記の問題に鑑み、本発明は、逐次近似解法アルゴリズムを用いた数値計算において、計算精度を保ちつつ、計算に必要な記憶領域を抑えるとともに、計算時間を短縮することを課題とする。
【０００６】
【課題を解決するための手段】
前記の問題を解決するために、請求項１の発明が講じた解決手段は、物理量Ｕが満たすべき偏微分方程式を離散化したときの係数行列（Ｎ行Ｎ列、Ｎは正の整数）をＡ、非斉次項（ソース項）をｆとしたとき、コンピュータを用いて、Ａ・Ｕ＝ｆを解き、物理量Ｕの数値計算を行う方法として、物理量Ｕの初期値Ｕ０を設定し、繰り返し回数ｍに初期値として０を、摂動量φの初期値として０を与え、残差ｒの初期値ｒ０として（ｆ−Ａ・Ｕ０）を設定し、繰り返し回数ｍをインクリメントしながら、Ａ・φ＝ｒｍの予測近似値Ψｍを内部ソルバを有する第１演算部によって反復計算によって求める第１のステップと、残差ｒｍのＬ２ノルムを最小とする修正近似値φｍを、第２演算部によって最適化ルーチンによって、予測近似値Ψｍから求め、近似解Ｕｍ＋１として（Ｕｍ＋φｍ）を与え、残差ｒｍ＋１として（ｒｍ−Ａ・φｍ）を与える第２のステップとを、近似解Ｕｍが収束するまで繰り返し実行するものであり、前記第２のステップは、求めたベクトル列Ａ・φｍ　の要素を、所定のサンプリング方法によってサンプリングして、メモリに保存するものであり、かつ、修正近似値φｍ　を求めるために用いる残差最小化係数αｌｍ（ｌ＝１，…，Ｌ）を、前記メモリに保存されているベクトル列Ａ・φｋ　（ｋ＝ｍ−Ｌ＋１，…，ｍ−１）の要素を用いて近似的に求めるものである。
【０００７】
そして、請求項２の発明では、前記請求項１の数値計算方法における第２のステップは、ベクトル列Ａ・φｍ　の要素ｂ１　，ｂ２　，…，ｂＮ　のサンプリングにおいて、要素ｂｉ（ｉ∈Ω）を選択する（ただし、
【数１】
2004118618.tif 000003
ｌｇは整数、βは実数、ｆｉ　はソース項ｆの要素、ａｉｉは行列Ａのｉ行ｉ列の対角項）ものとする。
【０００８】
また、請求項３の発明が講じた解決手段は、物理量Ｕが満たすべき偏微分方程式を離散化したときの係数行列（Ｎ行Ｎ列、Ｎは正の整数）をＡ、非斉次項（ソース項）をｆとしたとき、Ａ・Ｕ＝ｆを解き、物理量Ｕの数値計算を行う装置として、物理量Ｕの初期値Ｕ０を設定し、繰り返し回数ｍに初期値として０を、摂動量φの初期値として０を与え、残差ｒの初期値ｒ０として（ｆ−Ａ・Ｕ０）を設定し、繰り返し回数ｍをインクリメントしながら、Ａ・φ＝ｒｍの予測近似値Ψｍを内部ソルバを有する第１演算部によって反復計算によって求める第１のステップと、残差ｒｍのＬ２ノルムを最小とする修正近似値φｍを第２演算部によって最適化ルーチンによって予測近似値Ψｍから求め、近似解Ｕｍ＋１として（Ｕｍ＋φｍ）を与え、残差ｒｍ＋１として（ｒｍ−Ａ・φｍ）を与える第２のステップとを、近似解Ｕｍが収束するまで繰り返し実行するものであり、前記第２のステップは、求めたベクトル列Ａ・φｍ　の要素を、所定のサンプリング方法によってサンプリングしてメモリに保存するものであり、かつ、修正近似値φｍ　を求めるために用いる残差最小化係数αｌｍ（ｌ＝１，…，Ｌ）を、前記メモリに保存されているベクトル列Ａ・φｋ　（ｋ＝ｍ−Ｌ＋１，…，ｍ−１）の要素を用いて近似的に求めるものである。
【０００９】
請求項４の発明では、前記請求項３の数値計算装置における第２のステップは、ベクトル列Ａ・φｍ　の要素ｂ１　，ｂ２　，…，ｂＮ　のサンプリングにおいて、要素ｂｉ（ｉ∈Ω）を選択する（ただし、
【数１】
2004118618.tif 000004
ｌｇは整数、βは実数、ｆｉ　はソース項ｆの要素、ａｉｉは行列Ａのｉ行ｉ列の対角項）ものとする。
【００１０】
また、請求項５の発明が講じた解決手段は、コンピュータに、物理量Ｕが満たすべき偏微分方程式を離散化したときの係数行列（Ｎ行Ｎ列、Ｎは正の整数）をＡ、非斉次項（ソース項）をｆとしたとき、Ａ・Ｕ＝ｆを解かせ、物理量Ｕの数値計算を行う数値計算用プログラムとして、物理量Ｕの初期値Ｕ０を設定し、繰り返し回数ｍに初期値として０を、摂動量φの初期値として０を与え、残差ｒの初期値ｒ０として（ｆ−Ａ・Ｕ０）を設定し、繰り返し回数ｍをインクリメントしながら、Ａ・φ＝ｒｍの予測近似値Ψｍを内部ソルバを有する第１演算部によって反復計算によって求める第１のステップと、残差ｒｍのＬ２ノルムを最小とする修正近似値φｍを、第２演算部によって最適化ルーチンによって予測近似値Ψｍから求め、近似解Ｕｍ＋１として（Ｕｍ＋φｍ）を与え、残差ｒｍ＋１として（ｒｍ−Ａ・φｍ）を与える第２のステップとを、近似解Ｕｍが収束するまで繰り返し実行するものであり、前記第２のステップは、求めたベクトル列Ａ・φｍ　の要素を、所定のサンプリング方法によってサンプリングしてメモリに保存するものであり、かつ、修正近似値φｍ　を求めるために用いる残差最小化係数αｌｍ（ｌ＝１，…，Ｌ）を、前記メモリに保存されているベクトル列Ａ・φｋ　（ｋ＝ｍ−Ｌ＋１，…，ｍ−１）の要素を用いて近似的に求める処理をコンピュータに実行させるものである。
【００１１】
請求項６の発明では、前記請求項５の数値計算用プログラムにおける第２のステップは、ベクトル列Ａ・φｍ　の要素ｂ１　，ｂ２　，…，ｂＮ　のサンプリングにおいて、要素ｂｉ（ｉ∈Ω）を選択する（ただし、
【数１】
2004118618.tif 000005
ｌｇは整数、βは実数、ｆｉ　はソース項ｆの要素、ａｉｉは行列Ａのｉ行ｉ列の対角項）ものとする。
【００１２】
【発明の実施の形態】
以下、本発明の一実施形態について、図面を参照して説明する。
【００１３】
図１は本発明の一実施形態に係る数値計算方法を示すフローチャートである。
【００１４】
（残差方程式の定式化）
いま求めたい物理量をＵとし、その物理量Ｕが満たすべき偏微分方程式を離散化したときの係数行列（Ｎ行Ｎ列、Ｎは正の整数）をＡ、非斉次項（ソース項）をｆとすると、解くべき式は、
Ａ・Ｕ＝ｆ　…（１）
で表される。この式は、一般的には、ＳＯＲ法、ＡＤＩ法等の逐次近似法によって解くことができる。
【００１５】
本発明では、以下のような解法を用いる。
【００１６】
式（１）の解Ｕ∞を、近似解Ｕと摂動量（真の解との差）φとによって、次のように表す。
Ｕ∞＝Ｕ＋φ　…（２）
本発明に係る解法では、摂動量φを求めて近似解Ｕを修正していくことによって、式（１）の真の解Ｕ∞を求める。
【００１７】
ここで、近似解Ｕに対する残差ｒを次のように定義する。
ｒ＝ｆ−Ａ・Ｕ　…（３）
式（１）〜（３）から、
Ａ・（Ｕ＋φ）＝ｆ
∴　Ａφ＝ｆ−ＡＵ
＝ｒ
【００１８】
したがって、摂動量φを求めるためには、次の式を解けばよい。
Ａφ＝ｒ　…（４）
【００１９】
（残差切除法のアルゴリズム）
本発明に係るアルゴリズムでは、式（４）の収束解を求めるのではなく、ＡＤＩ法などによって最も収束勾配の急な最小単位の反復で予測近似値Ψを求め、求めた予測近似値Ψを、最適化制御ルーチンに供給する。そして、最適化制御ルーチンの実行結果が所定の条件を満たすまで、予測近似値Ψの算出を繰り返し実行する。
【００２０】
＜最適化制御ルーチン＞
いま、繰り返し回数をｍとしたとき、残差のＬ２ノルムを最小とする合成摂動量φｍと新しい近似解Ｕｍ＋１を次のように定義する。
【数２】
2004118618.tif 000006
Ｕｍ＋１＝Ｕｍ＋φｍ　…（６）
αｌｍ（ｌ＝１，２，３，…，Ｌ）は残差最小化係数であり、後述する計算方法によって求まる定数である。残差のＬ２ノルムがより小さくなるように、次のように残差の最小化を行う。
【００２１】
新しい近似解Ｕｍ＋１が式（６）で表されるとき、これに対する残差ｒｍ＋１は、式（３）から、次のように表される。
ｒｍ＋１＝ｆ−Ａ・Ｕｍ＋１　…（７）
式（７）に式（６）を代入し、式（３），（５）を用いると、
【数３】
2004118618.tif 000007
【００２２】
近似解Ｕｍ＋１に対する残差ｒｍ＋１のＬ２ノルム‖ｒｍ＋１‖は、（ｒｍ＋１）２の全点の総和の平方根として、式（９）で与えられる。
【数４】
2004118618.tif 000008
【００２３】
この残差ｒｍ＋１のＬ２ノルム‖ｒｍ＋１‖を最小にするように、式（９）の残差最小化係数αｌｍ（ｌ＝１〜Ｌ）を最小二乗法で定める。すなわち、
【数５】
2004118618.tif 000009
式（１０）はαｌｍのＬ元連立方程式となり、これを数値的に解くことによってαｌｍが定まる。αｌｍが定まると、式（５）からφｍが求まり、式（６）からＵｍ＋１が求まる。残差最小化係数αｌｍの計算方法については後述する。
【００２４】
ｍをインクリメントしながら、同様のルーチンを繰り返すことによって、式（９）の残差のＬ２ノルムを零または最小にする解Ｕ∞に収束する。
【００２５】
図１に示すフローチャートに従って、本実施形態に係る数値計算方法について説明する。この数値計算方法は、当該方法を実現するためのプログラムを、コンピュータに実行させることによって、実施することができる。
【００２６】
まず、ステップＳ１において、対象物理量Ｕの初期値Ｕ０を設定する。そして、ステップＳ２において、繰り返し回数ｍに初期値として０を、摂動量φの初期値として０を与え、残差ｒの初期値ｒ０として（ｆ−ＡＵ０）を設定する。
【００２７】
以下、近似解Ｕｍが収束するまで、繰り返し回数ｍをインクリメントしながら、ステップＳ３〜Ｓ７を繰り返し実行する。
【００２８】
ステップＳ３〜Ｓ６では、ＡＤＩ法、ＳＯＲ法、ＣＧ法などの逐次近似法を実行する内部ソルバを有する第１演算部によって、
Ａ・φ＝ｒｍ
の近似解（予測近似値）Ψｍ　を求める。ところがこの式は、既存の逐次近似法で解こうとすると、実用問題に対しては収束までに多くの繰り返しを必要とし、必然的に計算機使用時間も膨大となる。そこで本実施形態では、逐次計算の反復回数の上限Ｎを設定する（ステップＳ５）。
【００２９】
ステップＳ３〜Ｓ６において予測近似値Ψ０　が求められると、ステップＳ７において、この予測近似値Ψ０　を用いて、第２演算部による最適化ルーチンによって、次回の残差ｒ１が最小となるように、修正近似値φ０を求める。この修正近似値φ０によって、第１次の物理量の近似値Ｕ１および第１次の残差ｒ１は、式（２），（３）から、次式で求められる。
Ｕ１＝Ｕ０＋φ０　　…（１２）
ｒ１＝ｆ−ＡＵ１　　　　…（１３）
式（１３）に式（１２）を代入して、
ｒ１＝ｆ−Ａ（Ｕ０＋φ０）
＝ｆ−ＡＵ０−Ａφ０
＝ｒ０−Ａφ０　　…（１４）
よって、式（１４）により、第１次の残差ｒ１が求められる。
【００３０】
ステップＳ９でｍをインクリメントしてステップＳ３に戻り、同様の計算を行う。このような処理を繰り返すことによって、
（Ｕ２　，ｒ２　），（Ｕ３　，ｒ３　）…（Ｕｍ　，ｒｍ　）…
が順次求められ、残差ノルム‖ｒｍ‖が減少していく。
【００３１】
繰り返し回数がｍのとき、ステップＳ４において予測近似値Ψｍを求める方程式は、
Ａφ＝ｒｍ
となり、ステップＳ７では、最適化制御ルーチンによって最適化された修正近似値φｍによって、
Ｕｍ＋１＝Ｕｍ＋φｍ
ｒｍ＋１＝ｒｍ−Ａφｍ
と求められる。
【００３２】
（残差最小化係数αｌの求め方）
いま、表記を簡単にするために、次のようにδを導入する。
δｍ　　＝Ψｍ
δｍ−ｌ＋１＝φｍ−ｌ＋１　（ｌ＝２，３，…，Ｌ）
とおくと、修正近似値φｍは、
【数６】
2004118618.tif 000010
と表される。よって、
【数７】
2004118618.tif 000011
Σ記号を外して表記すると、
ｒｍ＋１＝ｒｍ−（α１ｍＡδｍ＋α２ｍＡδｍ−１＋α３ｍＡδｍ−２＋…＋αＬｍＡδｍ−Ｌ＋１）となる。
【００３３】
残差最小化係数αｌｍは、例えば最小自乗法を用いて、（ｍ＋１）次の残差ノルム‖ｒｍ＋１‖（自乗和の平方根）を最小にするという条件を与えることによって、求められる。
【数８】
2004118618.tif 000012
これをαｌｍで微分すると、
【数９】
2004118618.tif 000013
この式をαｌｍについて解くことによって、残差最小化係数αｌｍを求めることができる。微分を実行すれば、
【数１０】
2004118618.tif 000014
となる。ここで、ｌ’＝１，２，…，Ｌである。この連立一次方程式を解くことによって、残差最小化係数αｌｍを求めることができる。
【００３４】
例えば、Ｌ＝３のときは、上式は、
【数１１】
2004118618.tif 000015
のような３元一次連立方程式になる。そして、
2004118618.tif 000016
となる。
【００３５】
＜ベクトル列のサンプリング＞
すでに述べたとおり、残差最小化係数αｌｍを求めるためには、上で示したような連立一次方程式を解く必要がある。ところで、上式では、Ａ・φｍ−１　，Ａ・φｍ−２　，…といったベクトル列が頻繁に現れる。このため、本実施形態では、ステップＳ７において、求めたベクトル列Ａ・φｍ　の要素をメモリに保存しておくことにする。これにより、修正近似値φｍ　を求めるために用いる残差最小化係数αｌｍ（ｌ＝１，…，Ｌ）を、すでにメモリに保存されているベクトル列Ａ・φｋ　（ｋ＝ｍ−Ｌ＋１，…，ｍ−１）の要素を用いて求めることができるので、計算効率を向上させることができる。
【００３６】
さらに、本実施形態では、ベクトル列Ａ・φｍ　の保存方法を工夫することによって、メモリ容量と、計算時間とを節約する方法を提案する。すなわち、ベクトル列Ａ・φｍ　の全ての要素ｂ１　，ｂ２　，…，ｂＮ　をメモリに保存するのではなく、一部の要素をサンプリングして保存するようにする。
【００３７】
ここで、集合｛１，２，…，Ｎ｝の部分集合Ωを、次のように定義する。
【数１】
2004118618.tif 000017
ｌｇは整数パラメータであり、βは実数パラメータである。ａｉｉは行列Ａのｉ行ｉ列の要素すなわち対角項、ｆｉ　はソース項ｆの要素である。ｌｇ＝１のとき、または、β＝０のときには、Ω≡｛１，２，…，Ｎ｝となる。そして、ベクトル列Ａ・φｍ　の要素を保存する際には、全Ｎ要素ｂ１　，ｂ２　，…，ｂＮ　を保存するのではなく、行番号ｉが部分集合Ωに含まれる要素のみ、すなわち要素ｂｉ　（ｉ∈Ω）のみを選択して保存する。言い換えると、上のような内積計算を、サンプリング要素のみの和によって近似して行う。
【００３８】
図２はベクトル列Ａ・φｍ　の要素のサンプリング例を概念的に示す図である。同図中、（ａ）はソース項ｆの各要素ｆｉ　の分布を示すグラフであり、行列Ａの対角項ａｉｉによって正規化している。（ｂ）は（ａ）の分布に対するサンプリング例であり、例Ａは全要素を選択する場合、例Ｂは所定間隔おき（５要素おき）に選択する場合、そして例Ｃは部分集合Ωに属する要素を選択する本実施形態の場合である。ここでは、実数パラメータβの値は、｜ｆｉ　／ａｉｉ｜の最大値を基準にして設定しており、β＝０．９８ｍａｘ（｜ｆｉ　／ａｉｉ｜）としている。
【００３９】
図２を見ると分かるように、本実施形態のように部分集合Ωに属する要素を選択した場合（例Ｃ）では、空間的に均等なサンプリング（△印）と、物理現象の特徴を踏まえたいわば局所的なサンプリング（○印と◎印）とが、組み合わされている。なお、もちろん、サンプリングの手法はここで述べたものに限られるものではなく、例Ｂのような空間的なサンプリングのみを行ってもよいし、物理現象を踏まえた局所的なサンプリングのみを行ってもよい。例えば、部分集合Ωを次のように定めてもよい。
Ω　＝｛ｉ：ｍｏｄ［ｉ，ｌｇ］＝１｝
Ω　＝｛ｉ：｜ｆｉ　／ａｉｉ｜＞β｝
ただし、計算精度を保ちつつ、メモリ容量の節約や計算時間の短縮を実現するためには、空間的なサンプリングと、物理現象を踏まえた局所的なサンプリングとを組み合わせて行うのが好ましい。
【００４０】
＜本発明の評価＞
本発明の効果を評価するために、実際にコンピュータを用いて、共通の式の数値計算を、残差切除法を用いて、サンプリングを行わない場合（ＣＡＳＥ１）と、上述の部分集合Ωを用いたサンプリングを用いた場合（ＣＡＳＥ２）とで、実行した。対象とした式は、ポアソン方程式を三次元構造格子上で二次精度差分を用いて離散化した連立一次方程式であり、内部ソルバにはＳＯＲ法を用い、残差最小化係数αｌｍの個数は１５個（Ｌ＝１５）とした。また、ｌｇ＝７，β＝０．９８ｍａｘ（｜ｆｉ　／ａｉｉ｜）とした。
【００４１】
図３は評価結果を示している。図３（ｂ）において、縦軸は相対残差‖ｒ‖／‖ｆ‖であり、横軸はＣＰＵ時間である。図３から、本実施形態に係るサンプリングを用いた場合には、サンプリングを行わない場合に比べて、計算時間は約１５％短縮されており、使用されるメモリ容量が約２５％少なくなっていることが分かる。
【００４２】
なお、ここで説明した数値計算方法は、様々な分野に利用可能である。例えば、流体解析のような解析技術に利用できることはいうまでもなく、車両の走行制御のような制御技術にも利用できる。
【００４３】
【発明の効果】
以上のように本発明によると、修正近似値φｍ　を求めるために用いる残差最小化係数αｌｍ（ｌ＝１，…，Ｌ）が、メモリに保存されているベクトル列Ａ・φｋ　（ｋ＝ｍ−Ｌ＋１，…，ｍ−１）の要素を用いて近似的に求められるので、計算効率を向上させることができる。しかも、ベクトル列Ａ・φｍ　の要素がサンプリングされてメモリに保存されるので、計算に必要なメモリ容量を抑えるとともに、計算時間を短縮することができる。
【図面の簡単な説明】
【図１】残差切除法を用いた数値計算方法を示すフローチャートである。
【図２】本実施形態に係るベクトル列要素のサンプリングを概念的に示す図である。
【図３】本発明の評価結果を示す図である。
【図１】
2004118618.tif 000021
【図２】
2004118618.tif 000022
【図３】
2004118618.tif 000023
【手続補正書】
【提出日】平成15年11月17日(2003.11.17)
【手続補正１】
【補正対象書類名】明細書
【補正対象項目名】特許請求の範囲
【補正方法】変更
【補正の内容】
【特許請求の範囲】
【請求項１】
物理量Ｕが満たすべき偏微分方程式を離散化したときの係数行列（Ｎ行Ｎ列、Ｎは正の整数）をＡ、非斉次項（ソース項）をｆとしたとき、コンピュータを用いて、
Ａ・Ｕ＝ｆ
を解き、物理量Ｕの数値計算を行う方法であって、
物理量Ｕの初期値Ｕ０を設定し、
繰り返し回数ｍに初期値として０を、摂動量φの初期値として０を与え、残差ｒの初期値ｒ０として（ｆ−Ａ・Ｕ０）を設定し、
繰り返し回数ｍをインクリメントしながら、
Ａ・φ＝ｒｍの予測近似値Ψｍを、内部ソルバを有する第１演算部によって、反復計算によって求める第１のステップと、
残差ｒｍ＋１のＬ２ノルムを最小とする修正近似値φｍを、第２演算部によって、最適化ルーチンによって、予測近似値Ψｍから求め、近似解Ｕｍ＋１として（Ｕｍ＋φｍ）を与え、残差ｒｍ＋１として（ｒｍ−Ａ・φｍ）を与える第２のステップとを、近似解Ｕｍが収束するまで、繰り返し実行するものであり、
前記第２のステップは、
求めたベクトル列Ａ・φｍ　の要素を、所定のサンプリング方法によってサンプリングして、メモリに保存するものであり、かつ、
修正近似値φｍ　を求めるために用いる残差最小化係数αｌｍ（ｌ＝１，…，Ｌ）を、前記メモリに保存されているベクトル列Ａ・φｋ　（ｋ＝ｍ−Ｌ＋１，…，ｍ−１）の要素を用いて、近似的に求める
ことを特徴とする数値計算方法。
【請求項２】
請求項１記載の数値計算方法において、
前記第２のステップは、
ベクトル列Ａ・φｍ　の要素ｂ１　，ｂ２　，…，ｂＮ　のサンプリングにおいて、要素ｂｉ（ｉ∈Ω）を選択する（ただし、
【数１】
2004118618.tif 000024
ｌｇは整数、βは実数、ｆｉ　はソース項ｆの要素、ａｉｉは行列Ａのｉ行ｉ列の対角項）
ことを特徴とする数値計算方法。
【請求項３】
物理量Ｕが満たすべき偏微分方程式を離散化したときの係数行列（Ｎ行Ｎ列、Ｎは正の整数）をＡ、非斉次項（ソース項）をｆとしたとき、
Ａ・Ｕ＝ｆ
を解き、物理量Ｕの数値計算を行う装置であって、
物理量Ｕの初期値Ｕ０を設定し、
繰り返し回数ｍに初期値として０を、摂動量φの初期値として０を与え、残差ｒの初期値ｒ０として（ｆ−Ａ・Ｕ０）を設定し、
繰り返し回数ｍをインクリメントしながら、
Ａ・φ＝ｒｍの予測近似値Ψｍを、内部ソルバを有する第１演算部によって、反復計算によって求める第１のステップと、
残差ｒｍ＋１のＬ２ノルムを最小とする修正近似値φｍを、第２演算部によって、最適化ルーチンによって、予測近似値Ψｍから求め、近似解Ｕｍ＋１として（Ｕｍ＋φｍ）を与え、残差ｒｍ＋１として（ｒｍ−Ａ・φｍ）を与える第２のステップとを、近似解Ｕｍが収束するまで、繰り返し実行するものであり、
前記第２のステップは、
求めたベクトル列Ａ・φｍ　の要素を、所定のサンプリング方法によってサンプリングして、メモリに保存するものであり、かつ、
修正近似値φｍ　を求めるために用いる残差最小化係数αｌｍ（ｌ＝１，…，Ｌ）を、前記メモリに保存されているベクトル列Ａ・φｋ　（ｋ＝ｍ−Ｌ＋１，…，ｍ−１）の要素を用いて、近似的に求める
ことを特徴とする数値計算装置。
【請求項４】
請求項３記載の数値計算装置において、
前記第２のステップは、
ベクトル列Ａ・φｍ　の要素ｂ１　，ｂ２　，…，ｂＮ　のサンプリングにおいて、要素ｂｉ（ｉ∈Ω）を選択する（ただし、
【数１】
2004118618.tif 000025
ｌｇは整数、βは実数、ｆｉ　はソース項ｆの要素、ａｉｉは行列Ａのｉ行ｉ列の対角項）
ことを特徴とする数値計算装置。
【請求項５】
コンピュータに、物理量Ｕが満たすべき偏微分方程式を離散化したときの係数行列（Ｎ行Ｎ列、Ｎは正の整数）をＡ、非斉次項（ソース項）をｆとしたとき、
Ａ・Ｕ＝ｆ
を解かせ、物理量Ｕの数値計算を行う数値計算用プログラムであって、
物理量Ｕの初期値Ｕ０を設定し、
繰り返し回数ｍに初期値として０を、摂動量φの初期値として０を与え、残差ｒの初期値ｒ０として（ｆ−Ａ・Ｕ０）を設定し、
繰り返し回数ｍをインクリメントしながら、
Ａ・φ＝ｒｍの予測近似値Ψｍを、内部ソルバを有する第１演算部によって、反復計算によって求める第１のステップと、
残差ｒｍ＋１のＬ２ノルムを最小とする修正近似値φｍを、第２演算部によって、最適化ルーチンによって、予測近似値Ψｍから求め、近似解Ｕｍ＋１として（Ｕｍ＋φｍ）を与え、残差ｒｍ＋１として（ｒｍ−Ａ・φｍ）を与える第２のステップとを、近似解Ｕｍが収束するまで、繰り返し実行するものであり、
前記第２のステップは、
求めたベクトル列Ａ・φｍ　の要素を、所定のサンプリング方法によってサンプリングして、メモリに保存するものであり、かつ、
修正近似値φｍ　を求めるために用いる残差最小化係数αｌｍ（ｌ＝１，…，Ｌ）を、前記メモリに保存されているベクトル列Ａ・φｋ　（ｋ＝ｍ−Ｌ＋１，…，ｍ−１）の要素を用いて、近似的に求める
処理をコンピュータに実行させる数値計算用プログラム。
【請求項６】
請求項５記載の数値計算用プログラムにおいて、
前記第２のステップは、
ベクトル列Ａ・φｍ　の要素ｂ１　，ｂ２　，…，ｂＮ　のサンプリングにおいて、要素ｂｉ（ｉ∈Ω）を選択する（ただし、
【数１】
2004118618.tif 000026
ｌｇは整数、βは実数、ｆｉ　はソース項ｆの要素、ａｉｉは行列Ａのｉ行ｉ列の対角項）
ことを特徴とする数値計算用プログラム。
【手続補正２】
【補正対象書類名】明細書
【補正対象項目名】０００６
【補正方法】変更
【補正の内容】
【０００６】
【課題を解決するための手段】
前記の問題を解決するために、請求項１の発明が講じた解決手段は、物理量Ｕが満たすべき偏微分方程式を離散化したときの係数行列（Ｎ行Ｎ列、Ｎは正の整数）をＡ、非斉次項（ソース項）をｆとしたとき、コンピュータを用いて、Ａ・Ｕ＝ｆを解き、物理量Ｕの数値計算を行う方法として、物理量Ｕの初期値Ｕ０を設定し、繰り返し回数ｍに初期値として０を、摂動量φの初期値として０を与え、残差ｒの初期値ｒ０として（ｆ−Ａ・Ｕ０）を設定し、繰り返し回数ｍをインクリメントしながら、Ａ・φ＝ｒｍの予測近似値Ψｍを内部ソルバを有する第１演算部によって反復計算によって求める第１のステップと、残差ｒｍ＋１のＬ２ノルムを最小とする修正近似値φｍを、第２演算部によって最適化ルーチンによって、予測近似値Ψｍから求め、近似解Ｕｍ＋１として（Ｕｍ＋φｍ）を与え、残差ｒｍ＋１として（ｒｍ−Ａ・φｍ）を与える第２のステップとを、近似解Ｕｍが収束するまで繰り返し実行するものであり、前記第２のステップは、求めたベクトル列Ａ・φｍ　の要素を、所定のサンプリング方法によってサンプリングして、メモリに保存するものであり、かつ、修正近似値φｍ　を求めるために用いる残差最小化係数αｌｍ（ｌ＝１，…，Ｌ）を、前記メモリに保存されているベクトル列Ａ・φｋ　（ｋ＝ｍ−Ｌ＋１，…，ｍ−１）の要素を用いて近似的に求めるものである。


【手続補正３】
【補正対象書類名】明細書
【補正対象項目名】０００８
【補正方法】変更
【補正の内容】
【０００８】
また、請求項３の発明が講じた解決手段は、物理量Ｕが満たすべき偏微分方程式を離散化したときの係数行列（Ｎ行Ｎ列、Ｎは正の整数）をＡ、非斉次項（ソース項）をｆとしたとき、Ａ・Ｕ＝ｆを解き、物理量Ｕの数値計算を行う装置として、物理量Ｕの初期値Ｕ０を設定し、繰り返し回数ｍに初期値として０を、摂動量φの初期値として０を与え、残差ｒの初期値ｒ０として（ｆ−Ａ・Ｕ０）を設定し、繰り返し回数ｍをインクリメントしながら、Ａ・φ＝ｒｍの予測近似値Ψｍを内部ソルバを有する第１演算部によって反復計算によって求める第１のステップと、残差ｒｍ＋１のＬ２ノルムを最小とする修正近似値φｍを第２演算部によって最適化ルーチンによって予測近似値Ψｍから求め、近似解Ｕｍ＋１として（Ｕｍ＋φｍ）を与え、残差ｒｍ＋１として（ｒｍ−Ａ・φｍ）を与える第２のステップとを、近似解Ｕｍが収束するまで繰り返し実行するものであり、前記第２のステップは、求めたベクトル列Ａ・φｍ　の要素を、所定のサンプリング方法によってサンプリングしてメモリに保存するものであり、かつ、修正近似値φｍ　を求めるために用いる残差最小化係数αｌｍ（ｌ＝１，…，Ｌ）を、前記メモリに保存されているベクトル列Ａ・φｋ　（ｋ＝ｍ−Ｌ＋１，…，ｍ−１）の要素を用いて近似的に求めるものである。


【手続補正４】
【補正対象書類名】明細書
【補正対象項目名】００１０
【補正方法】変更
【補正の内容】
【００１０】
また、請求項５の発明が講じた解決手段は、コンピュータに、物理量Ｕが満たすべき偏微分方程式を離散化したときの係数行列（Ｎ行Ｎ列、Ｎは正の整数）をＡ、非斉次項（ソース項）をｆとしたとき、Ａ・Ｕ＝ｆを解かせ、物理量Ｕの数値計算を行う数値計算用プログラムとして、物理量Ｕの初期値Ｕ０を設定し、繰り返し回数ｍに初期値として０を、摂動量φの初期値として０を与え、残差ｒの初期値ｒ０として（ｆ−Ａ・Ｕ０）を設定し、繰り返し回数ｍをインクリメントしながら、Ａ・φ＝ｒｍの予測近似値Ψｍを内部ソルバを有する第１演算部によって反復計算によって求める第１のステップと、残差ｒｍ＋１のＬ２ノルムを最小とする修正近似値φｍを、第２演算部によって最適化ルーチンによって予測近似値Ψｍから求め、近似解Ｕｍ＋１として（Ｕｍ＋φｍ）を与え、残差ｒｍ＋１として（ｒｍ−Ａ・φｍ）を与える第２のステップとを、近似解Ｕｍが収束するまで繰り返し実行するものであり、前記第２のステップは、求めたベクトル列Ａ・φｍ　の要素を、所定のサンプリング方法によってサンプリングしてメモリに保存するものであり、かつ、修正近似値φｍ　を求めるために用いる残差最小化係数αｌｍ（ｌ＝１，…，Ｌ）を、前記メモリに保存されているベクトル列Ａ・φｋ　（ｋ＝ｍ−Ｌ＋１，…，ｍ−１）の要素を用いて近似的に求める処理をコンピュータに実行させるものである。


【手続補正５】
【補正対象書類名】明細書
【補正対象項目名】００２０
【補正方法】変更
【補正の内容】
【００２０】
＜最適化制御ルーチン＞
いま、繰り返し回数をｍとしたとき、（ｍ＋１）回目の残差ｒｍ＋１のＬ２ノルムを最小とする合成摂動量φｍと新しい近似解Ｕｍ＋１を次のように定義する。
【数式２】
2004118618.tif 000027
Ｕｍ＋１＝Ｕｍ＋φｍ　…（６）
αｌｍ（ｌ＝１，２，３，…，Ｌ）は残差最小化係数であり、後述する計算方法によって求まる定数である。残差ｒｍ＋１のＬ２ノルムがより小さくなるように、次のように残差の最小化を行う。


【手続補正６】
【補正対象書類名】明細書
【補正対象項目名】００３３
【補正方法】変更
【補正の内容】
【００３３】
残差最小化係数αｌｍは、例えば最小自乗法を用いて、（ｍ＋１）回目の残差ノルム‖ｒｍ＋１‖（自乗和の平方根）を最小にするという条件を与えることによって、求められる。
【数式８】
2004118618.tif 000028
これをαｌｍで微分すると、
【数式９】
2004118618.tif 000029
この式をαｌｍについて解くことによって、残差最小化係数αｌｍを求めることができる。微分を実行すれば、
【数式１０】
2004118618.tif 000030
となる。ここで、ｌ´＝１，２，…，Ｌである。この連立一次方程式を解くことによって、残差最小化係数αｌｍを求めることができる。