逆ラプラス変換プログラム、逆ラプラス変換のためのテーブル作成プログラム、逆ラプラス変換の数値解算出プログラム、および逆ラプラス変換装置
【課題】原像の導函数が増大する場合、または原像が原点で零でない場合においても、逆ラプラス変換の数値解を計算することができる逆ラプラス変換プログラム、逆ラプラス変換のためのテーブル作成プログラム、逆ラプラス変換の数値解算出プログラム、および逆ラプラス変換装置を提供する。
【解決手段】テーブル作成部4は、原点で零であり、かつ絶対連続な函数からなる重み付き再生核ヒルベルト空間上で、重み付き二乗可積分空間を観測空間としたチホノフ正則化法により導かれる第二種積分方程式を離散化して得られる連立方程式の解を求め、連立方程式の解に基づく第二種積分方程式の数値解を含む情報を記述したHテーブルを作成する。逆変換部5は、Hテーブルを参照して、第二種積分方程式の数値解と軟化子函数を乗じたラプラス変換像との重み付き二乗可積分空間での内積を数値計算で求める。
【解決手段】テーブル作成部4は、原点で零であり、かつ絶対連続な函数からなる重み付き再生核ヒルベルト空間上で、重み付き二乗可積分空間を観測空間としたチホノフ正則化法により導かれる第二種積分方程式を離散化して得られる連立方程式の解を求め、連立方程式の解に基づく第二種積分方程式の数値解を含む情報を記述したHテーブルを作成する。逆変換部5は、Hテーブルを参照して、第二種積分方程式の数値解と軟化子函数を乗じたラプラス変換像との重み付き二乗可積分空間での内積を数値計算で求める。
【発明の詳細な説明】
【技術分野】
【0001】
本発明は、逆ラプラス変換プログラム、逆ラプラス変換のためのテーブル作成プログラム、逆ラプラス変換の数値解算出プログラム、および逆ラプラス変換装置に関し、特に再生核と正則化法用いて逆ラプラス変換の数値解を算出する技術に関する。
【背景技術】
【0002】
逆ラプラス変換は、電気電子工学、油田の調査、画像処理をはじめとする様々な分野で広く応用されている。従来、コンピュータを用いて逆ラプラス変換の数値解を算出するためには、複素積分の数値計算による方法や、ラプラス変換表を用いた方法が用いられていた。しかしながら、このような従来の方法では、丸め誤差または離散化誤差などの計算誤差が生じ、数値的に不安定(ill-conditioned)であるという問題がある。
【0003】
このような数値的不安定性を解消する手法として、再生核と正則化法に基づいて逆ラプラス変換の数値解を求める方法が提案されている(たとえば、非特許文献1を参照)。この方法では、逆ラプラス変換の数値解が、再生核ヒルベルト空間での最良近似函数として求めることができる。
【非特許文献1】松浦他、"Numerical Real Inversion formulas of the Laplace Transform by Using a Fredholm Integral Equation of the Second Kind"、日本数学会、2007年度年会、応用数学分科会、講演アブストラクト、p.69-72
【発明の開示】
【発明が解決しようとする課題】
【0004】
しかしながら、上記の再生核と正則化法を用いる手法では、原像が原点で零、かつ原像の導函数が増大しない絶対連続函数からなる再生核空間を考えており、原像がこれら2つの条件を満たす場合に対してのみ適用可能である。
【0005】
それゆえに、本発明の目的は、原像の導函数が増大する場合、または原像が原点で零でない場合においても、再生核と正則化法を用いて、逆ラプラス変換の数値解を計算することができる逆ラプラス変換プログラム、逆ラプラス変換のためのテーブル作成プログラム、逆ラプラス変換の数値解算出プログラム、および逆ラプラス変換装置を提供することである。
【課題を解決するための手段】
【0006】
本発明のある局面に係る逆ラプラス変換プログラムは、与えられた正則化パラメータと軟化子パラメータのもとで、ラプラス変換像に対する原像の数値解を算出するために、コンピュータに、原点で零であり、かつ絶対連続な函数からなる重み付き再生核ヒルベルト空間上で、重み付き二乗可積分空間を観測空間としたチホノフ正則化法により導かれる第二種積分方程式を離散化して得られる連立方程式の解を求め、連立方程式の解に基づく第二種積分方程式の数値解を含む情報を記述したテーブルを作成するステップと、作成したテーブルを記憶部に保存するステップと、記憶部のテーブルを参照して、第二種積分方程式の数値解と軟化子函数を乗じたラプラス変換像との重み付き二乗可積分空間での内積を数値計算で求めることにより、原像の数値解を算出するステップと、原像の数値解を出力するステップとを実行させる。
【0007】
好ましくは、正則化パラメータをaとし、軟化子パラメータsとし、任意の関数gのラプラス変換像をLgとし、数値解を算出する原像に対するラプラス変換像をF(p)とし、軟化子函数をRs(p)としたときに、原点で零であり、かつ絶対連続な函数fからなる重み付き再生核ヒルベルト空間Hkのノルムは、式(B1)で表わされ、
【0008】
【数1】
【0009】
重み付き再生核ヒルベルト空間Hkでの再生核K(t,t′)が式(B2)で表わされ、
【0010】
【数2】
【0011】
重み付き二乗可積分空間は、L2(R+,u(p)dp)で表わされ、ρ(t)とu(p)の間には、式(B3)の条件が成立し、
【0012】
【数3】
【0013】
第二種積分方程式は、式(B4)で表わされ、
【0014】
【数4】
【0015】
内積は、式(B5)で表わされる。
【0016】
【数5】
【0017】
好ましくは、ρ(t)は、式(B6)で表わされ、u(p)は、式(B7)で表わされ、Rs(p)は、式(B8)で表わされる。
【0018】
【数6】
【0019】
好ましくは、コンピュータは、Kビット長(Kは自然数)の汎用レジスタと、直前の演算により繰上りまたは繰下がりが生じたか否かを示すフラグを格納するフラグレジスタと、フラグを参照して、直前の演算により繰上りまたは繰下がりが生じたか否かに応じた演算を行なう演算器とを備え、テーブルを作成するステップおよび原像の数値解を算出するステップは、数値計算の対象となる変数について、その仮数部をK×Nビット(Nは2以上の自然数)の多倍長で表わし、変数の仮数部をKビットごとに区切って、区切られた下位側から順次、汎用レジスタを用いて演算器に演算を行なわせるステップを含む。
【0020】
好ましくは、連立方程式は、Ax=bの形で表わされ、Aは各要素がtに依存しない係数行列であり、xは連立方程式の解を要素とするベクトルであり、bは各要素がtに依存するベクトルであり、テーブルを作成するステップは、行列Aの逆行列A-1を算出して、逆行列A-1の要素を記憶部に記憶するステップと、計算範囲内のすべてのtについて、記憶部に記憶されている逆行列A-1の要素の値に基づいて、連立方程式の解を算出するステップとを含む。
【0021】
好ましくは、連立方程式は、Ax=bの形で表わされ、Aは各要素がtに依存しない係数行列であり、xは連立方程式の解を要素とするベクトルであり、bは各要素がtに依存するベクトルであり、テーブルを作成するステップは、行列Aを上三角行列と下三角行列の積に分解して、上三角行列と下三角行列の要素を記憶部に記憶するステップと、計算範囲内のすべてのtについて、記憶部に記憶されている上三角行列と下三角行列の要素の値に基づいて、連立方程式の解を算出するステップとを含む。
【0022】
好ましくは、テーブルを作成するステップは、正則化パラメータaおよび軟化子パラメータsと、その他のパラメータn、U、Lとから、式(B9)〜(B14)に従って、それぞれ変数h、xj、pj、qj、aij、H(pj, t)を計算するステップと、
【0023】
【数7】
【0024】
連立方程式は、式(B15)で表わされ、式(B16)で表わされる行列Aをコレスキー分解し、式(B15)の解を算出するステップと、
【0025】
【数8】
【0026】
式(B17)に従って、式(B15)の解から第二種積分方程式の数値解{Hin,a,t:i=0,1,2,・・・,n}を算出し、第二種積分方程式の数値解を含む情報を記述したテーブルを作成するステップとを含み、
【0027】
【数9】
【0028】
原像の数値解を算出するステップは、テーブルを参照して、式(B18)に従って、ラプラス変換像F(pj)から、原像の数値解fa,s(n)(t)を算出するステップを含む。
【0029】
【数10】
【0030】
本発明のある局面に係るテーブル作成プログラムは、与えられた正則化パラメータと軟化子パラメータのもとで、ラプラス変換像に対する原像の数値解を算出するためのテーブルを作成するために、コンピュータに、原点で零であり、かつ絶対連続な函数からなる重み付き再生核ヒルベルト空間上で、重み付き二乗可積分空間を観測空間としたチホノフ正則化法により導かれる第二種積分方程式を離散化して得られる連立方程式の解を求め、連立方程式の解に基づく第二種積分方程式の数値解を含む情報を記述したテーブルを作成するステップと、作成したテーブルを記憶部に保存するステップとを実行させる。
【0031】
好ましくは、正則化パラメータをaとし、軟化子パラメータsとし、任意の関数gのラプラス変換像をLgとし、原点で零であり、かつ絶対連続な函数fからなる重み付き再生核ヒルベルト空間Hkのノルムは式(B19)で表わされ、
【0032】
【数11】
【0033】
重み付き再生核ヒルベルト空間Hkでの再生核K(t,t′)が式(B20)で表わされ、
【0034】
【数12】
【0035】
重み付き二乗可積分空間は、L2(R+,u(p)dp)で表わされ、ρ(t)とu(p)の間には、式(B21)の条件が成立し、
【0036】
【数13】
【0037】
第二種積分方程式は、式(B22)で表わされる。
【0038】
【数14】
【0039】
好ましくは、ρ(t)は、式(B23)で表わされ、u(p)は、式(B24)で表わされる。
【0040】
【数15】
【0041】
好ましくは、コンピュータは、Kビット長(Kは自然数)の汎用レジスタと、直前の演算により繰上りまたは繰下がりが生じたか否かを示すフラグを格納するフラグレジスタと、フラグを参照して、直前の演算により繰上りまたは繰下がりが生じたか否かに応じた演算を行なう演算器とを備え、テーブルを作成するステップは、数値計算の対象となる変数について、その仮数部をK×Nビット(Nは2以上の自然数)の多倍長で表わし、変数の仮数部をKビットごとに区切って、区切られた下位側から順次、汎用レジスタを用いて演算器に演算を行なわせるステップを含む。
【0042】
好ましくは、連立方程式は、Ax=bの形で表わされ、Aは各要素がtに依存しない係数行列であり、xは連立方程式の解を要素とするベクトルであり、bは各要素がtに依存するベクトルであり、テーブルを作成するステップは、行列Aの逆行列A-1を算出して、逆行列A-1の要素を記憶部に記憶するステップと、計算範囲内のすべてのtについて、記憶部に記憶されている逆行列A-1の要素の値に基づいて、連立方程式の解を算出するステップとを含む。
【0043】
好ましくは、連立方程式は、Ax=bの形で表わされ、Aは各要素がtに依存しない係数行列であり、xは連立方程式の解を要素とするベクトルであり、bは各要素がtに依存するベクトルであり、テーブルを作成するステップは、行列Aを上三角行列と下三角行列の積に分解して、上三角行列と下三角行列の要素を記憶部に記憶するステップと、計算範囲内のすべてのtについて、記憶部に記憶されている上三角行列と下三角行列の要素の値に基づいて、連立方程式の解を算出するステップとを含む。
【0044】
好ましくは、テーブルを作成するステップは、正則化パラメータaおよび軟化子パラメータsと、その他のパラメータn、U、Lとから、式(B25)〜(B30)に従って、それぞれ変数h、xj、pj、qj、aij、H(pj, t)を計算するステップと、
【0045】
【数16】
【0046】
連立方程式は、式(B31)で表わされ、式(B32)で表わされる行列Aをコレスキー分解し、式(B31)の解を算出するステップと、
【0047】
【数17】
【0048】
式(B33)に従って、式(B31)の解から第二種積分方程式の数値解{Hin,a,t:i=0,1,2,・・・,n}を算出し、第二種積分方程式の数値解を含む情報を記述したテーブルを作成するステップとを含む。
【0049】
【数18】
【0050】
本発明のある局面に係る逆ラプラス変換の数値解算出プログラムは、上記記載の逆ラプラス変換のためのテーブル作成プログラムによって作成されたテーブルを用いて、ラプラス変換像に対する原像の数値解を作成するために、コンピュータに、上記記載の逆ラプラス変換のためのテーブル作成プログラムで作成されたテーブルが入力され、入力されたテーブルを記憶部に保存するステップと、記憶部のテーブルを参照して、第二種積分方程式の数値解と軟化子函数を乗じたラプラス変換像との重み付き二乗可積分空間での内積を数値計算で求めることにより、原像の数値解を算出するステップと、原像の数値解を出力するステップとを実行させる。
【0051】
本発明のある局面に係る逆ラプラス変換の数値解算出プログラムは、上記記載の逆ラプラス変換のためのテーブル作成プログラムによって作成されたテーブルを用いて、ラプラス変換像に対する原像の数値解を作成するために、コンピュータに、上記記載の逆ラプラス変換のためのテーブル作成プログラムで作成されたテーブルが入力され、入力されたテーブルを記憶部に保存するステップと、記憶部のテーブルを参照して、第二種積分方程式の数値解と軟化子函数Rs(p)を乗じたラプラス変換像F(p)との重み付き二乗可積分空間での内積を数値計算で求めることにより、原像の数値解を算出するステップと、内積は、式(B34)で表わされ、
【0052】
【数19】
【0053】
原像の数値解を出力するステップとを実行させる。
本発明のある局面に係る逆ラプラス変換の数値解算出プログラムは、上記記載の逆ラプラス変換のためのテーブル作成プログラムによって作成されたテーブルを用いて、ラプラス変換像に対する原像の数値解を作成するために、コンピュータに、上記記載の逆ラプラス変換のためのテーブル作成プログラムで作成されたテーブルが入力され、入力されたテーブルを記憶部に保存するステップと、記憶部のテーブルを参照して、第二種積分方程式の数値解と軟化子函数Rs(p)を乗じたラプラス変換像F(p)との重み付き二乗可積分空間での内積を数値計算で求めることにより、原像の数値解を算出するステップと、軟化子函数Rs(p)は、式(B35)で表わされ、
【0054】
【数20】
【0055】
内積は、式(B36)で表わされ、
【0056】
【数21】
【0057】
原像の数値解を出力するステップとを実行させる。
本発明のある局面に係る逆ラプラス変換の数値解算出プログラムは、上記記載の逆ラプラス変換のためのテーブル作成プログラムによって作成されたテーブルを用いて、逆ラプラス変換の数値解を作成するために、Kビット長(Kは自然数)の汎用レジスタと、直前の演算により繰上りまたは繰下がりが生じたか否かを示すフラグを格納するフラグレジスタと、フラグを参照して、直前の演算により繰上りまたは繰下がりが生じたか否かに応じた演算を行なう演算器とを備えたコンピュータに、上記記載の逆ラプラス変換のためのテーブル作成プログラムで作成されたテーブルが入力され、入力されたテーブルを記憶部に保存するステップと、記憶部のテーブルを参照して、第二種積分方程式の数値解と軟化子函数Rs(p)を乗じたラプラス変換像F(p)との重み付き二乗可積分空間での内積を数値計算で求めることにより、原像の数値解を算出するステップと、軟化子函数Rs(p)は、式(B37)で表わされ、
【0058】
【数22】
【0059】
内積は、式(B38)で表わされ、
【0060】
【数23】
【0061】
原像の数値解を出力するステップとを実行させ、原像の数値解を算出するステップは、数値計算の対象となる変数について、その仮数部をK×Nビット(Nは2以上の自然数)の多倍長で表わし、変数の仮数部をKビットごとに区切って、区切られた下位側から順次、汎用レジスタを用いて演算器に演算を行なわせるステップを含む。
【0062】
本発明のある局面に係る逆ラプラス変換の数値解算出プログラムは、上記記載の逆ラプラス変換のためのテーブル作成プログラムによって作成されたテーブルを用いて、逆ラプラス変換の数値解を作成するために、コンピュータに、上記記載の逆ラプラス変換のためのテーブル作成プログラムで作成されたテーブルが入力され、入力されたテーブルを記憶部に保存するステップと、記憶部のテーブルを参照して、式(B39)に従って、第二種積分方程式の数値解とラプラス変換像F(pj)から原像の数値解fa,s(n)(t)を算出するステップと、
【0063】
【数24】
【0064】
原像の数値解を出力するステップとを実行させる。
本発明のある局面に係る逆ラプラス変換装置は、与えられた正則化パラメータと軟化子パラメータのもとで、ラプラス変換像に対する原像の数値解を算出する逆ラプラス変換装置であって、原点で零であり、かつ絶対連続な函数からなる重み付き再生核ヒルベルト空間上で、重み付き二乗可積分空間を観測空間としたチホノフ正則化法により導かれる第二種積分方程式を離散化して得られる連立方程式の解を求め、連立方程式の解に基づく第二種積分方程式の数値解を含む情報を記述したテーブルを作成するテーブル作成部と、作成したテーブルを記憶する記憶部と、記憶部のテーブルを参照して、第二種積分方程式の数値解と軟化子函数を乗じたラプラス変換像との重み付き二乗可積分空間での内積を数値計算で求めることにより、原像の数値解を算出する逆変換部と、原像の数値解を出力する出力部とを備える。
【0065】
好ましくは、正則化パラメータをaとし、軟化子パラメータsとし、任意の関数gのラプラス変換像をLgとし、数値解を算出する原像に対するラプラス変換像をF(p)とし、軟化子函数をRs(p)としたときに、原点で零であり、かつ絶対連続な函数fからなる重み付き再生核ヒルベルト空間Hkのノルムは式(B40)で表わされ、
【0066】
【数25】
【0067】
重み付き再生核ヒルベルト空間Hkでの再生核K(t,t′)が式(B41)で表わされ、
【0068】
【数26】
【0069】
重み付き二乗可積分空間は、L2(R+,u(p)dp)で表わされ、ρ(t)とu(p)の間には、式(B42)の条件が成立し、
【0070】
【数27】
【0071】
第二種積分方程式は、式(B43)で表わされ、
【0072】
【数28】
【0073】
内積は、式(B44)で表わされる。
【0074】
【数29】
【0075】
好ましくは、ρ(t)は、式(B45)で表わされ、u(p)は、式(B46)で表わされ、Rs(p)は、式(B47)で表わされる。
【0076】
【数30】
【0077】
テーブル作成部および逆変換部とは、共通して、Kビット長(Kは自然数)の汎用レジスタと、直前の演算により繰上りまたは繰下がりが生じたか否かを示すフラグを格納するフラグレジスタと、フラグを参照して、直前の演算により繰上りまたは繰下がりが生じたか否かに応じた演算を行なう演算器とを備え、演算器は、数値計算の対象となる変数について、その仮数部をK×Nビット(Nは2以上の自然数)の多倍長で表わし、変数の仮数部をKビットごとに区切って、区切られた下位側から順次、汎用レジスタを用いて演算を行なう。
【0078】
好ましくは、連立方程式は、Ax=bの形で表わされ、Aは各要素がtに依存しない係数行列であり、xは連立方程式の解を要素とするベクトルであり、bは各要素がtに依存するベクトルであり、テーブル作成部は、行列Aの逆行列A-1を算出して、逆行列A-1の要素を記憶部に記憶し、計算範囲内のすべてのtについて、記憶部に記憶されている逆行列A-1の要素の値に基づいて、連立方程式の解を算出する。
【0079】
好ましくは、連立方程式は、Ax=bの形で表わされ、Aは各要素がtに依存しない係数行列であり、xは連立方程式の解を要素とするベクトルであり、bは各要素がtに依存するベクトルであり、テーブル作成部は、行列Aを上三角行列と下三角行列の積に分解して、上三角行列と下三角行列の要素を記憶部に記憶し、計算範囲内のすべてのtについて、記憶部に記憶されている上三角行列と下三角行列の要素の値に基づいて、連立方程式の解を算出する。
【0080】
好ましくは、テーブル作成部は、正則化パラメータaおよび軟化子パラメータsと、その他のパラメータn、U、Lとから、式(B48)〜(B53)に従って、それぞれ変数h、xj、pj、qj、aij、H(pj, t)を計算し、
【0081】
【数31】
【0082】
連立方程式は、式(B54)で表わされ、テーブル作成部は、式(B55)で表わされる行列Aをコレスキー分解し、式(B54)の解を算出し、
【0083】
【数32】
【0084】
テーブル作成部は、さらに式(B56)に従って、式(B54)の解から第二種積分方程式の数値解{Hin,a,t:i=0,1,2,・・・,n}を算出し、第二種積分方程式の数値解を含む情報を記述したテーブルを作成し、
【0085】
【数33】
【0086】
逆変換部は、テーブルを参照して、式(B57)に従って、ラプラス変換像F(pj)から原像の数値解fa,s(n)(t)を算出する。
【0087】
【数34】
【発明の効果】
【0088】
本発明によれば、原像の導函数が増大する場合、または原像が原点で零でない場合においても、再生核と正則化法を用いて、逆ラプラス変換の数値解を計算することができる。
【発明を実施するための最良の形態】
【0089】
以下、本発明に係る実施の形態について図面を参照して説明する。
[第1の実施形態]
まず、第1の実施形態における逆ラプラス変換の基本的なアルゴリズムについて説明する。
【0090】
(再生核空間理論を用いた逆ラプラス変換)
ある自然な函数空間の函数f(t)のラプラス変換は、式(1)で表わされる。式(1)の逆変換を考えると、通常はBromwich積分による複素解析的な逆変換公式を考えるが、正の実軸上の値のみを用いた逆変換を求めたい場合がある。このような正の実軸上の値のみを用いた逆変換は実逆ラプラス変換と呼ばれる。また、f(t)は原像、F(P)をラプラス変換像と呼ばれる。
【0091】
【数35】
【0092】
式(2)で表わされる正の実軸R+上の有限ノルムを持ち、f(0)=0を満たす絶対連続な函数fからなる重み付き再生核ヒルベルト空間Hkを考える。この重み付き再生核ヒルベルト空間Hkは、原像の導函数f(t)′が増大する場合も包含している。
【0093】
【数36】
【0094】
この重み付き再生核ヒルベルト空間Hkは、式(3)で表わされる再生核K(t,t′)を有する。
【0095】
【数37】
【0096】
このとき、重み付き再生核ヒルベルト空間Hk上の重み付き二乗可積分空間L2(R+,u(p)dp)への線形写像(Lf)(p)pは有界であって、式(4)が成立する。
【0097】
【数38】
【0098】
式(4)から、一般論に従って、以下が成立する。
チホノフ正則化とは、任意のg∈L2(R+,u(p)dp)と任意の正則化パラメータa>0に対して、式(5)の右辺を最小化することである。重み付き再生核ヒルベルト空間Hk上で、重み付き二乗可積分空間L2(R+,u(p)dp)を観測空間としたチホノフ正則化によって、式(6)で表わされる最良近似函数fa(t)が導かれる。
【0099】
【数39】
【0100】
ここで、Ka(・,t)は、式(7)〜(9)で表わされる。
【0101】
【数40】
【0102】
上述の重み付き再生核ヒルベルト空間Hkは、f(0)=0を満たすことを条件としている。f(0)≠0なる函数f(t)についても、上述の重み付き再生核ヒルベルト空間Hkの特性が利用できるように、軟化子函数Rs(p)を用いて、式(6)の最良近似函数fa(t)を修正する。修正最良近似函数fa,s(t)は、式(10)で表わされる。ここで、sは、軟化子パラメータである。
【0103】
【数41】
【0104】
さて、式(9)のtについてラプラス変換をとり、t′を改めてtとおくと、式(11)が得られる。
【0105】
【数42】
【0106】
式(12)および(13)を用いると、式(11)から式(14)が成立する。式(14)は、書き下すと、式(15)のように表わされる。また、式(10)から、式(16)が成立する。
【0107】
【数43】
【0108】
式(14)は、未知の函数が積分の中に現れており、第二種積分方程式と呼ばれる。修正最良近似函数fa,ε(t)は、第二種積分方程式の解と、軟化子函数を乗じたラプラス変換像との重み付き二乗可積分空間での内積である式(16)で表わされる。本発明の実施形態の逆ラプラス変換では、修正最良近似函数fa,ε(t)を原像f(t)の近似解として求める。
【0109】
(具体的な函数の形)
さて、函数ρ(t)、u(p)、Rs(p)の具体例として、式(17)〜(19)をで表わされるものを用いることができる。
【0110】
【数44】
【0111】
さて、函数ρ(t)が、式(17)で表わされる場合には、式(3)は、式(20)のように表わされる。
【0112】
【数45】
【0113】
さらに、式(20)から、以下の式(21)が成立する。
【0114】
【数46】
【0115】
さらに、式(14)は、式(22)および(23)のように表わされ、式(16)は、式(24)のように表わされる。
【0116】
【数47】
【0117】
(離散化)
式(22)〜(24)をコンピュータで計算するために離散化を行なう。離散化の方法として、ここでは、計算を簡略化するためニストレーム法を用いる。
【0118】
積分区間の上限U、下限L、積分のサンプル数n、正則化パラメータaを用いて、式(25)〜(28)に従って、以下の変数h、xj、pj、wjを定義する。
【0119】
【数48】
【0120】
式(25)〜式(28)を用いて、式(22)を離散化すると式(29)が得られ、式(24)を離散化すると式(30)が得られる。
【0121】
【数49】
【0122】
したがって、式(29)の連立方程式の解{Hi n,a,t:i=0,1,2,・・・,n}を求めて、それらを式(30)に代入することによって、逆ラプラス変換による原像の数値解fa,s(n)(t)が得られる。
【0123】
(式(29)の連立方程式の解法)
式(29)の連立方程式の解を求める方法として、ここでは、以下の方法を用いる。
【0124】
まず、式(31)、(32)に示すように以下の変数qj、aijを用いる。このとき、式(33)、(34)が成り立つ。
【0125】
【数50】
【0126】
これらの変数を用いると、式(29)は、式(35)のように変形される。式(35)を行列の形に書き直すと、式(36)のようになり、係数行列Aを式(37)で定義することができる。
【0127】
【数51】
【0128】
係数行列AをLU分解し、前進代入および後退代入することによって、式(36)の解、つまり式(29)の解{Hin,a,t:i=0,1,2,・・・,n}を得ることができる。また、変数qjを用いると、式(30)は、式(38)のように変形することができる。
【0129】
【数52】
【0130】
(構成)
図1は、本発明の実施形態の逆ラプラス変換装置の構成を表わす図である。
【0131】
図1を参照して、この逆ラプラス変換装置1は、コンピュータを用いて実現されるものであって、入力部2と、メインメモリ6と、CPU3(Central Processing Unit)と、出力部13とを備える。メインメモリ6は、プログラム記憶部7と、変数記憶部8と、Hテーブル記憶部12と、三角行列記憶部9とを含む。CPU3は、テーブル作成部4および逆変換部5として機能する。
【0132】
入力部2は、サンプル数n、上限値U、下限値L、正則化パラメータa、軟化子パラメータs、ラプラス変換像F(pi)などを入力するものであり、たとえばキーボードなどで実現されている。
【0133】
出力部13は、逆ラプラス変換による原像の数値解fa,s(n)(t)を外部、たとえば表示装置14に出力する。表示装置14では、逆ラプラス変換による原像の数値解fa,s(n)(t)が文字またはグラフの形で表示される。
【0134】
プログラム記憶部7は、コンピュータを逆ラプラス変換装置1として機能させるための逆ラプラス変換プログラムが記憶されている。逆ラプラス変換プログラムは、USB(Universal Serial Bus)メモリ、CD−ROM(Compact Disk - Read Only Memory)などのコンピュータ読取り可能な記録媒体であるリムーマブルメディア99を通じて、外部からインストールすることができる。
【0135】
テーブル作成部4は、プログラム記憶部7内の逆ラプラス変換プログラムにしたがって、式(22)および式(23)で表わされる第二種積分方程式を離散化して得られる式(35)の連立方程式の解を数値計算で求め、数値計算に基づいて得られる第二種積分方程式の数値解{Hi n,a,t:i=0,1,2,・・・,n}をHテーブルに記述する。テーブル作成部4は、プログラム記憶部7内の逆ラプラス変換プログラムにしたがって、式(37)で示される係数行列Aを上三角行列Uと下三角行列Lの積に分解して、上三角行列Uと下三角行列Lの要素を三角行列記憶部9に記憶させ、計算範囲内のすべてのtについて、三角行列記憶部9に記憶されている上三角行列Uと下三角行列Lの要素の値に基づいて、式(35)の解を算出する。
【0136】
Hテーブル記憶部12は、テーブル作成部4で作成された図2で示されるようなHテーブルを記憶する。
【0137】
逆変換部5は、Hテーブルを参照して、式(24)を離散化した式(38)を数値計算することにより、逆ラプラス変換による原像の数値解fa,s(n)(t)を算出する。
【0138】
変数記憶部8は、テーブル作成部4によるテーブル作成処理および逆変換部5による逆ラプラス変換による原像の数値解を算出するのに必要な数値計算の変数の値を記憶する。
【0139】
三角行列記憶部9は、テーブル作成部4において、式(35)の連立方程式の解を計算する際に作成された下三角行列Lおよび上三角行列Uを記憶する。
【0140】
(Hテーブルの作成手順)
図3は、第1の実施形態におけるHテーブルを作成する手順を表わすフローチャートである。
【0141】
まず、外部から、入力部2にサンプル数n、上限値U、下限値L、正則化パラメータaが入力される(ステップS101)。
【0142】
次に、テーブル作成部4は、h=(U−L)/nを計算して、hの値を変数記憶部8に保存する(ステップS102)。
【0143】
次に、テーブル作成部4は、j=0、1、2、・・・、nについて、xj=L+j×hを計算して、xjの値を変数記憶部8に保存する(ステップS103)。
【0144】
次に、テーブル作成部4は、j=0、1、2、・・・、nについて、pj=exp{(π/2)×sinh(xj)}を計算して、pjの値を変数記憶部8に保存する(ステップS104)。
【0145】
次に、テーブル作成部4は、j=0、1、2、・・・、nについて、qj=pj×cosh(xj)を計算して、qjの値を変数記憶部8に保存する(ステップS105)。
【0146】
次に、テーブル作成部4は、i=0、1、2、・・・、n、およびj=0、1、2、・・・、nについて、aij=(π/2)×(2/(pi+pj)3)×{1+(pi+pj)+(pi+pj)2/2}×exp(−pj−1/pj)を計算して、aijの値を変数記憶部8に保存する(ステップS106)。
【0147】
次に、テーブル作成部4は、i=0、1、2、・・・、n、およびj=0、1、2、・・・、nについて係数行列Aの成分を計算する。テーブル作成部4は、i=jのときには、係数行列Aのij成分であるa+qi×aiiを計算し、i≠jのときには、係数行列Aのij成分であるqj×aijを計算する(ステップS107)。
【0148】
次に、テーブル作成部4は、係数行列AをLU分解して、上三角行列Uおよび下三角行列Lの各要素の値を三角行列記憶部9に保存する(ステップS108)。
【0149】
次に、テーブル作成部4は、t=t0(初期値)とする(ステップS109)。
次に、テーブル作成部4は、H(pi,t)=(2/pi3)×[1+pi+pi2/2−exp(−tpi)×{1+pi(t+1)+pi2(t+1)2/2}]を計算し、H(pi,t)を変数記憶部8に保存する(ステップS110)。
【0150】
次に、テーブル作成部4は、三角行列記憶部9に記憶されている上三角行列Uおよび下三角行列Lの要素の値をもとに、前進代入および後退代入によって、式(29)の解{Hin,a,t:i=0,1,2,・・・,n}を求めて、変数記憶部8内に記憶する(ステップS111)。
【0151】
次に、テーブル作成部4は、変数記憶部8内の式(29)の解{Hin,a,t:i=0,1,2,・・・,n}をHテーブルに記憶する(ステップS112)。
【0152】
次に、テーブル作成部4は、tがtm(終了値)と等しければ(ステップS113でYES)終了し、tがtm(終了値)と異なるならば(ステップS113でNO)、tをΔtだけ増加させて(ステップS114)、ステップS110からの処理を繰返す。
【0153】
(Hテーブルを用いた逆ラプラス変換の数値解の算出手順)
図4は、第1の実施形態におけるHテーブルを用いて逆ラプラス変換の数値解を算出する手順を表わすフローチャートである。
【0154】
まず、外部から、入力部2に軟化子パラメータsが入力される(ステップS200)。
次に、逆変換部5は、t=t0に設定する(ステップS201)。
【0155】
次に、逆変換部5は、j=0、SUM=0に設定する(ステップS202)。
次に、逆変換部5は、HテーブルからHjn,a,tを読出し、変数記憶部8からpj、qj、hを読み出す(ステップS203)。
【0156】
次に、外部から、入力部2にラプラス変換像F(pj)が入力される(ステップS204)。
【0157】
次に、逆変換部5は、SUM=SUM+(π/2)×h×F(pj)×(1/s2pj2){exp(−spj)−1}2×pj×Hjn,a,t ×qj×exp(−pj−1/pj)を計算し、計算結果を変数記憶部8に保存する(ステップS205)。
【0158】
次に、逆変換部5は、jがnと異なるならば(ステップS206でNO)、jを1だけ増加させて(ステップS207)、ステップS203からの処理を繰返す。
【0159】
逆変換部5は、jがnと等しければ(ステップS206でYES)、逆ラプラス変換値fa,s(n)(t)=SUMとする(ステップS208)。
【0160】
次に、逆変換部5は、tがtmと等しければ(ステップS209でYES)終了し、tがtm(終了値)と異なるならば(ステップS209でNO)、tをΔtだけ増加させて(ステップS210)、ステップS202からの処理を繰返す。
【0161】
以上のように、第1の実施形態の逆ラプラス変換装置によれば、原像の導函数が増大する場合、または原像が原点で零でない場合においても、逆ラプラス変換の数値解を計算することができる。また、式(36)の連立方程式の係数行列Aがtに依存しない性質を利用して、係数行列Aを一度だけLU分解し、分解して得られた下三角行列Lと上三角行列Uの各要素の値を記憶しておき、すべてのtについて、この記憶してある下三角行列Lと上三角行列Uの要素の値を用いて、連立方程式の解を算出するので、tが変わるごとに係数行列AをLU分解するのに比べて計算量を削減することができる。
【0162】
また、連立方程式の解{Hin,a,t:i=0,1,2,・・・,n}は、式(29)を参照すると、ラプラス変換像F(pi)に依存しない。したがって、一度だけ式(36)の解を計算しておき保存しておけば、どのようなF(pi)の逆ラプラス変換を求める場合にも利用することができるという特徴がある。
【0163】
[第1の実施形態の変形例]
本発明の実施形態では、第二種積分方程式を離散化して得られる連立方程式を構成する係数行列AをLU分解することによって連立方程式を解いたが、これに限定するものではなく、係数行列Aの逆行列A-1を計算して、連立方程式を解いてもよい。
【0164】
図5は、第1の実施形態の変形例の逆ラプラス変換装置の構成を表わす図である。
図5を参照して、この逆ラプラス変換装置81は、第1の実施形態のラプラス変換装置1と相違する点は、テーブル作成部84と、逆行列記憶部89である。
【0165】
テーブル作成部84は、プログラム記憶部7内の逆ラプラス変換プログラムにしたがって、式(22)および式(23)で表わされる第二種積分方程式を離散化して得られる式(35)の連立方程式の解を数値計算で求め、数値計算に基づいて得られる第二種積分方程式の数値解{Hi n,a,t:i=0,1,2,・・・,n}をHテーブルに書込む。テーブル作成部84は、プログラム記憶部7内の逆ラプラス変換プログラムにしたがって、式(37)で示される係数行列Aの逆行列A-1を計算して、逆行列A-1の要素の値を逆行列記憶部89に記憶させ、計算範囲内のすべてのtについて、逆行列記憶部89に記憶されている逆行列A-1の要素の値に基づいて、式(35)の解を算出する。
【0166】
逆行列記憶部89は、テーブル作成部4において、式(35)の連立方程式の解を計算する際に作成された逆行列A-1の各要素の値を記憶する。
【0167】
(Hテーブルの作成手順)
図6は、第1の実施形態の変形例おけるHテーブルを作成する手順を表わすフローチャートである。
【0168】
まず、外部から、入力部2にサンプル数n、上限値U、下限値L、正則化パラメータaが入力される(ステップS501)。
【0169】
次に、テーブル作成部4は、h=(U−L)/nを計算して、hの値を変数記憶部8に保存する(ステップS502)。
【0170】
次に、テーブル作成部4は、j=0、1、2、・・・、nについて、xj=L+j×hを計算して、xjの値を変数記憶部8に保存する(ステップS503)。
【0171】
次に、テーブル作成部4は、j=0、1、2、・・・、nについて、pj=exp{(π/2)×sinh(xj)}を計算して、pjの値を変数記憶部8に保存する(ステップS504)。
【0172】
次に、テーブル作成部4は、j=0、1、2、・・・、nについて、qj=pj×cosh(xj)を計算して、qjの値を変数記憶部8に保存する(ステップS505)。
【0173】
次に、テーブル作成部4は、i=0、1、2、・・・、n、およびj=0、1、2、・・・、nについて、aij=(π/2)×(2/(pi+pj)3)×{1+(pi+pj)+(pi+pj)2/2}×exp(−pj−1/pj)を計算して、aijの値を変数記憶部8に保存する(ステップS506)。
【0174】
次に、テーブル作成部4は、i=0、1、2、・・・、n、およびj=0、1、2、・・・、nについて係数行列Aの成分を計算する。テーブル作成部4は、i=jのときには、係数行列Aのij成分であるa+qi×aiiを計算し、i≠jのときには、係数行列Aのij成分であるqj×aijを計算する(ステップS507)。
【0175】
次に、テーブル作成部4は、係数行列Aの逆行列A-1を計算して、逆行列A-1の各要素の値を逆行列記憶部89に保存する(ステップS508)。
【0176】
次に、テーブル作成部4は、t=t0(初期値)とする(ステップS509)。
次に、テーブル作成部4は、H(pi,t)=(2/pi3)×[1+pi+pi2/2−exp(−tpi)×{1+pi(t+1)+pi2(t+1)2/2}]を計算し、H(pi,t)を変数記憶部8に保存する(ステップS510)。
【0177】
次に、テーブル作成部4は、逆行列記憶部89に記憶されている逆行列A-1の要素の値をもとに、式(29)の解{Hin,a,t:i=0,1,2,・・・,n}を求めて、変数記憶部8内に記憶する(ステップS511)。
【0178】
次に、テーブル作成部4は、変数記憶部8内の式(29)の解{Hin,a,t:i=0,1,2,・・・,n}をHテーブルに記憶する(ステップS512)。
【0179】
次に、テーブル作成部4は、tがtm(終了値)と等しければ(ステップS513でYES)終了し、tがtm(終了値)と異なるならば(ステップS513でNO)、tをΔtだけ増加させて(ステップS514)、ステップS510からの処理を繰返す。
【0180】
以上のように、第1の実施形態の変形例の逆ラプラス変換装置によれば、式(36)の連立方程式の係数行列Aがtに依存しない性質を利用して、係数行列Aの逆行列A-1を一度だけ計算して、逆行列A-1の要素の値を記憶しておき、すべてのtについて、この記憶してある逆行列A-1の要素の値を用いて、連立方程式の解を算出するので、tが変わるごとに係数行列Aの逆行列A-1を計算するのに比べて計算量を削減することができる。
【0181】
[第2の実施形態]
第2の実施形態は、式(29)を変形することによって、係数行列Aの代わりに対称行列A′を用いて、第二種積分方程式を離散化して連立方程式の解を算出する逆ラプラス変換装置に関する。
【0182】
(式(23)の連立方程式の解法)
第1の実施形態と同様に、式(31)〜(34)で表わされる変数qj、aijを用いる。
【0183】
第2の実施形態では、式(29)を式(39)のように変形する。
【0184】
【数53】
【0185】
式(40)のように、Hi′n,a,tを定義すると、式(39)は、式(41)のように変形される。式(41)を行列の形に書き直すと、式(42)のようになり、行列A′を式(43)のように定義することができる。
【0186】
【数54】
【0187】
行列A′は、対称行列であり、行列A′をコレスキー分解することによって、式(42)の解を得ることができる。式(42)の解{Hi′n,a,t:i=0,1,2,・・・,n}から、式(44)を用いて、式(29)の解{Hin,a,t:i=0,1,2,・・・,n}が得られる。
【0188】
【数55】
【0189】
(Hテーブルの作成手順)
図7は、第2の実施形態におけるHテーブルを作成する手順を表わすフローチャートである。
【0190】
まず、外部から、入力部2にサンプル数n、上限値U、下限値L、正則化パラメータaが入力される(ステップS301)。
【0191】
次に、テーブル作成部4は、h=(U−L)/nを計算して、hの値を変数記憶部8に保存する(ステップS302)。
【0192】
次に、テーブル作成部4は、j=0、1、2、・・・、nについて、xj=L+j×hを計算して、xjの値を変数記憶部8に保存する(ステップS303)。
【0193】
次に、テーブル作成部4は、j=0、1、2、・・・、nについて、pj=exp{(π/2)×sinh(xj)}を計算して、pjの値を変数記憶部8に保存する(ステップS304)。
【0194】
次に、テーブル作成部4は、j=0、1、2、・・・、nについて、qj=pj×cosh(xj)を計算して、qjの値を変数記憶部8に保存する(ステップS305)。
【0195】
次に、テーブル作成部4は、i=0、1、2、・・・、n、およびj=0、1、2、・・・、nについて、aij=(π/2)×(2/(pi+pj)3)×{1+(pi+pj)+(pi+pj)2/2}×exp(−pj−1/pj)を計算して、aijの値を変数記憶部8に保存する(ステップS306)。
【0196】
次に、テーブル作成部4は、i=0、1、2、・・・、n、およびj=0、1、2、・・・、nについて係数行列Aの成分を計算する。テーブル作成部4は、i=jのときには、係数行列Aのij成分であるa/qi×aiiを計算し、i≠jのときには、係数行列Aのij成分をaijとする(ステップS307)。
【0197】
次に、テーブル作成部4は、係数行列Aをコレスキー分解して、下三角行列L、および上三角行列LT(Tは、転置を表わす)の各要素を三角行列記憶部9に保存する(ステップS308)。
【0198】
次に、テーブル作成部4は、t=t0(初期値)とする(ステップS309)。
次に、テーブル作成部4は、H(pi,t)=(2/pi3)×[1+pi+pi2/2−exp(−tpi)×{1+pi(t+1)+pi2(t+1)2/2}]を計算し、H(pi,t)を変数記憶部8に保存する(ステップS310)。
【0199】
次に、テーブル作成部4は、三角行列記憶部9に記憶されている下三角行列Lおよび上三角行列LTをもとに、前進代入および後退代入によって、式(41)の解{Hi′n,a,t:i=0,1,2,・・・,n}を求めて、変数記憶部8に記憶する(ステップS311)。
【0200】
次に、テーブル作成部4は、i=0、1、2、・・・、nについて、Hi n,a,t=Hi′n,a,t/qiを計算して、変数記憶部8に記憶する(ステップS312)。
【0201】
次に、テーブル作成部4は、変数記憶部8内の式(29)の解である{Hin,a,t:i=0,1,2,・・・,n}をHテーブルに記憶する(ステップS313)。
【0202】
次に、テーブル作成部4は、tがtm(終了値)と等しければ(ステップS314でYES)終了し、tがtm(終了値)と異なるならば(ステップS314でNO)、tをΔtだけ増加させて(ステップS315)、ステップS310からの処理を繰返す。
【0203】
以上のように、第2の実施形態の逆ラプラス変換装置によれば、第1の実施形態と同様の効果に加えて、さらに第二種積分方程式を離散化して得られる連立方程式を構成する係数行列Aを対称行列A′に変換し、コレスキー分解を利用して連立方程式の解を算出することができるので、第1の実施形態よりもさらに計算量を削減できる。
【0204】
[第3の実施形態]
第3の実施形態は、第1の実施形態または第2の実施形態における、連立方程式の解{Hi n,a,t:i=0,1,2,・・・,n}および逆ラプラス変換の数値解fa,s(n)(t)を求める演算を多倍長で行なう装置に関する。以下では、第2の実施形態の演算を多倍長で行なう場合について説明するが、第1の実施形態の演算でも同様である。
【0205】
(構成)
図8は、第3の実施形態のCPU3の具体的な構成を表わす図である。
【0206】
図8を参照して、CPU3は、制御部53と、演算器54と、汎用レジスタ群56と、フラグレジスタ57とを備える。図1では、CPU3は、テーブル作成部4および逆変換部5として機能するとして説明した。図1と図8の関係を説明すると、テーブル作成部4および逆変換部5の機能を具体的に実現するのが、制御部53と、演算器54と、汎用レジスタ群56と、フラグレジスタ57である。つまり、テーブル作成部4および逆変換部5とは、共通に、制御部53と、演算器54と、汎用レジスタ群56と、フラグレジスタ57とを備える。言い換えると、制御部53と、演算器54と、汎用レジスタ群56と、フラグレジスタ57とは、プログラム記憶部7に格納された逆ラプラス変換プログラムにしたがって、第1または第2の実施形態で説明したHテーブルの作成と、逆ラプラス変換の数値解の算出を行なう。
【0207】
演算器54は、加算、減算、乗算、除算、ビットの論理演算、およびビットシフト演算などを行なう。
【0208】
汎用レジスタ群56は、16個の汎用レジスタRAX、RBX、RCX、RDX、RSI、RDI、RBP、RSP、R8、R9、R10、R11、R12、R13、R14、R15を含む。図8において、かっこ内は、アセンブラ言語において指定されるレジスタの名前を表わす。各レジスタは、64ビット長である。
【0209】
図9は、フラグレジスタ57の具体的な構成を表わす図である。
図9を参照して、フラグレジスタ57の最下位ビットに、直前の演算において、汎用レジスタ群56の各汎用レジスタの最上位ビット(64ビット目)において、次の上位ビットへの繰上りまたは繰下がりが生じたか否かを表わすキャリーフラグCFの値が保持される。キャリーフラグCFの値が「1」の場合には、繰上りまたは繰下がりが生じたことを示し、キャリーフラグCFの値が「0」の場合には、繰上りまたは繰下がりが生じなかったことを示す。
【0210】
演算器54による加算には、繰上がり付き加算と繰上りなし加算とがある。X+Yを繰上り加算する場合には、演算器54は、XとYとキャリーフラグCFの値とを加算し、加算結果に繰上りが生じた場合には、キャリーフラグCFの値を「1」にし、加算結果に繰上りが生じなかった場合には、キャリーフラグCFの値を「0」にする。一方、X+Yを繰上りなし加算する場合には、演算器54は、XとYとを加算し、加算結果に繰上りが生じた場合には、キャリーフラグCFの値を「1」にし、加算結果に繰上りが生じなかった場合には、キャリーフラグCFの値を「0」にする。
【0211】
演算器54による減算にも、繰下がり付き減算と繰下がりなし減算とがある。X−Yを繰下がり減算する場合には、演算器54は、XからYおよびキャリーフラグCFの値を減算し、減算結果に繰下がりが生じた場合には、キャリーフラグCFの値を「1」にし、減算結果に繰下がりが生じなかった場合には、キャリーフラグCFの値を「0」にする。一方、X−Yを繰下がりなし減算する場合には、演算器54は、XからYを減算し、減算結果に繰下がりが生じた場合には、キャリーフラグCFの値を「1」にし、減算結果に繰下がりが生じなかった場合には、キャリーフラグCFの値を「0」にする。
【0212】
図8の逆ラプラス変換装置で演算の対象となる数値、たとえば、n、U、L、a、h、s、xj、pj、qj、aij、係数行列Aの要素、係数行列Aをコレスキー分解したLおよびLT、H(pi,t)、Hi′n,a,t、Hin,a,t、F(pj)、S、fa,s(n)(t)などの数値は、浮動小数点形式で表わされる。
【0213】
図10(a)、(b)、(c)は、図8の逆ラプラス変換装置において演算の対象となる数値の表現を説明するための図である。
【0214】
図10(a)に示されるように、数値は、浮動小数点形式で表現される。符号部が「0」の場合は正の数を表わし、符号部が「1」の場合は負の数を表わす。仮数部のみが多倍長64nビット(nは2以上の自然数)で表わされる。指数部および符号部は、64ビットで表わされる。
【0215】
図10(b)、(c)に示されるように、仮数部の整数部の値は、通常は「1」であるが、加算および減算時には、小数点の位置を合わせるために仮数部のビットがシフトされる。仮数部の小数点の位置がKビットだけ左にシフトすると、指数部の値がKだけ減少する。逆に、仮数部の小数点の位置がKビットだけ右にシフトすると、指数部の値がKだけ増加する。
【0216】
図10(b)、(c)に示されるように、汎用レジスタ群56の64ビットの各汎用レジスタを用いて演算を行なうために、メインメモリ6に要素数nの配列aが確保される。仮数部は、64ビットごとに区切られて、区切られた各64ビットが配列aの各要素に格納される。たとえば、最上位から64個のビットが配列の要素a[0]に格納され、次の64個のビットが配列の要素a[1]に格納される。下位側のビットを格納する配列から、順次演算が行なわれる。
【0217】
(本発明の実施形態の加算処理)
本発明の実施形態における仮数部の加算処理について説明する。
【0218】
図11は、フラグレジスタ57を備えるCPU3の一例であるAMD64の加算処理の内容を表わす図である。同図は、C言語からadd(c,a,b,n)で呼び出されるアセンブラルーチンを表わす図である。
【0219】
add(c,a,b,n)は、メインメモリ6内の要素数nの配列aと、メインメモリ6内の要素数nの配列bとを加算して、加算結果をメインメモリ6内の要素数nの配列cに蓄積するルーチンである。
【0220】
add(c,a,b,n)が呼び出されると、制御部53は、カウンタiの値を保持するためのレジスタ%rcxにnの値を設定し、レジスタ%rsiに配列aへのポインタを設定し、レジスタ%rdxに配列bへのポインタを設定し、レジスタ%rdiに配列cへのポインタを設定する。
【0221】
2行目では、制御部53は、最上位ビットでの繰上りの有無を表わす情報を保持するためのレジスタ%raxの値を「0」(つまり、繰上りなし)に設定するとともに、フラグレジスタ57内のキャリーフラグCFの値を「0」(つまり、繰上りなし)に初期化する。
【0222】
6行目では、制御部53は、メインメモリ6内のa[i−1]の値をレジスタ%r8にロードする。
【0223】
7行目では、制御部53は、メインメモリ6内のb[i−1]の値をレジスタ%r10にロードする。
【0224】
8行目では、演算器54は、繰上り付き加算を実行する。すなわち、演算器54は、レジスタ%r8内のデータとレジスタ%r10内のデータと、フラグレジスタ57内のキャリーフラグCFの値とを加算し、加算結果をレジスタ%r10に保存するとともに、加算により繰上りがあった場合に、フラグレジスタ57内のキャリーフラグCFの値を「1」に設定し、繰上りがなかった場合に、フラグレジスタ57内のキャリーフラグCFの値を「0」に設定する。
【0225】
9行目では、制御部53は、レジスタ%r10内の加算結果を、メインメモリ6内のc[i−1]の領域に保存する。
【0226】
10行目では、演算器54は、カウンタiを「1」だけ減らす。
11行目では、制御部53は、カウンタiが「0」の場合には、ループを抜け、カウンタiが「0」以外の場合には、ループの先頭に戻るように制御する。
【0227】
13行目では、演算器54は、即値「0」と、レジスタ%rax内のデータ(=「0」)と、フラグレジスタ57内のキャリーフラグCFの値とを加算し、加算結果をレジスタ%raxに保存する。したがって、レジスタ%raxには、最上位ビットでの繰上りの有無を表わす情報が保持される。
【0228】
(参考:キャリーフラグを利用しない加算処理)
次に、比較のために、フラグレジスタ57を備えないCPU3の一例であるAlphaにおける仮数部の加算処理について説明する。
【0229】
図12は、Alphaの加算処理の内容を表わす図である。同図は、C言語からadd(c,a,b,n)で呼び出されるアセンブラルーチンを表わす図である。
【0230】
add(c,a,b,n)は、メインメモリ6内の要素数nの配列aと、メインメモリ6内の要素数nの配列bとを加算して、加算結果をメインメモリ6内の要素数nの配列cに蓄積するルーチンである。
【0231】
add(c,a,b,n)が呼び出されると、制御部53は、カウンタiの値を保持するためのレジスタ$19にnの値を設定し、レジスタ$17に配列aへのポインタを設定し、レジスタ$18に配列bへのポインタを設定し、レジスタ$16に配列cへのポインタを設定する。
【0232】
2行目では、制御部53は、配列a[n]のアドレスを計算して、それをレジスタ$17に格納する。
【0233】
3行目では、制御部53は、配列b[n]のアドレスを計算して、それをレジスタ$18に格納する。
【0234】
4行目では、制御部53は、配列c[n]のアドレスを計算して、それをレジスタ$16に格納する。
【0235】
6行目では、制御部53は、繰上りの有無を表わす情報を保持するためのレジスタ$0の値を「0」(つまり、繰上りなし)に設定する。
【0236】
9行目では、制御部53は、メインメモリ6内のa[i−1]の値をレジスタ$1にロードする。
【0237】
10行目では、制御部53は、メインメモリ6内のb[i−1]の値をレジスタ$2にロードする。
【0238】
12行目では、演算器54は、加算を実行する。すなわち、演算器54は、レジスタ$1内のデータとレジスタ$2内のデータとを加算し、加算結果をレジスタ$5に保存する。
【0239】
13行目では、演算器54は、12行目の加算により繰上りがあったか否かを調べる。すなわち、制御部53は、レジスタ$5内のデータ(=加算結果c)がレジスタ$2内のデータ(=b)よりも小さい場合には、繰上りがあったことがわかるので、レジスタ$23に「1」を格納し、それ以外の場合には、繰上りがないことがわかるので、レジスタ$23に「0」を格納する。
【0240】
15行目では、演算器54は、レジスタ$5内のデータ(=加算結果)と、レジスタ$0内に格納されている1回前のループでの繰上りの有無を表わすデータとを加算して、加算結果をレジスタ$6に格納する。
【0241】
16行目では、演算器54は、15行目の加算により繰上りがあったか否かを調べる。すなわち、制御部53は、レジスタ$6内のデータがレジスタ$0内のデータよりも小さい場合には、繰上りがあったことがわかるので、レジスタ$0に「1」を格納し、それ以外の場合には、繰上りがないことがわかるので、レジスタ$0に「0」を格納する。
【0242】
18行目では、制御部53は、レジスタ$6内の加算結果を、メインメモリ6内のc[i−1]の領域に保存する。
【0243】
19行目では、演算器54は、レジスタ$0内のデータと、レジスタ$23内のデータとの論理和をレジスタ$0に格納することにより、13行目の繰上り判定の結果と16行目の繰上り判定の結果とを統合する。このようにして判定結果が統合できるのは、レジスタ$0内のデータとレジスタ$23内のデータの両方が「1」となる場合は、ありえないからである。
【0244】
21行目では、制御部53は、配列a[i−2]のアドレスを計算して、それをレジスタ$17に格納する。
【0245】
22行目では、制御部53は、配列b[i−2]のアドレスを計算して、それをレジスタ$18に格納する。
【0246】
23行目では、制御部53は、配列c[i−2]のアドレスを計算して、それをレジスタ$16に格納する。
【0247】
25行目では、演算器54は、カウンタiを「1」だけ減らす。
26行目では、制御部53は、カウンタiが「0」の場合には、ループを抜け、カウンタiが「0」以外の場合には、ループの先頭に戻るように制御する。
【0248】
(加算処理の比較)
図11と図12を比較すると、図11の8行目の繰上り付き加算命令adcqと同等の機能を実現するために、図12では、12行目の加算命令addqと、13行目の比較命令cmplt、15行目の加算命令addqと、16行目の比較命令cmpltと、19行目の論理和命令orの5つの命令が必要となる。以上のことから、本発明の実施形態のようにキャリーフラグCFを利用することにより、それを利用しない場合と比べて、仮数部の加算処理の命令数を削減することができることがわかる。
【0249】
(減算処理)
本発明の実施形態における仮数部の減算処理について説明する。
【0250】
図13は、フラグレジスタ57を備えるCPU3の一例であるAMD64の減算処理の内容を表わす図である。同図は、C言語からsub(c,a,b,n)で呼び出されるアセンブラルーチンを表わす図である。
【0251】
sub(c,a,b,n)は、メインメモリ6内の要素数nの配列aから、メインメモリ6内の要素数nの配列bを減算して、減算結果をメインメモリ6内の要素数nの配列cに蓄積するルーチンである。
【0252】
sub(c,a,b,n)が呼び出されると、制御部53は、カウンタiの値を保持するためのレジスタ%rcxにnの値を設定し、レジスタ%rsiに配列aへのポインタを設定し、レジスタ%rdxに配列bへのポインタを設定し、レジスタ%rdiに配列cへのポインタを設定する。
【0253】
2行目では、制御部53は、最上位ビットでの繰下がりの有無を表わす情報を保持するためのレジスタ%raxの値を「0」(つまり、繰下がりなし)に設定するとともに、フラグレジスタ57内のキャリーフラグCFの値を「0」(つまり、繰下がりなし)に初期化する。
【0254】
6行目では、制御部53は、メインメモリ6内のa[i−1]の値をレジスタ%r8にロードする。
【0255】
7行目では、制御部53は、メインメモリ6内のb[i−1]の値をレジスタ%r10にロードする。
【0256】
8行目では、演算器54は、繰下がり付き減算を実行する。すなわち、演算器54は、レジスタ%r8内のデータから、レジスタ%r10内のデータおよびフラグレジスタ57内のキャリーフラグCFの値とを減算し、減算結果をレジスタ%r10に保存するとともに、減算により繰下がりがあった場合に、フラグレジスタ57内のキャリーフラグCFの値を「1」に設定し、繰下がりがなかった場合に、フラグレジスタ57内のキャリーフラグCFの値を「0」に設定する。
【0257】
9行目では、制御部53は、レジスタ%r10内の加算結果を、メインメモリ6内のc[i−1]の領域に保存する。
【0258】
10行目では、演算器54は、カウンタiを「1」だけ減らす。
11行目では、制御部53は、カウンタiが「0」の場合には、ループを抜け、カウンタiが「0」以外の場合には、ループの先頭に戻るように制御する。
【0259】
13行目では、演算器54は、即値「0」と、レジスタ%rax内のデータ(=「0」)と、フラグレジスタ57内のキャリーフラグCFの値とを加算し、加算結果をレジスタ%raxに保存する。したがって、レジスタ%raxには、最上位ビットでの繰り下がりの有無を表わす情報が保持される。
【0260】
(参考:キャリーフラグを利用しない減算処理)
次に、比較のために、フラグレジスタ57を備えないCPU3の一例であるAlphaにおける仮数部の減算処理について説明する。
【0261】
図14は、Alphaの減算処理の内容を表わす図である。同図は、C言語からsub(c,a,b,n)で呼び出されるアセンブラルーチンを表わす図である。
【0262】
sub(c,a,b,n)は、メインメモリ6内の要素数nの配列aから、メインメモリ6内の要素数nの配列bと減算して、減算結果をメインメモリ6内の要素数nの配列cに蓄積するルーチンである。
【0263】
sub(c,a,b,n)が呼び出されると、制御部53は、カウンタiの値を保持するためのレジスタ$19にnの値を設定し、レジスタ$17に配列aへのポインタを設定し、レジスタ$18に配列bへのポインタを設定し、レジスタ$16に配列cへのポインタを設定する。
【0264】
2行目では、制御部53は、配列a[n]のアドレスを計算して、それをレジスタ$17に格納する。
【0265】
3行目では、制御部53は、配列b[n]のアドレスを計算して、それをレジスタ$18に格納する。
【0266】
4行目では、制御部53は、配列c[n]のアドレスを計算して、それをレジスタ$16に格納する。
【0267】
6行目では、制御部53は、繰下がりの有無を表わす情報を保持するためのレジスタ$0の値を「0」(つまり、繰下がりなし)に設定する。
【0268】
9行目では、制御部53は、メインメモリ6内のa[i−1]の値をレジスタ$1にロードする。
【0269】
10行目では、制御部53は、メインメモリ6内のb[i−1]の値をレジスタ$2にロードする。
【0270】
12行目では、演算器54は、減算を実行する。すなわち、演算器54は、レジスタ$1内のデータからレジスタ$2内のデータを減算し、減算結果をレジスタ$5に保存する。
【0271】
13行目では、演算器54は、12行目の減算により繰下がりがあったか否かを調べる。すなわち、制御部53は、レジスタ$1内のデータ(=a)がレジスタ$2内のデータ(=b)よりも小さい場合には、繰下がりがあったことがわかるので、レジスタ$6に「1」を格納し、それ以外の場合には、繰下がりがないことがわかるので、レジスタ$6に「0」を格納する。
【0272】
15行目では、演算器54は、レジスタ$5内のデータ(=減算結果)から、レジスタ$0内に格納されている1回前のループでの繰下がりの有無を表わすデータを減算して、減算結果をレジスタ$7に格納する。
【0273】
16行目では、演算器54は、15行目の減算により繰下がりがあったか否かを調べる。すなわち、制御部53は、レジスタ$5内のデータがレジスタ$0内のデータよりも小さい場合には、繰下がりがあったことがわかるので、レジスタ$0に「1」を格納し、それ以外の場合には、繰下がりがないことがわかるので、レジスタ$0に「0」を格納する。
【0274】
18行目では、制御部53は、レジスタ$7内の減算結果を、メインメモリ6内のc[i−1]の領域に保存する。
【0275】
19行目では、演算器54は、レジスタ$0内のデータと、レジスタ$6内のデータとの論理和をレジスタ$0に格納することにより、13行目の繰下がり判定の結果と16行目の繰下がり判定の結果とを統合する。このようにして判定結果が統合できるのは、レジスタ$0内のデータとレジスタ$6内のデータの両方が「1」となる場合は、ありえないからである。
【0276】
21行目では、制御部53は、配列a[i−2]のアドレスを計算して、それをレジスタ$17に格納する。
【0277】
22行目では、制御部53は、配列b[i−2]のアドレスを計算して、それをレジスタ$18に格納する。
【0278】
23行目では、制御部53は、配列c[i−2]のアドレスを計算して、それをレジスタ$16に格納する。
【0279】
25行目では、演算器54は、カウンタiを「1」だけ減らす。
26行目では、制御部53は、カウンタiが「0」の場合には、ループを抜け、カウンタiが「0」以外の場合には、ループの先頭に戻るように制御する。
【0280】
(減算処理の比較)
図13と図14を比較すると、図13の8行目の繰下がり付き減算命令sbbqと同等の機能を実現するために、図14では、12行目の減算命令subqと、13行目の比較命令cmplt、15行目の減算命令subqと、16行目の比較命令cmpltと、19行目の論理和命令orの5つの命令が必要となる。以上のことから、本発明の実施形態のようにキャリーフラグCFを利用することにより、それを利用しない場合と比べて、加算部の減算処理の命令数を削減することができることがわかる。
【0281】
(乗算処理)
まず、本発明の実施形態で用いる乗算の基本的な処理内容を説明する。図15に示すように、要素数4の配列Aと、要素数4の配列Bとの乗算を考える。配列A[0]のデータがa0、配列A[1]のデータがa1、配列A[2]のデータがa2、配列A[3]のデータがa3とする。配列B[0]のデータがb0、配列B[1]のデータがb1、配列B[2]のデータがb2、配列B[3]のデータがb3とする。a0、a1、a2、a3、b0、b1、b2、b3のデータ長は、それぞれ64ビットである。図15に示すように、配列Aと、b0、b1、b2、b3それぞれとの乗算結果は、64×5ビットとなる。配列Aと配列Bとを乗算すると、乗算結果は64×8ビットとなる。
【0282】
浮動小数点表現の仮数部の積を求めることを考えた場合、このような64×8ビットがすべて必要になるわけではない。下位のビットは、実際上無視しても、計算精度に大きな影響はない。
【0283】
本発明の実施形態では、要素数nの配列A(つまり64×nビット)と、要素数nの配列B(つまり64×nビット)との乗算では、上位半分の64×nビットのみを求めて、要素数nの配列Cに格納する。
【0284】
図15で示される配列Aと配列Bの乗算の場合には、上位から64×4ビットだけを求める。この場合、A×b3の計算では、a0×b3の計算は必要だが、a1×b3、a2×b3、a3×b3の計算は省略することができる。一方、A×b0の計算では、a0×b0、a1×b0、a2×b0、a3×b0のすべての計算が必要となる。
【0285】
図16は、本発明の実施形態における乗算処理の内容を表わす図である。
図16を参照して、この乗算処理は、要素数nの配列A(各要素は64ビット)と、要素数nの配列B(各要素は64ビット)とを乗算して、乗算結果のうち上位半分の64×nビットを、要素数nの配列Cに格納する。
【0286】
8行目のmul_addは、配列Aの先頭から(n−i)個の要素(A[0]、・・・、A[n−i−1])と、配列B[i]とを乗算して、乗算結果を配列C′[i]から後ろに加えるルーチンである。たとえば、図17に示すように、nが4で、iが3の場合には、B[3](=b3)とA[0](=a0)の乗算だけが行なわれて、乗算結果が配列C′[3]から後ろに加えられる。
【0287】
(本発明の実施形態のmul_add)
図18は、フラグレジスタ57を備えるCPU3の一例であるAMD64のmul_addの内容を表わす図である。同図は、C言語から、mul_add(c,a,b,n)で呼び出されるアセンブラルーチンを表わす図である。
【0288】
mul_add(c,a,b,n)は、メインメモリ6内の要素数nの配列aと、メインメモリ6内の変数bとを乗算して、乗算結果をメインメモリ6内の要素数n+1の配列cに加算する命令である。
【0289】
mul_add(c,a,b,n)が呼び出されると、制御部53は、レジスタ%rdiに配列cへのポインタを設定し、レジスタ%rsiに配列aへのポインタを設定し、レジスタ%rdxに乗数bの値を設定し、カウンタiの値を保持するためのレジスタ%rcxにnの値を設定する。
【0290】
1行目では、制御部53は、繰上りの有無を表わす情報を保持するためのレジスタ%r9の値を「0」(つまり、繰上りなし)に設定するとともに、フラグレジスタ57内のキャリーフラグCFの値を「0」(つまり、繰上りなし)に初期化する。
【0291】
2行目では、制御部53は、レジスタ%rdx内の乗数bの値をレジスタ%r8に格納する。
【0292】
5行目では、制御部53は、メインメモリ6内のa[i−1]の値をレジスタ%raxにロードする。
【0293】
6行目では、制御部53は、レジスタ%rax内のデータ(=a[i−1])と、レジスタ%r8内のデータ(=b)とを乗算する。乗算結果の上位ビットはレジスタ%rdxに格納され、乗算結果の下位ビットはレジスタ%raxに格納される。
【0294】
9行目では、演算器54は、繰上りなし加算を実行する。すなわち、演算器54は、レジスタ%r9に格納されている前回のループによる繰上りの有無を表わすデータと、レジスタ%rdxに格納されている乗算結果の上位ビットとを加算して、加算結果をレジスタ%rdxに格納する。この加算では、繰上りが生じないので、フラグレジスタ57内のキャリーフラグCFは「0」となる。
【0295】
10行目では、演算器54は、繰上りなし加算を実行する。すなわち、演算器54は、レジスタ%raxに格納されている乗算結果の下位ビットと、メインメモリ6内の配列c[i]とを加算して、加算結果をメインメモリ6内の配列c[i]に格納する。この加算によって、繰上りが生じた場合は、キャリーフラグCFは「1」となり、繰上りが生じなかった場合は、キャリーフラグCFは「0」となる。
【0296】
11行目では、制御部53は、繰上りの有無を表わす情報を保持するためのレジスタ%r9の値を「0」(つまり、繰上りなし)に設定する。
【0297】
12行目では、演算器54は、繰上り付き加算を実行する。すなわち、演算器54は、フラグレジスタ57内のキャリーフラグCFの値と、レジスタ%rdxに格納されている乗算結果の上位ビットと、メインメモリ6内の配列c[i−1]とを加算して、加算結果をメインメモリ6内の配列c[i−1]に格納するとともに、加算により繰上りがあった場合に、フラグレジスタ57内のキャリーフラグCFの値を「1」に設定し、繰上りがなかった場合に、フラグレジスタ57内のキャリーフラグCFの値を「0」に設定する。
【0298】
13行目では、演算器54は、繰上り付き加算を実行する。すなわち、演算器54は、即値「0」と、レジスタ%r9内の繰上りの有無を表わす情報と、フラグレジスタ57内のキャリーフラグCFの値とを加算して、加算結果をレジスタ%r9に格納する。この加算では、繰上りが生じないので、フラグレジスタ57内のキャリーフラグCFは「0」となる。
【0299】
15行目では、演算器54は、カウンタiを「1」だけ減らす。
16行目では、制御部53は、カウンタiが「0」の場合には、ループを抜け、カウンタiが「0」以外の場合には、ループの先頭に戻るように制御する。
【0300】
(参考:キャリーフラグを利用しないmul_add)
次に、比較のために、フラグレジスタ57を備えないCPU3の一例であるAlphaにおけるmul_addについて説明する。
【0301】
図19は、Alphaのmul_addの内容を表わす図である。同図は、C言語から、mul_add(c,a,b,n)で呼び出されるアセンブラルーチンを表わす図である。
【0302】
mul_add(c,a,b,n)は、メインメモリ6内の要素数nの配列aと、メインメモリ6内の変数bとを乗算して、乗算結果をメインメモリ6内の要素数n+1の配列cに加算するルーチンである。
【0303】
mul_add(c,a,b,n)が呼び出されると、制御部53は、カウンタiの値を保持するためのレジスタ$19にnの値を設定し、レジスタ$17に配列aへのポインタを設定し、レジスタ$18に乗数bの値を設定し、レジスタ$16に配列cへのポインタを設定する。
【0304】
1行目では、制御部53は、繰上がりの有無を表わす情報を保持するためのレジスタ$23の値を「0」(つまり、繰上がりなし)に設定する。
【0305】
2行目では、制御部53は、配列c[n]のアドレスを計算して、それをレジスタ$16に格納する。
【0306】
3行目では、制御部53は、配列a[n]のアドレスを計算して、それをレジスタ$17に格納する。
【0307】
6行目では、制御部53は、メインメモリ6内のa[i−1]の値をレジスタ$1にロードする。
【0308】
7行目では、制御部53は、メインメモリ6内のc[i]をレジスタ$21にロードする。
【0309】
8行目では、制御部53は、メインメモリ6内のc[i−1]をレジスタ$20にロードする。
【0310】
10行目では、演算器54は、レジスタ$1内のa[i−1]とレジスタ$18内のbとを乗算して、乗算結果の下位64ビットの値をレジスタ$2に格納する。
【0311】
11行目では、演算器54は、レジスタ$1内のa[i−1]とレジスタ$18内のbとを乗算して、乗算結果の上位64ビットの値をレジスタ$1に格納する。
【0312】
13行目では、演算器54は、レジスタ$1内の乗算結果の上位64ビットと、レジスタ$23に格納されている1回前のループでの繰上りの有無を表わすデータとを加算して、加算結果をレジスタ$1に格納する。
【0313】
15行目では、演算器54は、レジスタ$2内の乗算結果の下位64ビットと、レジスタ$21内のc[i]とを加算して、レジスタ$21に格納する。
【0314】
16行目では、制御部53は、15行目の加算により繰上がりがあったか否かを調べる。すなわち、制御部53は、レジスタ$21内のデータ(=c[i])がレジスタ$2内のデータよりも小さい場合には、繰上がりがあったことがわかるので、レジスタ$2に「1」を格納し、それ以外の場合には、繰上がりがないことがわかるので、レジスタ$2に「0」を格納する。
【0315】
18行目では、演算器54は、レジスタ$1内の乗算結果の上位64ビット(13行目で繰上り処理されたもの)と、レジスタ$20内のc[i−1]とを加算して、レジスタ$20に格納する。
【0316】
19行目では、制御部53は、18行目の加算により繰上がりがあったか否かを調べる。すなわち、制御部53は、レジスタ$20内のデータ(=c[i−1])がレジスタ$1内のデータよりも小さい場合には、繰上がりがあったことがわかるので、レジスタ$1に「1」を格納し、それ以外の場合には、繰上がりがないことがわかるので、レジスタ$1に「0」を格納する。
【0317】
21行目では、演算器54は、レジスタ$2内の15行目の加算の繰上りの有無を表わすデータと、レジスタ$20内のc[i−1](上位64ビットが加算後のもの)とを加算して、レジスタ$20に格納する。
【0318】
22行目では、制御部53は、21行目の加算により繰上がりがあったか否かを調べる。すなわち、制御部53は、レジスタ$20内のデータ(=c[i−1])がレジスタ$2内のデータよりも小さい場合には、繰上がりがあったことがわかるので、レジスタ$2に「1」を格納し、それ以外の場合には、繰上がりがないことがわかるので、レジスタ$2に「0」を格納する。
【0319】
24行目では、制御部53は、レジスタ$21内のc[i]をメインメモリ6内の領域に保存する。
【0320】
25行目では、制御部53は、レジスタ$20内のc[i−1]をメインメモリ6内の領域に保存する。
【0321】
26行目では、演算器54は、レジスタ$1内のデータと、レジスタ$2内のデータとの論理和をレジスタ$23に格納することにより、19行目の繰上がり判定の結果と22行目の繰上がり判定の結果とを統合する。このようにして判定結果が統合できるのは、レジスタ$1内のデータとレジスタ$2内のデータの両方が「1」となる場合は、ありえないからである。
【0322】
28行目では、制御部53は、配列a[i−2]のアドレスを計算して、それをレジスタ$17に格納する。
【0323】
29行目では、制御部53は、配列c[i−2]のアドレスを計算して、それをレジスタ$16に格納する。
【0324】
31行目では、演算器54は、カウンタiを「1」だけ減らす。
32行目では、制御部53は、カウンタiが「0」の場合には、ループを抜け、カウンタiが「0」以外の場合には、ループの先頭に戻るように制御する。
【0325】
(mul_addおよび乗算処理の比較)
図18と図19とを比較すると、図18では、ループ内部の命令数が7であるのに対して、図19では、ループ内部の命令数が18である。以上のことから、本発明の実施形態のようにキャリーフラグCFを利用することにより、それを利用しない場合と比べて、配列と配列の1要素との乗算であるmul_addの命令数を削減することができることがわかる。その結果、キャリーフラグCFを利用することにより、仮数部を表わす配列と配列の乗算の命令数も削減できる。
【0326】
(除算処理)
本発明の実施形態では、古典的なニュートン法により除算を行なう。x=a/bを求めるために、まず、1/bを求め、その結果にaを乗ずる。1/bを求めるために、f(x)=(1/x)−bについてニュートン法を適用する。このとき、ニュートン法の漸化式は、たとえば、以下の式(45)、(46)で与えられる。
【0327】
x0=1 ・・・(45)
xn+1=xn−f(xn)/f′(xn)=xn×(2−b×xn) ・・・(46)
図20は、ニュートン法を用いた除算処理の内容を説明するための図である。図21は、図20中のサブルーチンcmpの内容を表わす図である。
【0328】
図20において、21行目は、q=b×xnの乗算を実行するルーチンである。
22行目は、q=(2−q)の減算を実行するルーチンである。
【0329】
23行目は、q=xn×qの乗算を実行するルーチンである。
25行目および26行目において、xn+1(=q)がxnと等しければ、20行目〜30行目のループを抜ける。
【0330】
32行目は、xn、すなわち(1/b)にaを乗ずるルーチンである。
(除算処理の比較)
図20の21行目、23行目、32行目の乗算ルーチンmulを実行するCPU3がAMD64の場合では、図16および図18に示されるものであるのに対して、CPU3がAlphaの場合では、図16および図19に示されるものである。両者を比較すると、CPU3がAMD64の方が、命令数が少ないことがわかる。
【0331】
また、図20の22行目の減算ルーチンsubを実行するCPU3がAMD64の場合では、図13に示されるものであるのに対して、CPU3がAlphaの場合では、図14に示されるものである。両者を比較すると、CPU3がAMD64の方が、命令数が少ないことがわかる。以上のことから、本発明の実施形態のようにキャリーフラグCFを利用することにより、それを利用しない場合と比べて、仮数部の除算処理の命令数を削減することができることがわかる。
【0332】
(計算の精度についての実験結果)
本発明の実施の形態では、仮数部が64×nビットで表わされるが、nを大きくすることで、fa,s(n)(t)の計算精度に影響する正則化パラメータaの値を小さくすることができる。以下では、正則化パラメータaの大きさによって、fa,s(n)(t)がどのように変化するかの実験結果を説明する。f(t)がステップ函数の場合に、ラプラス変換像F(pi)が解析的に求まるが、これを用いて、実験を行なった。
【0333】
図22(a)は、正則化パラメータa=10-1、10-4、10-8、10-12のときの、計算結果fa,s(n)(t)を示す図である。図22(b)は、正則化パラメータa=10-100、10-400のときの、計算結果fa,s(n)(t)を示す図である。
【0334】
これらの図から、正則化パラメータaを小さくするほど、fa,s(n)(t)の計算精度が増加することがわかる。図22内のステップ函数の場合には、正則化パラメータaは、少なくとも10-100ぐらい小さくなくては有効な逆ラプラス変換像が得られないことがわかる。正則化パラメータaの値をこのような小さな値にして、かつ現実的な時間内で逆ラプラス変換像を算出するためには、本発明の実施形態で説明したように、キャリーフラグCFを用いて命令数を大幅に削減した多倍長演算が必要不可欠となる。
【0335】
また、本発明の実施の形態では、式(6)を修正した式(10)を用いることによって、f(0)が0以外の函数についても、逆ラプラス変換の値を高精度で求めることができるという特徴がある。以下では、軟化子パラメータsの大きさによって、fa,s(n)(t)がどのように変化するかの実験結果を説明する。ここでは、式(47)、(48)で表される函数f(t)を実験に用いた。
【0336】
0<t<1において、f(t)=1 ・・・(47)
1<tにおいて、f(t)=0 (48)
図23(a)は、s=0.1での計算結果fa,s(n)(t)を示す図であり、図23(b)は、s=0.01での計算結果fa,s(n)(t)を示す図である。
【0337】
これらの図から、軟化子パラメータsを小さくするほど、fa,s(n)(t)の計算精度が増加することがわかる。軟化子パラメータsを目的に応じて適当な値に設定することによって、f(0)が0以外の函数について、逆ラプラス変換の値を高精度で求めることができる。
【0338】
[第4の実施形態]
第4の実施形態は、Hテーブルを作成する処理と、Hテーブルを用いて逆ラプラス変換の数値解を算出する処理とを別個の装置で行なう構成に関するものである。
【0339】
図24は、Hテーブルを作成する逆ラプラス変換用テーブル作成装置の構成を表わす図である。図24の逆ラプラス変換用テーブル作成装置は、コンピュータで実現されるものであり、図24において、図1と同一の構成要素については同一の符号を付している。
【0340】
図24を参照して、CPU63は、第1〜第3の実施形態と異なり、テーブル作成部4としてのみ機能する。
【0341】
入力部62は、第1〜第3の実施形態と異なり、Hテーブルの作成に必要なサンプル数n、上限値U、下限値L、正則化パラメータaが入力される。
【0342】
プログラム記憶部65には、第1〜第3の実施形態と異なり、コンピュータをテーブル作成部4として機能させるための逆ラプラス変換のためのテーブル作成プログラムが記憶されている。逆ラプラス変換のためのテーブル作成プログラムは、USB(Universal Serial Bus)メモリ、CD−ROM(Compact Disk - Read Only Memory)などのコンピュータ読取り可能な記録媒体であるリムーマブルメディア67を通じて、外部からインストールすることができる。
【0343】
変数記憶部66には、第1〜第3の実施形態と異なり、Hテーブル作成にのみ必要な変数が記憶される。
【0344】
メインメモリ64内のHテーブル記憶部12に記憶されているHテーブルと、変数記憶部66に記憶されている変数h、pi、qi(i=0,1,2,・・・,n)は、リムーマブルメディア67に出力することができる。これにより、リムーマブルメディア67に記憶されたHテーブルおよび変数を他の装置のメインメモリに送り記憶させることで、他の装置で、Hテーブルおよび変数h、pi、qiを利用することができる。
【0345】
図25は、逆ラプラス変換の数値解を算出する逆ラプラス変換の数値解算出装置の構成を表わす図である。図25の逆ラプラス変換の数値解算出装置は、コンピュータで実現されるものであり、図25において、図1と同一の構成要素については同一の符号を付している。
【0346】
図25を参照して、CPU73は、第1〜第3の実施形態と個なり、逆変換部5としてのみ機能する。
【0347】
入力部72は、第1〜第3の実施形態と異なり、逆ラプラス変換の数値解算出のために必要なラプラス変換像F(pi)と、軟化子パラメータsとが入力される。
【0348】
プログラム記憶部75には、第1〜第3の実施形態と異なり、コンピュータを逆変換部5および出力部13として機能させるための逆ラプラス変換の数値解算出プログラムが記憶されている。逆ラプラス変換の数値解算出プログラムは、USB(Universal Serial Bus)メモリ、CD−ROM(Compact Disk - Read Only Memory)などのコンピュータ読取り可能な記録媒体であるリムーマブルメディア77を通じて、外部からインストールすることができる。
【0349】
リムーマブルメディア77に記憶されているHテーブルと、変数h、pi、qi(i=0,1,2,・・・,n)は、それぞれメインメモリ74内のHテーブル記憶部12と、変数記憶部76に出力することができる。これにより、図24の装置で作成されてリムーマブルメディア67に記憶されたHテーブルおよび変数h、pi、qi(i=0,1,2,・・・,n)を、図25の装置で利用することができる。
【0350】
以上のように、第4の実施形態では、Hテーブルの作成と、逆ラプラス変換の数値解の算出とを別個の装置で行なうことができるので、Hテーブルの作成を高速の計算機で行なってユーザに提供し、ユーザは、自身のパソコン上で提供を受けたHテーブルを用いて逆ラプラス変換の数値解を算出を行なうようにすることができる。
【0351】
(変形例)
本発明は、上記の実施の形態に限定されるものではなく、たとえば以下のような変形例を含む。
【0352】
(1) テーブルに記述する情報
本発明の第4の実施形態では、図24の逆ラプラス変換用テーブル作成装置が、連立方程式の解{Hi n,a,t:i=0,1,2,・・・,n}を記述したHテーブルと、h、pi、qi(i=0,1,2,・・・,n)とをリムーバブルメディアに出力して、図25の逆ラプラス変換の数値解算出装置でこれを利用したが、これに限定するものではない。たとえば、逆ラプラス変換用テーブル作成装置が、Hテーブルだけをリムーバブルメディアに出力し、逆ラプラス変換の数値解算出装置では、n、U、Lを入力して、h、pi、qjを再計算するものとしてもよい。また、逆ラプラス変換用テーブル作成装置が、h×pi×Hi n,a,t×qiを計算して、計算結果を記述したテーブルとpiとをリムーバブルメディアに出力し、逆ラプラス変換の数値解算出装置で、これらを利用して、式(38)で表わされる計算をより簡単に計算するようにしてもよい。
【0353】
(2) 離散化の手法
本発明の実施形態で説明した式(22)〜(24)を離散化する手法は、様々な方法があり、式(25)〜(36)はその一例にすぎず、その他の方法で行なってもよい。
【0354】
(3) 64ビットレジスタ
本発明の第3の実施形態では、CPUが64ビット長のレジスタを備える場合を例にして説明したが、これは一例であって、Kビット長(Kは自然数)のレジスタ、たとえば16ビット長、32ビット長、128ビット長などのレジスタであってもよい。
【0355】
(4) LU分解
本発明の第1の実施形態では、第二種積分方程式を離散化して得られる連立方程式を構成する係数行列AをLU分解したが、これに限定するものではなく、下三角行列と上三角行列の積に分解できるものなら、どのような分解方法であってもよい。
【0356】
今回開示された実施の形態はすべての点で例示であって制限的なものではないと考えられるべきである。本発明の範囲は上記した説明ではなくて特許請求の範囲によって示され、特許請求の範囲と均等の意味および範囲内でのすべての変更が含まれることが意図される。
【図面の簡単な説明】
【0357】
【図1】第1の実施形態の逆ラプラス変換装置の構成を表わす図である。
【図2】Hテーブルを表わす図である。
【図3】第1の実施形態におけるHテーブルを作成する手順を表わすフローチャートである。
【図4】第1の実施形態におけるHテーブルを用いて逆ラプラス変換の数値解を算出する手順を表わすフローチャートである。
【図5】第1の実施形態の変形例の逆ラプラス変換装置の構成を表わす図である。
【図6】第1の実施形態の変形例におけるHテーブルを作成する手順を表わすフローチャートである。
【図7】第2の実施形態におけるHテーブルを作成する手順を表わすフローチャートである。
【図8】第3の実施形態のCPUの具体的な構成を表わす図である。
【図9】フラグレジスタの具体的な構成を表わす図である。
【図10】(a)、(b)、(c)は、図8の逆ラプラス変換装置において演算の対象となる数値の表現を説明するための図である。
【図11】フラグレジスタを備えるCPUの一例であるAMD64の加算処理の内容を表わす図である。
【図12】Alphaの加算処理の内容を表わす図である。
【図13】フラグレジスタを備えるCPUの一例であるAMD64の減算処理の内容を表わす図である。
【図14】Alphaの減算処理の内容を表わす図である。
【図15】本発明の実施形態で用いる乗算処理の内容を説明するための図である。
【図16】本発明の実施形態における乗算処理の内容を表わす図である。
【図17】本発明の実施形態で用いる乗算処理の内容を説明するための図である。
【図18】フラグレジスタを備えるCPUの一例であるAMD64のmul_addの内容を表わす図である。
【図19】Alphaのmul_addの内容を表わす図である。
【図20】ニュートン法を用いた除算処理の内容を説明するための図である。
【図21】図20中のサブルーチンcmpの内容を表わす図である。
【図22】(a)は、正則化パラメータa=10-1、10-4、10-8、10-12のときの、計算結果fa,s(n)(t)を示す図である。(b)は、正則化パラメータa=10-100、10-400のときの、計算結果fa,s(n)(t)を示す図である。
【図23】(a)は、s=0.1での計算結果fa,s(n)(t)を示す図であり、(b)は、s=0.01での計算結果fa,s(n)(t)を示す図である。
【図24】Hテーブルを作成する逆ラプラス変換用テーブル作成装置の構成を表わす図である。
【図25】逆ラプラス変換の数値解を算出する逆ラプラス変換の数値解算出装置の構成を表わす図である。
【符号の説明】
【0358】
1,81 逆ラプラス変換装置、61 逆ラプラス変換用テーブル作成装置、71 逆ラプラス変換数値解算出装置、2,62,72 入力部、3,63 CPU、4,84 テーブル作成部、5 逆変換部、6,64,74 メインメモリ、7,65,75 プログラム記憶部、8,66,76 変数記憶部、9 三角行列記憶部、12 Hテーブル記憶部、13 出力部、14 表示装置、67,77,99 リムーバブルメディア、53 制御部、54 演算器、56 汎用レジスタ群、57 フラグレジスタ、89 逆行列記憶部。
【技術分野】
【0001】
本発明は、逆ラプラス変換プログラム、逆ラプラス変換のためのテーブル作成プログラム、逆ラプラス変換の数値解算出プログラム、および逆ラプラス変換装置に関し、特に再生核と正則化法用いて逆ラプラス変換の数値解を算出する技術に関する。
【背景技術】
【0002】
逆ラプラス変換は、電気電子工学、油田の調査、画像処理をはじめとする様々な分野で広く応用されている。従来、コンピュータを用いて逆ラプラス変換の数値解を算出するためには、複素積分の数値計算による方法や、ラプラス変換表を用いた方法が用いられていた。しかしながら、このような従来の方法では、丸め誤差または離散化誤差などの計算誤差が生じ、数値的に不安定(ill-conditioned)であるという問題がある。
【0003】
このような数値的不安定性を解消する手法として、再生核と正則化法に基づいて逆ラプラス変換の数値解を求める方法が提案されている(たとえば、非特許文献1を参照)。この方法では、逆ラプラス変換の数値解が、再生核ヒルベルト空間での最良近似函数として求めることができる。
【非特許文献1】松浦他、"Numerical Real Inversion formulas of the Laplace Transform by Using a Fredholm Integral Equation of the Second Kind"、日本数学会、2007年度年会、応用数学分科会、講演アブストラクト、p.69-72
【発明の開示】
【発明が解決しようとする課題】
【0004】
しかしながら、上記の再生核と正則化法を用いる手法では、原像が原点で零、かつ原像の導函数が増大しない絶対連続函数からなる再生核空間を考えており、原像がこれら2つの条件を満たす場合に対してのみ適用可能である。
【0005】
それゆえに、本発明の目的は、原像の導函数が増大する場合、または原像が原点で零でない場合においても、再生核と正則化法を用いて、逆ラプラス変換の数値解を計算することができる逆ラプラス変換プログラム、逆ラプラス変換のためのテーブル作成プログラム、逆ラプラス変換の数値解算出プログラム、および逆ラプラス変換装置を提供することである。
【課題を解決するための手段】
【0006】
本発明のある局面に係る逆ラプラス変換プログラムは、与えられた正則化パラメータと軟化子パラメータのもとで、ラプラス変換像に対する原像の数値解を算出するために、コンピュータに、原点で零であり、かつ絶対連続な函数からなる重み付き再生核ヒルベルト空間上で、重み付き二乗可積分空間を観測空間としたチホノフ正則化法により導かれる第二種積分方程式を離散化して得られる連立方程式の解を求め、連立方程式の解に基づく第二種積分方程式の数値解を含む情報を記述したテーブルを作成するステップと、作成したテーブルを記憶部に保存するステップと、記憶部のテーブルを参照して、第二種積分方程式の数値解と軟化子函数を乗じたラプラス変換像との重み付き二乗可積分空間での内積を数値計算で求めることにより、原像の数値解を算出するステップと、原像の数値解を出力するステップとを実行させる。
【0007】
好ましくは、正則化パラメータをaとし、軟化子パラメータsとし、任意の関数gのラプラス変換像をLgとし、数値解を算出する原像に対するラプラス変換像をF(p)とし、軟化子函数をRs(p)としたときに、原点で零であり、かつ絶対連続な函数fからなる重み付き再生核ヒルベルト空間Hkのノルムは、式(B1)で表わされ、
【0008】
【数1】
【0009】
重み付き再生核ヒルベルト空間Hkでの再生核K(t,t′)が式(B2)で表わされ、
【0010】
【数2】
【0011】
重み付き二乗可積分空間は、L2(R+,u(p)dp)で表わされ、ρ(t)とu(p)の間には、式(B3)の条件が成立し、
【0012】
【数3】
【0013】
第二種積分方程式は、式(B4)で表わされ、
【0014】
【数4】
【0015】
内積は、式(B5)で表わされる。
【0016】
【数5】
【0017】
好ましくは、ρ(t)は、式(B6)で表わされ、u(p)は、式(B7)で表わされ、Rs(p)は、式(B8)で表わされる。
【0018】
【数6】
【0019】
好ましくは、コンピュータは、Kビット長(Kは自然数)の汎用レジスタと、直前の演算により繰上りまたは繰下がりが生じたか否かを示すフラグを格納するフラグレジスタと、フラグを参照して、直前の演算により繰上りまたは繰下がりが生じたか否かに応じた演算を行なう演算器とを備え、テーブルを作成するステップおよび原像の数値解を算出するステップは、数値計算の対象となる変数について、その仮数部をK×Nビット(Nは2以上の自然数)の多倍長で表わし、変数の仮数部をKビットごとに区切って、区切られた下位側から順次、汎用レジスタを用いて演算器に演算を行なわせるステップを含む。
【0020】
好ましくは、連立方程式は、Ax=bの形で表わされ、Aは各要素がtに依存しない係数行列であり、xは連立方程式の解を要素とするベクトルであり、bは各要素がtに依存するベクトルであり、テーブルを作成するステップは、行列Aの逆行列A-1を算出して、逆行列A-1の要素を記憶部に記憶するステップと、計算範囲内のすべてのtについて、記憶部に記憶されている逆行列A-1の要素の値に基づいて、連立方程式の解を算出するステップとを含む。
【0021】
好ましくは、連立方程式は、Ax=bの形で表わされ、Aは各要素がtに依存しない係数行列であり、xは連立方程式の解を要素とするベクトルであり、bは各要素がtに依存するベクトルであり、テーブルを作成するステップは、行列Aを上三角行列と下三角行列の積に分解して、上三角行列と下三角行列の要素を記憶部に記憶するステップと、計算範囲内のすべてのtについて、記憶部に記憶されている上三角行列と下三角行列の要素の値に基づいて、連立方程式の解を算出するステップとを含む。
【0022】
好ましくは、テーブルを作成するステップは、正則化パラメータaおよび軟化子パラメータsと、その他のパラメータn、U、Lとから、式(B9)〜(B14)に従って、それぞれ変数h、xj、pj、qj、aij、H(pj, t)を計算するステップと、
【0023】
【数7】
【0024】
連立方程式は、式(B15)で表わされ、式(B16)で表わされる行列Aをコレスキー分解し、式(B15)の解を算出するステップと、
【0025】
【数8】
【0026】
式(B17)に従って、式(B15)の解から第二種積分方程式の数値解{Hin,a,t:i=0,1,2,・・・,n}を算出し、第二種積分方程式の数値解を含む情報を記述したテーブルを作成するステップとを含み、
【0027】
【数9】
【0028】
原像の数値解を算出するステップは、テーブルを参照して、式(B18)に従って、ラプラス変換像F(pj)から、原像の数値解fa,s(n)(t)を算出するステップを含む。
【0029】
【数10】
【0030】
本発明のある局面に係るテーブル作成プログラムは、与えられた正則化パラメータと軟化子パラメータのもとで、ラプラス変換像に対する原像の数値解を算出するためのテーブルを作成するために、コンピュータに、原点で零であり、かつ絶対連続な函数からなる重み付き再生核ヒルベルト空間上で、重み付き二乗可積分空間を観測空間としたチホノフ正則化法により導かれる第二種積分方程式を離散化して得られる連立方程式の解を求め、連立方程式の解に基づく第二種積分方程式の数値解を含む情報を記述したテーブルを作成するステップと、作成したテーブルを記憶部に保存するステップとを実行させる。
【0031】
好ましくは、正則化パラメータをaとし、軟化子パラメータsとし、任意の関数gのラプラス変換像をLgとし、原点で零であり、かつ絶対連続な函数fからなる重み付き再生核ヒルベルト空間Hkのノルムは式(B19)で表わされ、
【0032】
【数11】
【0033】
重み付き再生核ヒルベルト空間Hkでの再生核K(t,t′)が式(B20)で表わされ、
【0034】
【数12】
【0035】
重み付き二乗可積分空間は、L2(R+,u(p)dp)で表わされ、ρ(t)とu(p)の間には、式(B21)の条件が成立し、
【0036】
【数13】
【0037】
第二種積分方程式は、式(B22)で表わされる。
【0038】
【数14】
【0039】
好ましくは、ρ(t)は、式(B23)で表わされ、u(p)は、式(B24)で表わされる。
【0040】
【数15】
【0041】
好ましくは、コンピュータは、Kビット長(Kは自然数)の汎用レジスタと、直前の演算により繰上りまたは繰下がりが生じたか否かを示すフラグを格納するフラグレジスタと、フラグを参照して、直前の演算により繰上りまたは繰下がりが生じたか否かに応じた演算を行なう演算器とを備え、テーブルを作成するステップは、数値計算の対象となる変数について、その仮数部をK×Nビット(Nは2以上の自然数)の多倍長で表わし、変数の仮数部をKビットごとに区切って、区切られた下位側から順次、汎用レジスタを用いて演算器に演算を行なわせるステップを含む。
【0042】
好ましくは、連立方程式は、Ax=bの形で表わされ、Aは各要素がtに依存しない係数行列であり、xは連立方程式の解を要素とするベクトルであり、bは各要素がtに依存するベクトルであり、テーブルを作成するステップは、行列Aの逆行列A-1を算出して、逆行列A-1の要素を記憶部に記憶するステップと、計算範囲内のすべてのtについて、記憶部に記憶されている逆行列A-1の要素の値に基づいて、連立方程式の解を算出するステップとを含む。
【0043】
好ましくは、連立方程式は、Ax=bの形で表わされ、Aは各要素がtに依存しない係数行列であり、xは連立方程式の解を要素とするベクトルであり、bは各要素がtに依存するベクトルであり、テーブルを作成するステップは、行列Aを上三角行列と下三角行列の積に分解して、上三角行列と下三角行列の要素を記憶部に記憶するステップと、計算範囲内のすべてのtについて、記憶部に記憶されている上三角行列と下三角行列の要素の値に基づいて、連立方程式の解を算出するステップとを含む。
【0044】
好ましくは、テーブルを作成するステップは、正則化パラメータaおよび軟化子パラメータsと、その他のパラメータn、U、Lとから、式(B25)〜(B30)に従って、それぞれ変数h、xj、pj、qj、aij、H(pj, t)を計算するステップと、
【0045】
【数16】
【0046】
連立方程式は、式(B31)で表わされ、式(B32)で表わされる行列Aをコレスキー分解し、式(B31)の解を算出するステップと、
【0047】
【数17】
【0048】
式(B33)に従って、式(B31)の解から第二種積分方程式の数値解{Hin,a,t:i=0,1,2,・・・,n}を算出し、第二種積分方程式の数値解を含む情報を記述したテーブルを作成するステップとを含む。
【0049】
【数18】
【0050】
本発明のある局面に係る逆ラプラス変換の数値解算出プログラムは、上記記載の逆ラプラス変換のためのテーブル作成プログラムによって作成されたテーブルを用いて、ラプラス変換像に対する原像の数値解を作成するために、コンピュータに、上記記載の逆ラプラス変換のためのテーブル作成プログラムで作成されたテーブルが入力され、入力されたテーブルを記憶部に保存するステップと、記憶部のテーブルを参照して、第二種積分方程式の数値解と軟化子函数を乗じたラプラス変換像との重み付き二乗可積分空間での内積を数値計算で求めることにより、原像の数値解を算出するステップと、原像の数値解を出力するステップとを実行させる。
【0051】
本発明のある局面に係る逆ラプラス変換の数値解算出プログラムは、上記記載の逆ラプラス変換のためのテーブル作成プログラムによって作成されたテーブルを用いて、ラプラス変換像に対する原像の数値解を作成するために、コンピュータに、上記記載の逆ラプラス変換のためのテーブル作成プログラムで作成されたテーブルが入力され、入力されたテーブルを記憶部に保存するステップと、記憶部のテーブルを参照して、第二種積分方程式の数値解と軟化子函数Rs(p)を乗じたラプラス変換像F(p)との重み付き二乗可積分空間での内積を数値計算で求めることにより、原像の数値解を算出するステップと、内積は、式(B34)で表わされ、
【0052】
【数19】
【0053】
原像の数値解を出力するステップとを実行させる。
本発明のある局面に係る逆ラプラス変換の数値解算出プログラムは、上記記載の逆ラプラス変換のためのテーブル作成プログラムによって作成されたテーブルを用いて、ラプラス変換像に対する原像の数値解を作成するために、コンピュータに、上記記載の逆ラプラス変換のためのテーブル作成プログラムで作成されたテーブルが入力され、入力されたテーブルを記憶部に保存するステップと、記憶部のテーブルを参照して、第二種積分方程式の数値解と軟化子函数Rs(p)を乗じたラプラス変換像F(p)との重み付き二乗可積分空間での内積を数値計算で求めることにより、原像の数値解を算出するステップと、軟化子函数Rs(p)は、式(B35)で表わされ、
【0054】
【数20】
【0055】
内積は、式(B36)で表わされ、
【0056】
【数21】
【0057】
原像の数値解を出力するステップとを実行させる。
本発明のある局面に係る逆ラプラス変換の数値解算出プログラムは、上記記載の逆ラプラス変換のためのテーブル作成プログラムによって作成されたテーブルを用いて、逆ラプラス変換の数値解を作成するために、Kビット長(Kは自然数)の汎用レジスタと、直前の演算により繰上りまたは繰下がりが生じたか否かを示すフラグを格納するフラグレジスタと、フラグを参照して、直前の演算により繰上りまたは繰下がりが生じたか否かに応じた演算を行なう演算器とを備えたコンピュータに、上記記載の逆ラプラス変換のためのテーブル作成プログラムで作成されたテーブルが入力され、入力されたテーブルを記憶部に保存するステップと、記憶部のテーブルを参照して、第二種積分方程式の数値解と軟化子函数Rs(p)を乗じたラプラス変換像F(p)との重み付き二乗可積分空間での内積を数値計算で求めることにより、原像の数値解を算出するステップと、軟化子函数Rs(p)は、式(B37)で表わされ、
【0058】
【数22】
【0059】
内積は、式(B38)で表わされ、
【0060】
【数23】
【0061】
原像の数値解を出力するステップとを実行させ、原像の数値解を算出するステップは、数値計算の対象となる変数について、その仮数部をK×Nビット(Nは2以上の自然数)の多倍長で表わし、変数の仮数部をKビットごとに区切って、区切られた下位側から順次、汎用レジスタを用いて演算器に演算を行なわせるステップを含む。
【0062】
本発明のある局面に係る逆ラプラス変換の数値解算出プログラムは、上記記載の逆ラプラス変換のためのテーブル作成プログラムによって作成されたテーブルを用いて、逆ラプラス変換の数値解を作成するために、コンピュータに、上記記載の逆ラプラス変換のためのテーブル作成プログラムで作成されたテーブルが入力され、入力されたテーブルを記憶部に保存するステップと、記憶部のテーブルを参照して、式(B39)に従って、第二種積分方程式の数値解とラプラス変換像F(pj)から原像の数値解fa,s(n)(t)を算出するステップと、
【0063】
【数24】
【0064】
原像の数値解を出力するステップとを実行させる。
本発明のある局面に係る逆ラプラス変換装置は、与えられた正則化パラメータと軟化子パラメータのもとで、ラプラス変換像に対する原像の数値解を算出する逆ラプラス変換装置であって、原点で零であり、かつ絶対連続な函数からなる重み付き再生核ヒルベルト空間上で、重み付き二乗可積分空間を観測空間としたチホノフ正則化法により導かれる第二種積分方程式を離散化して得られる連立方程式の解を求め、連立方程式の解に基づく第二種積分方程式の数値解を含む情報を記述したテーブルを作成するテーブル作成部と、作成したテーブルを記憶する記憶部と、記憶部のテーブルを参照して、第二種積分方程式の数値解と軟化子函数を乗じたラプラス変換像との重み付き二乗可積分空間での内積を数値計算で求めることにより、原像の数値解を算出する逆変換部と、原像の数値解を出力する出力部とを備える。
【0065】
好ましくは、正則化パラメータをaとし、軟化子パラメータsとし、任意の関数gのラプラス変換像をLgとし、数値解を算出する原像に対するラプラス変換像をF(p)とし、軟化子函数をRs(p)としたときに、原点で零であり、かつ絶対連続な函数fからなる重み付き再生核ヒルベルト空間Hkのノルムは式(B40)で表わされ、
【0066】
【数25】
【0067】
重み付き再生核ヒルベルト空間Hkでの再生核K(t,t′)が式(B41)で表わされ、
【0068】
【数26】
【0069】
重み付き二乗可積分空間は、L2(R+,u(p)dp)で表わされ、ρ(t)とu(p)の間には、式(B42)の条件が成立し、
【0070】
【数27】
【0071】
第二種積分方程式は、式(B43)で表わされ、
【0072】
【数28】
【0073】
内積は、式(B44)で表わされる。
【0074】
【数29】
【0075】
好ましくは、ρ(t)は、式(B45)で表わされ、u(p)は、式(B46)で表わされ、Rs(p)は、式(B47)で表わされる。
【0076】
【数30】
【0077】
テーブル作成部および逆変換部とは、共通して、Kビット長(Kは自然数)の汎用レジスタと、直前の演算により繰上りまたは繰下がりが生じたか否かを示すフラグを格納するフラグレジスタと、フラグを参照して、直前の演算により繰上りまたは繰下がりが生じたか否かに応じた演算を行なう演算器とを備え、演算器は、数値計算の対象となる変数について、その仮数部をK×Nビット(Nは2以上の自然数)の多倍長で表わし、変数の仮数部をKビットごとに区切って、区切られた下位側から順次、汎用レジスタを用いて演算を行なう。
【0078】
好ましくは、連立方程式は、Ax=bの形で表わされ、Aは各要素がtに依存しない係数行列であり、xは連立方程式の解を要素とするベクトルであり、bは各要素がtに依存するベクトルであり、テーブル作成部は、行列Aの逆行列A-1を算出して、逆行列A-1の要素を記憶部に記憶し、計算範囲内のすべてのtについて、記憶部に記憶されている逆行列A-1の要素の値に基づいて、連立方程式の解を算出する。
【0079】
好ましくは、連立方程式は、Ax=bの形で表わされ、Aは各要素がtに依存しない係数行列であり、xは連立方程式の解を要素とするベクトルであり、bは各要素がtに依存するベクトルであり、テーブル作成部は、行列Aを上三角行列と下三角行列の積に分解して、上三角行列と下三角行列の要素を記憶部に記憶し、計算範囲内のすべてのtについて、記憶部に記憶されている上三角行列と下三角行列の要素の値に基づいて、連立方程式の解を算出する。
【0080】
好ましくは、テーブル作成部は、正則化パラメータaおよび軟化子パラメータsと、その他のパラメータn、U、Lとから、式(B48)〜(B53)に従って、それぞれ変数h、xj、pj、qj、aij、H(pj, t)を計算し、
【0081】
【数31】
【0082】
連立方程式は、式(B54)で表わされ、テーブル作成部は、式(B55)で表わされる行列Aをコレスキー分解し、式(B54)の解を算出し、
【0083】
【数32】
【0084】
テーブル作成部は、さらに式(B56)に従って、式(B54)の解から第二種積分方程式の数値解{Hin,a,t:i=0,1,2,・・・,n}を算出し、第二種積分方程式の数値解を含む情報を記述したテーブルを作成し、
【0085】
【数33】
【0086】
逆変換部は、テーブルを参照して、式(B57)に従って、ラプラス変換像F(pj)から原像の数値解fa,s(n)(t)を算出する。
【0087】
【数34】
【発明の効果】
【0088】
本発明によれば、原像の導函数が増大する場合、または原像が原点で零でない場合においても、再生核と正則化法を用いて、逆ラプラス変換の数値解を計算することができる。
【発明を実施するための最良の形態】
【0089】
以下、本発明に係る実施の形態について図面を参照して説明する。
[第1の実施形態]
まず、第1の実施形態における逆ラプラス変換の基本的なアルゴリズムについて説明する。
【0090】
(再生核空間理論を用いた逆ラプラス変換)
ある自然な函数空間の函数f(t)のラプラス変換は、式(1)で表わされる。式(1)の逆変換を考えると、通常はBromwich積分による複素解析的な逆変換公式を考えるが、正の実軸上の値のみを用いた逆変換を求めたい場合がある。このような正の実軸上の値のみを用いた逆変換は実逆ラプラス変換と呼ばれる。また、f(t)は原像、F(P)をラプラス変換像と呼ばれる。
【0091】
【数35】
【0092】
式(2)で表わされる正の実軸R+上の有限ノルムを持ち、f(0)=0を満たす絶対連続な函数fからなる重み付き再生核ヒルベルト空間Hkを考える。この重み付き再生核ヒルベルト空間Hkは、原像の導函数f(t)′が増大する場合も包含している。
【0093】
【数36】
【0094】
この重み付き再生核ヒルベルト空間Hkは、式(3)で表わされる再生核K(t,t′)を有する。
【0095】
【数37】
【0096】
このとき、重み付き再生核ヒルベルト空間Hk上の重み付き二乗可積分空間L2(R+,u(p)dp)への線形写像(Lf)(p)pは有界であって、式(4)が成立する。
【0097】
【数38】
【0098】
式(4)から、一般論に従って、以下が成立する。
チホノフ正則化とは、任意のg∈L2(R+,u(p)dp)と任意の正則化パラメータa>0に対して、式(5)の右辺を最小化することである。重み付き再生核ヒルベルト空間Hk上で、重み付き二乗可積分空間L2(R+,u(p)dp)を観測空間としたチホノフ正則化によって、式(6)で表わされる最良近似函数fa(t)が導かれる。
【0099】
【数39】
【0100】
ここで、Ka(・,t)は、式(7)〜(9)で表わされる。
【0101】
【数40】
【0102】
上述の重み付き再生核ヒルベルト空間Hkは、f(0)=0を満たすことを条件としている。f(0)≠0なる函数f(t)についても、上述の重み付き再生核ヒルベルト空間Hkの特性が利用できるように、軟化子函数Rs(p)を用いて、式(6)の最良近似函数fa(t)を修正する。修正最良近似函数fa,s(t)は、式(10)で表わされる。ここで、sは、軟化子パラメータである。
【0103】
【数41】
【0104】
さて、式(9)のtについてラプラス変換をとり、t′を改めてtとおくと、式(11)が得られる。
【0105】
【数42】
【0106】
式(12)および(13)を用いると、式(11)から式(14)が成立する。式(14)は、書き下すと、式(15)のように表わされる。また、式(10)から、式(16)が成立する。
【0107】
【数43】
【0108】
式(14)は、未知の函数が積分の中に現れており、第二種積分方程式と呼ばれる。修正最良近似函数fa,ε(t)は、第二種積分方程式の解と、軟化子函数を乗じたラプラス変換像との重み付き二乗可積分空間での内積である式(16)で表わされる。本発明の実施形態の逆ラプラス変換では、修正最良近似函数fa,ε(t)を原像f(t)の近似解として求める。
【0109】
(具体的な函数の形)
さて、函数ρ(t)、u(p)、Rs(p)の具体例として、式(17)〜(19)をで表わされるものを用いることができる。
【0110】
【数44】
【0111】
さて、函数ρ(t)が、式(17)で表わされる場合には、式(3)は、式(20)のように表わされる。
【0112】
【数45】
【0113】
さらに、式(20)から、以下の式(21)が成立する。
【0114】
【数46】
【0115】
さらに、式(14)は、式(22)および(23)のように表わされ、式(16)は、式(24)のように表わされる。
【0116】
【数47】
【0117】
(離散化)
式(22)〜(24)をコンピュータで計算するために離散化を行なう。離散化の方法として、ここでは、計算を簡略化するためニストレーム法を用いる。
【0118】
積分区間の上限U、下限L、積分のサンプル数n、正則化パラメータaを用いて、式(25)〜(28)に従って、以下の変数h、xj、pj、wjを定義する。
【0119】
【数48】
【0120】
式(25)〜式(28)を用いて、式(22)を離散化すると式(29)が得られ、式(24)を離散化すると式(30)が得られる。
【0121】
【数49】
【0122】
したがって、式(29)の連立方程式の解{Hi n,a,t:i=0,1,2,・・・,n}を求めて、それらを式(30)に代入することによって、逆ラプラス変換による原像の数値解fa,s(n)(t)が得られる。
【0123】
(式(29)の連立方程式の解法)
式(29)の連立方程式の解を求める方法として、ここでは、以下の方法を用いる。
【0124】
まず、式(31)、(32)に示すように以下の変数qj、aijを用いる。このとき、式(33)、(34)が成り立つ。
【0125】
【数50】
【0126】
これらの変数を用いると、式(29)は、式(35)のように変形される。式(35)を行列の形に書き直すと、式(36)のようになり、係数行列Aを式(37)で定義することができる。
【0127】
【数51】
【0128】
係数行列AをLU分解し、前進代入および後退代入することによって、式(36)の解、つまり式(29)の解{Hin,a,t:i=0,1,2,・・・,n}を得ることができる。また、変数qjを用いると、式(30)は、式(38)のように変形することができる。
【0129】
【数52】
【0130】
(構成)
図1は、本発明の実施形態の逆ラプラス変換装置の構成を表わす図である。
【0131】
図1を参照して、この逆ラプラス変換装置1は、コンピュータを用いて実現されるものであって、入力部2と、メインメモリ6と、CPU3(Central Processing Unit)と、出力部13とを備える。メインメモリ6は、プログラム記憶部7と、変数記憶部8と、Hテーブル記憶部12と、三角行列記憶部9とを含む。CPU3は、テーブル作成部4および逆変換部5として機能する。
【0132】
入力部2は、サンプル数n、上限値U、下限値L、正則化パラメータa、軟化子パラメータs、ラプラス変換像F(pi)などを入力するものであり、たとえばキーボードなどで実現されている。
【0133】
出力部13は、逆ラプラス変換による原像の数値解fa,s(n)(t)を外部、たとえば表示装置14に出力する。表示装置14では、逆ラプラス変換による原像の数値解fa,s(n)(t)が文字またはグラフの形で表示される。
【0134】
プログラム記憶部7は、コンピュータを逆ラプラス変換装置1として機能させるための逆ラプラス変換プログラムが記憶されている。逆ラプラス変換プログラムは、USB(Universal Serial Bus)メモリ、CD−ROM(Compact Disk - Read Only Memory)などのコンピュータ読取り可能な記録媒体であるリムーマブルメディア99を通じて、外部からインストールすることができる。
【0135】
テーブル作成部4は、プログラム記憶部7内の逆ラプラス変換プログラムにしたがって、式(22)および式(23)で表わされる第二種積分方程式を離散化して得られる式(35)の連立方程式の解を数値計算で求め、数値計算に基づいて得られる第二種積分方程式の数値解{Hi n,a,t:i=0,1,2,・・・,n}をHテーブルに記述する。テーブル作成部4は、プログラム記憶部7内の逆ラプラス変換プログラムにしたがって、式(37)で示される係数行列Aを上三角行列Uと下三角行列Lの積に分解して、上三角行列Uと下三角行列Lの要素を三角行列記憶部9に記憶させ、計算範囲内のすべてのtについて、三角行列記憶部9に記憶されている上三角行列Uと下三角行列Lの要素の値に基づいて、式(35)の解を算出する。
【0136】
Hテーブル記憶部12は、テーブル作成部4で作成された図2で示されるようなHテーブルを記憶する。
【0137】
逆変換部5は、Hテーブルを参照して、式(24)を離散化した式(38)を数値計算することにより、逆ラプラス変換による原像の数値解fa,s(n)(t)を算出する。
【0138】
変数記憶部8は、テーブル作成部4によるテーブル作成処理および逆変換部5による逆ラプラス変換による原像の数値解を算出するのに必要な数値計算の変数の値を記憶する。
【0139】
三角行列記憶部9は、テーブル作成部4において、式(35)の連立方程式の解を計算する際に作成された下三角行列Lおよび上三角行列Uを記憶する。
【0140】
(Hテーブルの作成手順)
図3は、第1の実施形態におけるHテーブルを作成する手順を表わすフローチャートである。
【0141】
まず、外部から、入力部2にサンプル数n、上限値U、下限値L、正則化パラメータaが入力される(ステップS101)。
【0142】
次に、テーブル作成部4は、h=(U−L)/nを計算して、hの値を変数記憶部8に保存する(ステップS102)。
【0143】
次に、テーブル作成部4は、j=0、1、2、・・・、nについて、xj=L+j×hを計算して、xjの値を変数記憶部8に保存する(ステップS103)。
【0144】
次に、テーブル作成部4は、j=0、1、2、・・・、nについて、pj=exp{(π/2)×sinh(xj)}を計算して、pjの値を変数記憶部8に保存する(ステップS104)。
【0145】
次に、テーブル作成部4は、j=0、1、2、・・・、nについて、qj=pj×cosh(xj)を計算して、qjの値を変数記憶部8に保存する(ステップS105)。
【0146】
次に、テーブル作成部4は、i=0、1、2、・・・、n、およびj=0、1、2、・・・、nについて、aij=(π/2)×(2/(pi+pj)3)×{1+(pi+pj)+(pi+pj)2/2}×exp(−pj−1/pj)を計算して、aijの値を変数記憶部8に保存する(ステップS106)。
【0147】
次に、テーブル作成部4は、i=0、1、2、・・・、n、およびj=0、1、2、・・・、nについて係数行列Aの成分を計算する。テーブル作成部4は、i=jのときには、係数行列Aのij成分であるa+qi×aiiを計算し、i≠jのときには、係数行列Aのij成分であるqj×aijを計算する(ステップS107)。
【0148】
次に、テーブル作成部4は、係数行列AをLU分解して、上三角行列Uおよび下三角行列Lの各要素の値を三角行列記憶部9に保存する(ステップS108)。
【0149】
次に、テーブル作成部4は、t=t0(初期値)とする(ステップS109)。
次に、テーブル作成部4は、H(pi,t)=(2/pi3)×[1+pi+pi2/2−exp(−tpi)×{1+pi(t+1)+pi2(t+1)2/2}]を計算し、H(pi,t)を変数記憶部8に保存する(ステップS110)。
【0150】
次に、テーブル作成部4は、三角行列記憶部9に記憶されている上三角行列Uおよび下三角行列Lの要素の値をもとに、前進代入および後退代入によって、式(29)の解{Hin,a,t:i=0,1,2,・・・,n}を求めて、変数記憶部8内に記憶する(ステップS111)。
【0151】
次に、テーブル作成部4は、変数記憶部8内の式(29)の解{Hin,a,t:i=0,1,2,・・・,n}をHテーブルに記憶する(ステップS112)。
【0152】
次に、テーブル作成部4は、tがtm(終了値)と等しければ(ステップS113でYES)終了し、tがtm(終了値)と異なるならば(ステップS113でNO)、tをΔtだけ増加させて(ステップS114)、ステップS110からの処理を繰返す。
【0153】
(Hテーブルを用いた逆ラプラス変換の数値解の算出手順)
図4は、第1の実施形態におけるHテーブルを用いて逆ラプラス変換の数値解を算出する手順を表わすフローチャートである。
【0154】
まず、外部から、入力部2に軟化子パラメータsが入力される(ステップS200)。
次に、逆変換部5は、t=t0に設定する(ステップS201)。
【0155】
次に、逆変換部5は、j=0、SUM=0に設定する(ステップS202)。
次に、逆変換部5は、HテーブルからHjn,a,tを読出し、変数記憶部8からpj、qj、hを読み出す(ステップS203)。
【0156】
次に、外部から、入力部2にラプラス変換像F(pj)が入力される(ステップS204)。
【0157】
次に、逆変換部5は、SUM=SUM+(π/2)×h×F(pj)×(1/s2pj2){exp(−spj)−1}2×pj×Hjn,a,t ×qj×exp(−pj−1/pj)を計算し、計算結果を変数記憶部8に保存する(ステップS205)。
【0158】
次に、逆変換部5は、jがnと異なるならば(ステップS206でNO)、jを1だけ増加させて(ステップS207)、ステップS203からの処理を繰返す。
【0159】
逆変換部5は、jがnと等しければ(ステップS206でYES)、逆ラプラス変換値fa,s(n)(t)=SUMとする(ステップS208)。
【0160】
次に、逆変換部5は、tがtmと等しければ(ステップS209でYES)終了し、tがtm(終了値)と異なるならば(ステップS209でNO)、tをΔtだけ増加させて(ステップS210)、ステップS202からの処理を繰返す。
【0161】
以上のように、第1の実施形態の逆ラプラス変換装置によれば、原像の導函数が増大する場合、または原像が原点で零でない場合においても、逆ラプラス変換の数値解を計算することができる。また、式(36)の連立方程式の係数行列Aがtに依存しない性質を利用して、係数行列Aを一度だけLU分解し、分解して得られた下三角行列Lと上三角行列Uの各要素の値を記憶しておき、すべてのtについて、この記憶してある下三角行列Lと上三角行列Uの要素の値を用いて、連立方程式の解を算出するので、tが変わるごとに係数行列AをLU分解するのに比べて計算量を削減することができる。
【0162】
また、連立方程式の解{Hin,a,t:i=0,1,2,・・・,n}は、式(29)を参照すると、ラプラス変換像F(pi)に依存しない。したがって、一度だけ式(36)の解を計算しておき保存しておけば、どのようなF(pi)の逆ラプラス変換を求める場合にも利用することができるという特徴がある。
【0163】
[第1の実施形態の変形例]
本発明の実施形態では、第二種積分方程式を離散化して得られる連立方程式を構成する係数行列AをLU分解することによって連立方程式を解いたが、これに限定するものではなく、係数行列Aの逆行列A-1を計算して、連立方程式を解いてもよい。
【0164】
図5は、第1の実施形態の変形例の逆ラプラス変換装置の構成を表わす図である。
図5を参照して、この逆ラプラス変換装置81は、第1の実施形態のラプラス変換装置1と相違する点は、テーブル作成部84と、逆行列記憶部89である。
【0165】
テーブル作成部84は、プログラム記憶部7内の逆ラプラス変換プログラムにしたがって、式(22)および式(23)で表わされる第二種積分方程式を離散化して得られる式(35)の連立方程式の解を数値計算で求め、数値計算に基づいて得られる第二種積分方程式の数値解{Hi n,a,t:i=0,1,2,・・・,n}をHテーブルに書込む。テーブル作成部84は、プログラム記憶部7内の逆ラプラス変換プログラムにしたがって、式(37)で示される係数行列Aの逆行列A-1を計算して、逆行列A-1の要素の値を逆行列記憶部89に記憶させ、計算範囲内のすべてのtについて、逆行列記憶部89に記憶されている逆行列A-1の要素の値に基づいて、式(35)の解を算出する。
【0166】
逆行列記憶部89は、テーブル作成部4において、式(35)の連立方程式の解を計算する際に作成された逆行列A-1の各要素の値を記憶する。
【0167】
(Hテーブルの作成手順)
図6は、第1の実施形態の変形例おけるHテーブルを作成する手順を表わすフローチャートである。
【0168】
まず、外部から、入力部2にサンプル数n、上限値U、下限値L、正則化パラメータaが入力される(ステップS501)。
【0169】
次に、テーブル作成部4は、h=(U−L)/nを計算して、hの値を変数記憶部8に保存する(ステップS502)。
【0170】
次に、テーブル作成部4は、j=0、1、2、・・・、nについて、xj=L+j×hを計算して、xjの値を変数記憶部8に保存する(ステップS503)。
【0171】
次に、テーブル作成部4は、j=0、1、2、・・・、nについて、pj=exp{(π/2)×sinh(xj)}を計算して、pjの値を変数記憶部8に保存する(ステップS504)。
【0172】
次に、テーブル作成部4は、j=0、1、2、・・・、nについて、qj=pj×cosh(xj)を計算して、qjの値を変数記憶部8に保存する(ステップS505)。
【0173】
次に、テーブル作成部4は、i=0、1、2、・・・、n、およびj=0、1、2、・・・、nについて、aij=(π/2)×(2/(pi+pj)3)×{1+(pi+pj)+(pi+pj)2/2}×exp(−pj−1/pj)を計算して、aijの値を変数記憶部8に保存する(ステップS506)。
【0174】
次に、テーブル作成部4は、i=0、1、2、・・・、n、およびj=0、1、2、・・・、nについて係数行列Aの成分を計算する。テーブル作成部4は、i=jのときには、係数行列Aのij成分であるa+qi×aiiを計算し、i≠jのときには、係数行列Aのij成分であるqj×aijを計算する(ステップS507)。
【0175】
次に、テーブル作成部4は、係数行列Aの逆行列A-1を計算して、逆行列A-1の各要素の値を逆行列記憶部89に保存する(ステップS508)。
【0176】
次に、テーブル作成部4は、t=t0(初期値)とする(ステップS509)。
次に、テーブル作成部4は、H(pi,t)=(2/pi3)×[1+pi+pi2/2−exp(−tpi)×{1+pi(t+1)+pi2(t+1)2/2}]を計算し、H(pi,t)を変数記憶部8に保存する(ステップS510)。
【0177】
次に、テーブル作成部4は、逆行列記憶部89に記憶されている逆行列A-1の要素の値をもとに、式(29)の解{Hin,a,t:i=0,1,2,・・・,n}を求めて、変数記憶部8内に記憶する(ステップS511)。
【0178】
次に、テーブル作成部4は、変数記憶部8内の式(29)の解{Hin,a,t:i=0,1,2,・・・,n}をHテーブルに記憶する(ステップS512)。
【0179】
次に、テーブル作成部4は、tがtm(終了値)と等しければ(ステップS513でYES)終了し、tがtm(終了値)と異なるならば(ステップS513でNO)、tをΔtだけ増加させて(ステップS514)、ステップS510からの処理を繰返す。
【0180】
以上のように、第1の実施形態の変形例の逆ラプラス変換装置によれば、式(36)の連立方程式の係数行列Aがtに依存しない性質を利用して、係数行列Aの逆行列A-1を一度だけ計算して、逆行列A-1の要素の値を記憶しておき、すべてのtについて、この記憶してある逆行列A-1の要素の値を用いて、連立方程式の解を算出するので、tが変わるごとに係数行列Aの逆行列A-1を計算するのに比べて計算量を削減することができる。
【0181】
[第2の実施形態]
第2の実施形態は、式(29)を変形することによって、係数行列Aの代わりに対称行列A′を用いて、第二種積分方程式を離散化して連立方程式の解を算出する逆ラプラス変換装置に関する。
【0182】
(式(23)の連立方程式の解法)
第1の実施形態と同様に、式(31)〜(34)で表わされる変数qj、aijを用いる。
【0183】
第2の実施形態では、式(29)を式(39)のように変形する。
【0184】
【数53】
【0185】
式(40)のように、Hi′n,a,tを定義すると、式(39)は、式(41)のように変形される。式(41)を行列の形に書き直すと、式(42)のようになり、行列A′を式(43)のように定義することができる。
【0186】
【数54】
【0187】
行列A′は、対称行列であり、行列A′をコレスキー分解することによって、式(42)の解を得ることができる。式(42)の解{Hi′n,a,t:i=0,1,2,・・・,n}から、式(44)を用いて、式(29)の解{Hin,a,t:i=0,1,2,・・・,n}が得られる。
【0188】
【数55】
【0189】
(Hテーブルの作成手順)
図7は、第2の実施形態におけるHテーブルを作成する手順を表わすフローチャートである。
【0190】
まず、外部から、入力部2にサンプル数n、上限値U、下限値L、正則化パラメータaが入力される(ステップS301)。
【0191】
次に、テーブル作成部4は、h=(U−L)/nを計算して、hの値を変数記憶部8に保存する(ステップS302)。
【0192】
次に、テーブル作成部4は、j=0、1、2、・・・、nについて、xj=L+j×hを計算して、xjの値を変数記憶部8に保存する(ステップS303)。
【0193】
次に、テーブル作成部4は、j=0、1、2、・・・、nについて、pj=exp{(π/2)×sinh(xj)}を計算して、pjの値を変数記憶部8に保存する(ステップS304)。
【0194】
次に、テーブル作成部4は、j=0、1、2、・・・、nについて、qj=pj×cosh(xj)を計算して、qjの値を変数記憶部8に保存する(ステップS305)。
【0195】
次に、テーブル作成部4は、i=0、1、2、・・・、n、およびj=0、1、2、・・・、nについて、aij=(π/2)×(2/(pi+pj)3)×{1+(pi+pj)+(pi+pj)2/2}×exp(−pj−1/pj)を計算して、aijの値を変数記憶部8に保存する(ステップS306)。
【0196】
次に、テーブル作成部4は、i=0、1、2、・・・、n、およびj=0、1、2、・・・、nについて係数行列Aの成分を計算する。テーブル作成部4は、i=jのときには、係数行列Aのij成分であるa/qi×aiiを計算し、i≠jのときには、係数行列Aのij成分をaijとする(ステップS307)。
【0197】
次に、テーブル作成部4は、係数行列Aをコレスキー分解して、下三角行列L、および上三角行列LT(Tは、転置を表わす)の各要素を三角行列記憶部9に保存する(ステップS308)。
【0198】
次に、テーブル作成部4は、t=t0(初期値)とする(ステップS309)。
次に、テーブル作成部4は、H(pi,t)=(2/pi3)×[1+pi+pi2/2−exp(−tpi)×{1+pi(t+1)+pi2(t+1)2/2}]を計算し、H(pi,t)を変数記憶部8に保存する(ステップS310)。
【0199】
次に、テーブル作成部4は、三角行列記憶部9に記憶されている下三角行列Lおよび上三角行列LTをもとに、前進代入および後退代入によって、式(41)の解{Hi′n,a,t:i=0,1,2,・・・,n}を求めて、変数記憶部8に記憶する(ステップS311)。
【0200】
次に、テーブル作成部4は、i=0、1、2、・・・、nについて、Hi n,a,t=Hi′n,a,t/qiを計算して、変数記憶部8に記憶する(ステップS312)。
【0201】
次に、テーブル作成部4は、変数記憶部8内の式(29)の解である{Hin,a,t:i=0,1,2,・・・,n}をHテーブルに記憶する(ステップS313)。
【0202】
次に、テーブル作成部4は、tがtm(終了値)と等しければ(ステップS314でYES)終了し、tがtm(終了値)と異なるならば(ステップS314でNO)、tをΔtだけ増加させて(ステップS315)、ステップS310からの処理を繰返す。
【0203】
以上のように、第2の実施形態の逆ラプラス変換装置によれば、第1の実施形態と同様の効果に加えて、さらに第二種積分方程式を離散化して得られる連立方程式を構成する係数行列Aを対称行列A′に変換し、コレスキー分解を利用して連立方程式の解を算出することができるので、第1の実施形態よりもさらに計算量を削減できる。
【0204】
[第3の実施形態]
第3の実施形態は、第1の実施形態または第2の実施形態における、連立方程式の解{Hi n,a,t:i=0,1,2,・・・,n}および逆ラプラス変換の数値解fa,s(n)(t)を求める演算を多倍長で行なう装置に関する。以下では、第2の実施形態の演算を多倍長で行なう場合について説明するが、第1の実施形態の演算でも同様である。
【0205】
(構成)
図8は、第3の実施形態のCPU3の具体的な構成を表わす図である。
【0206】
図8を参照して、CPU3は、制御部53と、演算器54と、汎用レジスタ群56と、フラグレジスタ57とを備える。図1では、CPU3は、テーブル作成部4および逆変換部5として機能するとして説明した。図1と図8の関係を説明すると、テーブル作成部4および逆変換部5の機能を具体的に実現するのが、制御部53と、演算器54と、汎用レジスタ群56と、フラグレジスタ57である。つまり、テーブル作成部4および逆変換部5とは、共通に、制御部53と、演算器54と、汎用レジスタ群56と、フラグレジスタ57とを備える。言い換えると、制御部53と、演算器54と、汎用レジスタ群56と、フラグレジスタ57とは、プログラム記憶部7に格納された逆ラプラス変換プログラムにしたがって、第1または第2の実施形態で説明したHテーブルの作成と、逆ラプラス変換の数値解の算出を行なう。
【0207】
演算器54は、加算、減算、乗算、除算、ビットの論理演算、およびビットシフト演算などを行なう。
【0208】
汎用レジスタ群56は、16個の汎用レジスタRAX、RBX、RCX、RDX、RSI、RDI、RBP、RSP、R8、R9、R10、R11、R12、R13、R14、R15を含む。図8において、かっこ内は、アセンブラ言語において指定されるレジスタの名前を表わす。各レジスタは、64ビット長である。
【0209】
図9は、フラグレジスタ57の具体的な構成を表わす図である。
図9を参照して、フラグレジスタ57の最下位ビットに、直前の演算において、汎用レジスタ群56の各汎用レジスタの最上位ビット(64ビット目)において、次の上位ビットへの繰上りまたは繰下がりが生じたか否かを表わすキャリーフラグCFの値が保持される。キャリーフラグCFの値が「1」の場合には、繰上りまたは繰下がりが生じたことを示し、キャリーフラグCFの値が「0」の場合には、繰上りまたは繰下がりが生じなかったことを示す。
【0210】
演算器54による加算には、繰上がり付き加算と繰上りなし加算とがある。X+Yを繰上り加算する場合には、演算器54は、XとYとキャリーフラグCFの値とを加算し、加算結果に繰上りが生じた場合には、キャリーフラグCFの値を「1」にし、加算結果に繰上りが生じなかった場合には、キャリーフラグCFの値を「0」にする。一方、X+Yを繰上りなし加算する場合には、演算器54は、XとYとを加算し、加算結果に繰上りが生じた場合には、キャリーフラグCFの値を「1」にし、加算結果に繰上りが生じなかった場合には、キャリーフラグCFの値を「0」にする。
【0211】
演算器54による減算にも、繰下がり付き減算と繰下がりなし減算とがある。X−Yを繰下がり減算する場合には、演算器54は、XからYおよびキャリーフラグCFの値を減算し、減算結果に繰下がりが生じた場合には、キャリーフラグCFの値を「1」にし、減算結果に繰下がりが生じなかった場合には、キャリーフラグCFの値を「0」にする。一方、X−Yを繰下がりなし減算する場合には、演算器54は、XからYを減算し、減算結果に繰下がりが生じた場合には、キャリーフラグCFの値を「1」にし、減算結果に繰下がりが生じなかった場合には、キャリーフラグCFの値を「0」にする。
【0212】
図8の逆ラプラス変換装置で演算の対象となる数値、たとえば、n、U、L、a、h、s、xj、pj、qj、aij、係数行列Aの要素、係数行列Aをコレスキー分解したLおよびLT、H(pi,t)、Hi′n,a,t、Hin,a,t、F(pj)、S、fa,s(n)(t)などの数値は、浮動小数点形式で表わされる。
【0213】
図10(a)、(b)、(c)は、図8の逆ラプラス変換装置において演算の対象となる数値の表現を説明するための図である。
【0214】
図10(a)に示されるように、数値は、浮動小数点形式で表現される。符号部が「0」の場合は正の数を表わし、符号部が「1」の場合は負の数を表わす。仮数部のみが多倍長64nビット(nは2以上の自然数)で表わされる。指数部および符号部は、64ビットで表わされる。
【0215】
図10(b)、(c)に示されるように、仮数部の整数部の値は、通常は「1」であるが、加算および減算時には、小数点の位置を合わせるために仮数部のビットがシフトされる。仮数部の小数点の位置がKビットだけ左にシフトすると、指数部の値がKだけ減少する。逆に、仮数部の小数点の位置がKビットだけ右にシフトすると、指数部の値がKだけ増加する。
【0216】
図10(b)、(c)に示されるように、汎用レジスタ群56の64ビットの各汎用レジスタを用いて演算を行なうために、メインメモリ6に要素数nの配列aが確保される。仮数部は、64ビットごとに区切られて、区切られた各64ビットが配列aの各要素に格納される。たとえば、最上位から64個のビットが配列の要素a[0]に格納され、次の64個のビットが配列の要素a[1]に格納される。下位側のビットを格納する配列から、順次演算が行なわれる。
【0217】
(本発明の実施形態の加算処理)
本発明の実施形態における仮数部の加算処理について説明する。
【0218】
図11は、フラグレジスタ57を備えるCPU3の一例であるAMD64の加算処理の内容を表わす図である。同図は、C言語からadd(c,a,b,n)で呼び出されるアセンブラルーチンを表わす図である。
【0219】
add(c,a,b,n)は、メインメモリ6内の要素数nの配列aと、メインメモリ6内の要素数nの配列bとを加算して、加算結果をメインメモリ6内の要素数nの配列cに蓄積するルーチンである。
【0220】
add(c,a,b,n)が呼び出されると、制御部53は、カウンタiの値を保持するためのレジスタ%rcxにnの値を設定し、レジスタ%rsiに配列aへのポインタを設定し、レジスタ%rdxに配列bへのポインタを設定し、レジスタ%rdiに配列cへのポインタを設定する。
【0221】
2行目では、制御部53は、最上位ビットでの繰上りの有無を表わす情報を保持するためのレジスタ%raxの値を「0」(つまり、繰上りなし)に設定するとともに、フラグレジスタ57内のキャリーフラグCFの値を「0」(つまり、繰上りなし)に初期化する。
【0222】
6行目では、制御部53は、メインメモリ6内のa[i−1]の値をレジスタ%r8にロードする。
【0223】
7行目では、制御部53は、メインメモリ6内のb[i−1]の値をレジスタ%r10にロードする。
【0224】
8行目では、演算器54は、繰上り付き加算を実行する。すなわち、演算器54は、レジスタ%r8内のデータとレジスタ%r10内のデータと、フラグレジスタ57内のキャリーフラグCFの値とを加算し、加算結果をレジスタ%r10に保存するとともに、加算により繰上りがあった場合に、フラグレジスタ57内のキャリーフラグCFの値を「1」に設定し、繰上りがなかった場合に、フラグレジスタ57内のキャリーフラグCFの値を「0」に設定する。
【0225】
9行目では、制御部53は、レジスタ%r10内の加算結果を、メインメモリ6内のc[i−1]の領域に保存する。
【0226】
10行目では、演算器54は、カウンタiを「1」だけ減らす。
11行目では、制御部53は、カウンタiが「0」の場合には、ループを抜け、カウンタiが「0」以外の場合には、ループの先頭に戻るように制御する。
【0227】
13行目では、演算器54は、即値「0」と、レジスタ%rax内のデータ(=「0」)と、フラグレジスタ57内のキャリーフラグCFの値とを加算し、加算結果をレジスタ%raxに保存する。したがって、レジスタ%raxには、最上位ビットでの繰上りの有無を表わす情報が保持される。
【0228】
(参考:キャリーフラグを利用しない加算処理)
次に、比較のために、フラグレジスタ57を備えないCPU3の一例であるAlphaにおける仮数部の加算処理について説明する。
【0229】
図12は、Alphaの加算処理の内容を表わす図である。同図は、C言語からadd(c,a,b,n)で呼び出されるアセンブラルーチンを表わす図である。
【0230】
add(c,a,b,n)は、メインメモリ6内の要素数nの配列aと、メインメモリ6内の要素数nの配列bとを加算して、加算結果をメインメモリ6内の要素数nの配列cに蓄積するルーチンである。
【0231】
add(c,a,b,n)が呼び出されると、制御部53は、カウンタiの値を保持するためのレジスタ$19にnの値を設定し、レジスタ$17に配列aへのポインタを設定し、レジスタ$18に配列bへのポインタを設定し、レジスタ$16に配列cへのポインタを設定する。
【0232】
2行目では、制御部53は、配列a[n]のアドレスを計算して、それをレジスタ$17に格納する。
【0233】
3行目では、制御部53は、配列b[n]のアドレスを計算して、それをレジスタ$18に格納する。
【0234】
4行目では、制御部53は、配列c[n]のアドレスを計算して、それをレジスタ$16に格納する。
【0235】
6行目では、制御部53は、繰上りの有無を表わす情報を保持するためのレジスタ$0の値を「0」(つまり、繰上りなし)に設定する。
【0236】
9行目では、制御部53は、メインメモリ6内のa[i−1]の値をレジスタ$1にロードする。
【0237】
10行目では、制御部53は、メインメモリ6内のb[i−1]の値をレジスタ$2にロードする。
【0238】
12行目では、演算器54は、加算を実行する。すなわち、演算器54は、レジスタ$1内のデータとレジスタ$2内のデータとを加算し、加算結果をレジスタ$5に保存する。
【0239】
13行目では、演算器54は、12行目の加算により繰上りがあったか否かを調べる。すなわち、制御部53は、レジスタ$5内のデータ(=加算結果c)がレジスタ$2内のデータ(=b)よりも小さい場合には、繰上りがあったことがわかるので、レジスタ$23に「1」を格納し、それ以外の場合には、繰上りがないことがわかるので、レジスタ$23に「0」を格納する。
【0240】
15行目では、演算器54は、レジスタ$5内のデータ(=加算結果)と、レジスタ$0内に格納されている1回前のループでの繰上りの有無を表わすデータとを加算して、加算結果をレジスタ$6に格納する。
【0241】
16行目では、演算器54は、15行目の加算により繰上りがあったか否かを調べる。すなわち、制御部53は、レジスタ$6内のデータがレジスタ$0内のデータよりも小さい場合には、繰上りがあったことがわかるので、レジスタ$0に「1」を格納し、それ以外の場合には、繰上りがないことがわかるので、レジスタ$0に「0」を格納する。
【0242】
18行目では、制御部53は、レジスタ$6内の加算結果を、メインメモリ6内のc[i−1]の領域に保存する。
【0243】
19行目では、演算器54は、レジスタ$0内のデータと、レジスタ$23内のデータとの論理和をレジスタ$0に格納することにより、13行目の繰上り判定の結果と16行目の繰上り判定の結果とを統合する。このようにして判定結果が統合できるのは、レジスタ$0内のデータとレジスタ$23内のデータの両方が「1」となる場合は、ありえないからである。
【0244】
21行目では、制御部53は、配列a[i−2]のアドレスを計算して、それをレジスタ$17に格納する。
【0245】
22行目では、制御部53は、配列b[i−2]のアドレスを計算して、それをレジスタ$18に格納する。
【0246】
23行目では、制御部53は、配列c[i−2]のアドレスを計算して、それをレジスタ$16に格納する。
【0247】
25行目では、演算器54は、カウンタiを「1」だけ減らす。
26行目では、制御部53は、カウンタiが「0」の場合には、ループを抜け、カウンタiが「0」以外の場合には、ループの先頭に戻るように制御する。
【0248】
(加算処理の比較)
図11と図12を比較すると、図11の8行目の繰上り付き加算命令adcqと同等の機能を実現するために、図12では、12行目の加算命令addqと、13行目の比較命令cmplt、15行目の加算命令addqと、16行目の比較命令cmpltと、19行目の論理和命令orの5つの命令が必要となる。以上のことから、本発明の実施形態のようにキャリーフラグCFを利用することにより、それを利用しない場合と比べて、仮数部の加算処理の命令数を削減することができることがわかる。
【0249】
(減算処理)
本発明の実施形態における仮数部の減算処理について説明する。
【0250】
図13は、フラグレジスタ57を備えるCPU3の一例であるAMD64の減算処理の内容を表わす図である。同図は、C言語からsub(c,a,b,n)で呼び出されるアセンブラルーチンを表わす図である。
【0251】
sub(c,a,b,n)は、メインメモリ6内の要素数nの配列aから、メインメモリ6内の要素数nの配列bを減算して、減算結果をメインメモリ6内の要素数nの配列cに蓄積するルーチンである。
【0252】
sub(c,a,b,n)が呼び出されると、制御部53は、カウンタiの値を保持するためのレジスタ%rcxにnの値を設定し、レジスタ%rsiに配列aへのポインタを設定し、レジスタ%rdxに配列bへのポインタを設定し、レジスタ%rdiに配列cへのポインタを設定する。
【0253】
2行目では、制御部53は、最上位ビットでの繰下がりの有無を表わす情報を保持するためのレジスタ%raxの値を「0」(つまり、繰下がりなし)に設定するとともに、フラグレジスタ57内のキャリーフラグCFの値を「0」(つまり、繰下がりなし)に初期化する。
【0254】
6行目では、制御部53は、メインメモリ6内のa[i−1]の値をレジスタ%r8にロードする。
【0255】
7行目では、制御部53は、メインメモリ6内のb[i−1]の値をレジスタ%r10にロードする。
【0256】
8行目では、演算器54は、繰下がり付き減算を実行する。すなわち、演算器54は、レジスタ%r8内のデータから、レジスタ%r10内のデータおよびフラグレジスタ57内のキャリーフラグCFの値とを減算し、減算結果をレジスタ%r10に保存するとともに、減算により繰下がりがあった場合に、フラグレジスタ57内のキャリーフラグCFの値を「1」に設定し、繰下がりがなかった場合に、フラグレジスタ57内のキャリーフラグCFの値を「0」に設定する。
【0257】
9行目では、制御部53は、レジスタ%r10内の加算結果を、メインメモリ6内のc[i−1]の領域に保存する。
【0258】
10行目では、演算器54は、カウンタiを「1」だけ減らす。
11行目では、制御部53は、カウンタiが「0」の場合には、ループを抜け、カウンタiが「0」以外の場合には、ループの先頭に戻るように制御する。
【0259】
13行目では、演算器54は、即値「0」と、レジスタ%rax内のデータ(=「0」)と、フラグレジスタ57内のキャリーフラグCFの値とを加算し、加算結果をレジスタ%raxに保存する。したがって、レジスタ%raxには、最上位ビットでの繰り下がりの有無を表わす情報が保持される。
【0260】
(参考:キャリーフラグを利用しない減算処理)
次に、比較のために、フラグレジスタ57を備えないCPU3の一例であるAlphaにおける仮数部の減算処理について説明する。
【0261】
図14は、Alphaの減算処理の内容を表わす図である。同図は、C言語からsub(c,a,b,n)で呼び出されるアセンブラルーチンを表わす図である。
【0262】
sub(c,a,b,n)は、メインメモリ6内の要素数nの配列aから、メインメモリ6内の要素数nの配列bと減算して、減算結果をメインメモリ6内の要素数nの配列cに蓄積するルーチンである。
【0263】
sub(c,a,b,n)が呼び出されると、制御部53は、カウンタiの値を保持するためのレジスタ$19にnの値を設定し、レジスタ$17に配列aへのポインタを設定し、レジスタ$18に配列bへのポインタを設定し、レジスタ$16に配列cへのポインタを設定する。
【0264】
2行目では、制御部53は、配列a[n]のアドレスを計算して、それをレジスタ$17に格納する。
【0265】
3行目では、制御部53は、配列b[n]のアドレスを計算して、それをレジスタ$18に格納する。
【0266】
4行目では、制御部53は、配列c[n]のアドレスを計算して、それをレジスタ$16に格納する。
【0267】
6行目では、制御部53は、繰下がりの有無を表わす情報を保持するためのレジスタ$0の値を「0」(つまり、繰下がりなし)に設定する。
【0268】
9行目では、制御部53は、メインメモリ6内のa[i−1]の値をレジスタ$1にロードする。
【0269】
10行目では、制御部53は、メインメモリ6内のb[i−1]の値をレジスタ$2にロードする。
【0270】
12行目では、演算器54は、減算を実行する。すなわち、演算器54は、レジスタ$1内のデータからレジスタ$2内のデータを減算し、減算結果をレジスタ$5に保存する。
【0271】
13行目では、演算器54は、12行目の減算により繰下がりがあったか否かを調べる。すなわち、制御部53は、レジスタ$1内のデータ(=a)がレジスタ$2内のデータ(=b)よりも小さい場合には、繰下がりがあったことがわかるので、レジスタ$6に「1」を格納し、それ以外の場合には、繰下がりがないことがわかるので、レジスタ$6に「0」を格納する。
【0272】
15行目では、演算器54は、レジスタ$5内のデータ(=減算結果)から、レジスタ$0内に格納されている1回前のループでの繰下がりの有無を表わすデータを減算して、減算結果をレジスタ$7に格納する。
【0273】
16行目では、演算器54は、15行目の減算により繰下がりがあったか否かを調べる。すなわち、制御部53は、レジスタ$5内のデータがレジスタ$0内のデータよりも小さい場合には、繰下がりがあったことがわかるので、レジスタ$0に「1」を格納し、それ以外の場合には、繰下がりがないことがわかるので、レジスタ$0に「0」を格納する。
【0274】
18行目では、制御部53は、レジスタ$7内の減算結果を、メインメモリ6内のc[i−1]の領域に保存する。
【0275】
19行目では、演算器54は、レジスタ$0内のデータと、レジスタ$6内のデータとの論理和をレジスタ$0に格納することにより、13行目の繰下がり判定の結果と16行目の繰下がり判定の結果とを統合する。このようにして判定結果が統合できるのは、レジスタ$0内のデータとレジスタ$6内のデータの両方が「1」となる場合は、ありえないからである。
【0276】
21行目では、制御部53は、配列a[i−2]のアドレスを計算して、それをレジスタ$17に格納する。
【0277】
22行目では、制御部53は、配列b[i−2]のアドレスを計算して、それをレジスタ$18に格納する。
【0278】
23行目では、制御部53は、配列c[i−2]のアドレスを計算して、それをレジスタ$16に格納する。
【0279】
25行目では、演算器54は、カウンタiを「1」だけ減らす。
26行目では、制御部53は、カウンタiが「0」の場合には、ループを抜け、カウンタiが「0」以外の場合には、ループの先頭に戻るように制御する。
【0280】
(減算処理の比較)
図13と図14を比較すると、図13の8行目の繰下がり付き減算命令sbbqと同等の機能を実現するために、図14では、12行目の減算命令subqと、13行目の比較命令cmplt、15行目の減算命令subqと、16行目の比較命令cmpltと、19行目の論理和命令orの5つの命令が必要となる。以上のことから、本発明の実施形態のようにキャリーフラグCFを利用することにより、それを利用しない場合と比べて、加算部の減算処理の命令数を削減することができることがわかる。
【0281】
(乗算処理)
まず、本発明の実施形態で用いる乗算の基本的な処理内容を説明する。図15に示すように、要素数4の配列Aと、要素数4の配列Bとの乗算を考える。配列A[0]のデータがa0、配列A[1]のデータがa1、配列A[2]のデータがa2、配列A[3]のデータがa3とする。配列B[0]のデータがb0、配列B[1]のデータがb1、配列B[2]のデータがb2、配列B[3]のデータがb3とする。a0、a1、a2、a3、b0、b1、b2、b3のデータ長は、それぞれ64ビットである。図15に示すように、配列Aと、b0、b1、b2、b3それぞれとの乗算結果は、64×5ビットとなる。配列Aと配列Bとを乗算すると、乗算結果は64×8ビットとなる。
【0282】
浮動小数点表現の仮数部の積を求めることを考えた場合、このような64×8ビットがすべて必要になるわけではない。下位のビットは、実際上無視しても、計算精度に大きな影響はない。
【0283】
本発明の実施形態では、要素数nの配列A(つまり64×nビット)と、要素数nの配列B(つまり64×nビット)との乗算では、上位半分の64×nビットのみを求めて、要素数nの配列Cに格納する。
【0284】
図15で示される配列Aと配列Bの乗算の場合には、上位から64×4ビットだけを求める。この場合、A×b3の計算では、a0×b3の計算は必要だが、a1×b3、a2×b3、a3×b3の計算は省略することができる。一方、A×b0の計算では、a0×b0、a1×b0、a2×b0、a3×b0のすべての計算が必要となる。
【0285】
図16は、本発明の実施形態における乗算処理の内容を表わす図である。
図16を参照して、この乗算処理は、要素数nの配列A(各要素は64ビット)と、要素数nの配列B(各要素は64ビット)とを乗算して、乗算結果のうち上位半分の64×nビットを、要素数nの配列Cに格納する。
【0286】
8行目のmul_addは、配列Aの先頭から(n−i)個の要素(A[0]、・・・、A[n−i−1])と、配列B[i]とを乗算して、乗算結果を配列C′[i]から後ろに加えるルーチンである。たとえば、図17に示すように、nが4で、iが3の場合には、B[3](=b3)とA[0](=a0)の乗算だけが行なわれて、乗算結果が配列C′[3]から後ろに加えられる。
【0287】
(本発明の実施形態のmul_add)
図18は、フラグレジスタ57を備えるCPU3の一例であるAMD64のmul_addの内容を表わす図である。同図は、C言語から、mul_add(c,a,b,n)で呼び出されるアセンブラルーチンを表わす図である。
【0288】
mul_add(c,a,b,n)は、メインメモリ6内の要素数nの配列aと、メインメモリ6内の変数bとを乗算して、乗算結果をメインメモリ6内の要素数n+1の配列cに加算する命令である。
【0289】
mul_add(c,a,b,n)が呼び出されると、制御部53は、レジスタ%rdiに配列cへのポインタを設定し、レジスタ%rsiに配列aへのポインタを設定し、レジスタ%rdxに乗数bの値を設定し、カウンタiの値を保持するためのレジスタ%rcxにnの値を設定する。
【0290】
1行目では、制御部53は、繰上りの有無を表わす情報を保持するためのレジスタ%r9の値を「0」(つまり、繰上りなし)に設定するとともに、フラグレジスタ57内のキャリーフラグCFの値を「0」(つまり、繰上りなし)に初期化する。
【0291】
2行目では、制御部53は、レジスタ%rdx内の乗数bの値をレジスタ%r8に格納する。
【0292】
5行目では、制御部53は、メインメモリ6内のa[i−1]の値をレジスタ%raxにロードする。
【0293】
6行目では、制御部53は、レジスタ%rax内のデータ(=a[i−1])と、レジスタ%r8内のデータ(=b)とを乗算する。乗算結果の上位ビットはレジスタ%rdxに格納され、乗算結果の下位ビットはレジスタ%raxに格納される。
【0294】
9行目では、演算器54は、繰上りなし加算を実行する。すなわち、演算器54は、レジスタ%r9に格納されている前回のループによる繰上りの有無を表わすデータと、レジスタ%rdxに格納されている乗算結果の上位ビットとを加算して、加算結果をレジスタ%rdxに格納する。この加算では、繰上りが生じないので、フラグレジスタ57内のキャリーフラグCFは「0」となる。
【0295】
10行目では、演算器54は、繰上りなし加算を実行する。すなわち、演算器54は、レジスタ%raxに格納されている乗算結果の下位ビットと、メインメモリ6内の配列c[i]とを加算して、加算結果をメインメモリ6内の配列c[i]に格納する。この加算によって、繰上りが生じた場合は、キャリーフラグCFは「1」となり、繰上りが生じなかった場合は、キャリーフラグCFは「0」となる。
【0296】
11行目では、制御部53は、繰上りの有無を表わす情報を保持するためのレジスタ%r9の値を「0」(つまり、繰上りなし)に設定する。
【0297】
12行目では、演算器54は、繰上り付き加算を実行する。すなわち、演算器54は、フラグレジスタ57内のキャリーフラグCFの値と、レジスタ%rdxに格納されている乗算結果の上位ビットと、メインメモリ6内の配列c[i−1]とを加算して、加算結果をメインメモリ6内の配列c[i−1]に格納するとともに、加算により繰上りがあった場合に、フラグレジスタ57内のキャリーフラグCFの値を「1」に設定し、繰上りがなかった場合に、フラグレジスタ57内のキャリーフラグCFの値を「0」に設定する。
【0298】
13行目では、演算器54は、繰上り付き加算を実行する。すなわち、演算器54は、即値「0」と、レジスタ%r9内の繰上りの有無を表わす情報と、フラグレジスタ57内のキャリーフラグCFの値とを加算して、加算結果をレジスタ%r9に格納する。この加算では、繰上りが生じないので、フラグレジスタ57内のキャリーフラグCFは「0」となる。
【0299】
15行目では、演算器54は、カウンタiを「1」だけ減らす。
16行目では、制御部53は、カウンタiが「0」の場合には、ループを抜け、カウンタiが「0」以外の場合には、ループの先頭に戻るように制御する。
【0300】
(参考:キャリーフラグを利用しないmul_add)
次に、比較のために、フラグレジスタ57を備えないCPU3の一例であるAlphaにおけるmul_addについて説明する。
【0301】
図19は、Alphaのmul_addの内容を表わす図である。同図は、C言語から、mul_add(c,a,b,n)で呼び出されるアセンブラルーチンを表わす図である。
【0302】
mul_add(c,a,b,n)は、メインメモリ6内の要素数nの配列aと、メインメモリ6内の変数bとを乗算して、乗算結果をメインメモリ6内の要素数n+1の配列cに加算するルーチンである。
【0303】
mul_add(c,a,b,n)が呼び出されると、制御部53は、カウンタiの値を保持するためのレジスタ$19にnの値を設定し、レジスタ$17に配列aへのポインタを設定し、レジスタ$18に乗数bの値を設定し、レジスタ$16に配列cへのポインタを設定する。
【0304】
1行目では、制御部53は、繰上がりの有無を表わす情報を保持するためのレジスタ$23の値を「0」(つまり、繰上がりなし)に設定する。
【0305】
2行目では、制御部53は、配列c[n]のアドレスを計算して、それをレジスタ$16に格納する。
【0306】
3行目では、制御部53は、配列a[n]のアドレスを計算して、それをレジスタ$17に格納する。
【0307】
6行目では、制御部53は、メインメモリ6内のa[i−1]の値をレジスタ$1にロードする。
【0308】
7行目では、制御部53は、メインメモリ6内のc[i]をレジスタ$21にロードする。
【0309】
8行目では、制御部53は、メインメモリ6内のc[i−1]をレジスタ$20にロードする。
【0310】
10行目では、演算器54は、レジスタ$1内のa[i−1]とレジスタ$18内のbとを乗算して、乗算結果の下位64ビットの値をレジスタ$2に格納する。
【0311】
11行目では、演算器54は、レジスタ$1内のa[i−1]とレジスタ$18内のbとを乗算して、乗算結果の上位64ビットの値をレジスタ$1に格納する。
【0312】
13行目では、演算器54は、レジスタ$1内の乗算結果の上位64ビットと、レジスタ$23に格納されている1回前のループでの繰上りの有無を表わすデータとを加算して、加算結果をレジスタ$1に格納する。
【0313】
15行目では、演算器54は、レジスタ$2内の乗算結果の下位64ビットと、レジスタ$21内のc[i]とを加算して、レジスタ$21に格納する。
【0314】
16行目では、制御部53は、15行目の加算により繰上がりがあったか否かを調べる。すなわち、制御部53は、レジスタ$21内のデータ(=c[i])がレジスタ$2内のデータよりも小さい場合には、繰上がりがあったことがわかるので、レジスタ$2に「1」を格納し、それ以外の場合には、繰上がりがないことがわかるので、レジスタ$2に「0」を格納する。
【0315】
18行目では、演算器54は、レジスタ$1内の乗算結果の上位64ビット(13行目で繰上り処理されたもの)と、レジスタ$20内のc[i−1]とを加算して、レジスタ$20に格納する。
【0316】
19行目では、制御部53は、18行目の加算により繰上がりがあったか否かを調べる。すなわち、制御部53は、レジスタ$20内のデータ(=c[i−1])がレジスタ$1内のデータよりも小さい場合には、繰上がりがあったことがわかるので、レジスタ$1に「1」を格納し、それ以外の場合には、繰上がりがないことがわかるので、レジスタ$1に「0」を格納する。
【0317】
21行目では、演算器54は、レジスタ$2内の15行目の加算の繰上りの有無を表わすデータと、レジスタ$20内のc[i−1](上位64ビットが加算後のもの)とを加算して、レジスタ$20に格納する。
【0318】
22行目では、制御部53は、21行目の加算により繰上がりがあったか否かを調べる。すなわち、制御部53は、レジスタ$20内のデータ(=c[i−1])がレジスタ$2内のデータよりも小さい場合には、繰上がりがあったことがわかるので、レジスタ$2に「1」を格納し、それ以外の場合には、繰上がりがないことがわかるので、レジスタ$2に「0」を格納する。
【0319】
24行目では、制御部53は、レジスタ$21内のc[i]をメインメモリ6内の領域に保存する。
【0320】
25行目では、制御部53は、レジスタ$20内のc[i−1]をメインメモリ6内の領域に保存する。
【0321】
26行目では、演算器54は、レジスタ$1内のデータと、レジスタ$2内のデータとの論理和をレジスタ$23に格納することにより、19行目の繰上がり判定の結果と22行目の繰上がり判定の結果とを統合する。このようにして判定結果が統合できるのは、レジスタ$1内のデータとレジスタ$2内のデータの両方が「1」となる場合は、ありえないからである。
【0322】
28行目では、制御部53は、配列a[i−2]のアドレスを計算して、それをレジスタ$17に格納する。
【0323】
29行目では、制御部53は、配列c[i−2]のアドレスを計算して、それをレジスタ$16に格納する。
【0324】
31行目では、演算器54は、カウンタiを「1」だけ減らす。
32行目では、制御部53は、カウンタiが「0」の場合には、ループを抜け、カウンタiが「0」以外の場合には、ループの先頭に戻るように制御する。
【0325】
(mul_addおよび乗算処理の比較)
図18と図19とを比較すると、図18では、ループ内部の命令数が7であるのに対して、図19では、ループ内部の命令数が18である。以上のことから、本発明の実施形態のようにキャリーフラグCFを利用することにより、それを利用しない場合と比べて、配列と配列の1要素との乗算であるmul_addの命令数を削減することができることがわかる。その結果、キャリーフラグCFを利用することにより、仮数部を表わす配列と配列の乗算の命令数も削減できる。
【0326】
(除算処理)
本発明の実施形態では、古典的なニュートン法により除算を行なう。x=a/bを求めるために、まず、1/bを求め、その結果にaを乗ずる。1/bを求めるために、f(x)=(1/x)−bについてニュートン法を適用する。このとき、ニュートン法の漸化式は、たとえば、以下の式(45)、(46)で与えられる。
【0327】
x0=1 ・・・(45)
xn+1=xn−f(xn)/f′(xn)=xn×(2−b×xn) ・・・(46)
図20は、ニュートン法を用いた除算処理の内容を説明するための図である。図21は、図20中のサブルーチンcmpの内容を表わす図である。
【0328】
図20において、21行目は、q=b×xnの乗算を実行するルーチンである。
22行目は、q=(2−q)の減算を実行するルーチンである。
【0329】
23行目は、q=xn×qの乗算を実行するルーチンである。
25行目および26行目において、xn+1(=q)がxnと等しければ、20行目〜30行目のループを抜ける。
【0330】
32行目は、xn、すなわち(1/b)にaを乗ずるルーチンである。
(除算処理の比較)
図20の21行目、23行目、32行目の乗算ルーチンmulを実行するCPU3がAMD64の場合では、図16および図18に示されるものであるのに対して、CPU3がAlphaの場合では、図16および図19に示されるものである。両者を比較すると、CPU3がAMD64の方が、命令数が少ないことがわかる。
【0331】
また、図20の22行目の減算ルーチンsubを実行するCPU3がAMD64の場合では、図13に示されるものであるのに対して、CPU3がAlphaの場合では、図14に示されるものである。両者を比較すると、CPU3がAMD64の方が、命令数が少ないことがわかる。以上のことから、本発明の実施形態のようにキャリーフラグCFを利用することにより、それを利用しない場合と比べて、仮数部の除算処理の命令数を削減することができることがわかる。
【0332】
(計算の精度についての実験結果)
本発明の実施の形態では、仮数部が64×nビットで表わされるが、nを大きくすることで、fa,s(n)(t)の計算精度に影響する正則化パラメータaの値を小さくすることができる。以下では、正則化パラメータaの大きさによって、fa,s(n)(t)がどのように変化するかの実験結果を説明する。f(t)がステップ函数の場合に、ラプラス変換像F(pi)が解析的に求まるが、これを用いて、実験を行なった。
【0333】
図22(a)は、正則化パラメータa=10-1、10-4、10-8、10-12のときの、計算結果fa,s(n)(t)を示す図である。図22(b)は、正則化パラメータa=10-100、10-400のときの、計算結果fa,s(n)(t)を示す図である。
【0334】
これらの図から、正則化パラメータaを小さくするほど、fa,s(n)(t)の計算精度が増加することがわかる。図22内のステップ函数の場合には、正則化パラメータaは、少なくとも10-100ぐらい小さくなくては有効な逆ラプラス変換像が得られないことがわかる。正則化パラメータaの値をこのような小さな値にして、かつ現実的な時間内で逆ラプラス変換像を算出するためには、本発明の実施形態で説明したように、キャリーフラグCFを用いて命令数を大幅に削減した多倍長演算が必要不可欠となる。
【0335】
また、本発明の実施の形態では、式(6)を修正した式(10)を用いることによって、f(0)が0以外の函数についても、逆ラプラス変換の値を高精度で求めることができるという特徴がある。以下では、軟化子パラメータsの大きさによって、fa,s(n)(t)がどのように変化するかの実験結果を説明する。ここでは、式(47)、(48)で表される函数f(t)を実験に用いた。
【0336】
0<t<1において、f(t)=1 ・・・(47)
1<tにおいて、f(t)=0 (48)
図23(a)は、s=0.1での計算結果fa,s(n)(t)を示す図であり、図23(b)は、s=0.01での計算結果fa,s(n)(t)を示す図である。
【0337】
これらの図から、軟化子パラメータsを小さくするほど、fa,s(n)(t)の計算精度が増加することがわかる。軟化子パラメータsを目的に応じて適当な値に設定することによって、f(0)が0以外の函数について、逆ラプラス変換の値を高精度で求めることができる。
【0338】
[第4の実施形態]
第4の実施形態は、Hテーブルを作成する処理と、Hテーブルを用いて逆ラプラス変換の数値解を算出する処理とを別個の装置で行なう構成に関するものである。
【0339】
図24は、Hテーブルを作成する逆ラプラス変換用テーブル作成装置の構成を表わす図である。図24の逆ラプラス変換用テーブル作成装置は、コンピュータで実現されるものであり、図24において、図1と同一の構成要素については同一の符号を付している。
【0340】
図24を参照して、CPU63は、第1〜第3の実施形態と異なり、テーブル作成部4としてのみ機能する。
【0341】
入力部62は、第1〜第3の実施形態と異なり、Hテーブルの作成に必要なサンプル数n、上限値U、下限値L、正則化パラメータaが入力される。
【0342】
プログラム記憶部65には、第1〜第3の実施形態と異なり、コンピュータをテーブル作成部4として機能させるための逆ラプラス変換のためのテーブル作成プログラムが記憶されている。逆ラプラス変換のためのテーブル作成プログラムは、USB(Universal Serial Bus)メモリ、CD−ROM(Compact Disk - Read Only Memory)などのコンピュータ読取り可能な記録媒体であるリムーマブルメディア67を通じて、外部からインストールすることができる。
【0343】
変数記憶部66には、第1〜第3の実施形態と異なり、Hテーブル作成にのみ必要な変数が記憶される。
【0344】
メインメモリ64内のHテーブル記憶部12に記憶されているHテーブルと、変数記憶部66に記憶されている変数h、pi、qi(i=0,1,2,・・・,n)は、リムーマブルメディア67に出力することができる。これにより、リムーマブルメディア67に記憶されたHテーブルおよび変数を他の装置のメインメモリに送り記憶させることで、他の装置で、Hテーブルおよび変数h、pi、qiを利用することができる。
【0345】
図25は、逆ラプラス変換の数値解を算出する逆ラプラス変換の数値解算出装置の構成を表わす図である。図25の逆ラプラス変換の数値解算出装置は、コンピュータで実現されるものであり、図25において、図1と同一の構成要素については同一の符号を付している。
【0346】
図25を参照して、CPU73は、第1〜第3の実施形態と個なり、逆変換部5としてのみ機能する。
【0347】
入力部72は、第1〜第3の実施形態と異なり、逆ラプラス変換の数値解算出のために必要なラプラス変換像F(pi)と、軟化子パラメータsとが入力される。
【0348】
プログラム記憶部75には、第1〜第3の実施形態と異なり、コンピュータを逆変換部5および出力部13として機能させるための逆ラプラス変換の数値解算出プログラムが記憶されている。逆ラプラス変換の数値解算出プログラムは、USB(Universal Serial Bus)メモリ、CD−ROM(Compact Disk - Read Only Memory)などのコンピュータ読取り可能な記録媒体であるリムーマブルメディア77を通じて、外部からインストールすることができる。
【0349】
リムーマブルメディア77に記憶されているHテーブルと、変数h、pi、qi(i=0,1,2,・・・,n)は、それぞれメインメモリ74内のHテーブル記憶部12と、変数記憶部76に出力することができる。これにより、図24の装置で作成されてリムーマブルメディア67に記憶されたHテーブルおよび変数h、pi、qi(i=0,1,2,・・・,n)を、図25の装置で利用することができる。
【0350】
以上のように、第4の実施形態では、Hテーブルの作成と、逆ラプラス変換の数値解の算出とを別個の装置で行なうことができるので、Hテーブルの作成を高速の計算機で行なってユーザに提供し、ユーザは、自身のパソコン上で提供を受けたHテーブルを用いて逆ラプラス変換の数値解を算出を行なうようにすることができる。
【0351】
(変形例)
本発明は、上記の実施の形態に限定されるものではなく、たとえば以下のような変形例を含む。
【0352】
(1) テーブルに記述する情報
本発明の第4の実施形態では、図24の逆ラプラス変換用テーブル作成装置が、連立方程式の解{Hi n,a,t:i=0,1,2,・・・,n}を記述したHテーブルと、h、pi、qi(i=0,1,2,・・・,n)とをリムーバブルメディアに出力して、図25の逆ラプラス変換の数値解算出装置でこれを利用したが、これに限定するものではない。たとえば、逆ラプラス変換用テーブル作成装置が、Hテーブルだけをリムーバブルメディアに出力し、逆ラプラス変換の数値解算出装置では、n、U、Lを入力して、h、pi、qjを再計算するものとしてもよい。また、逆ラプラス変換用テーブル作成装置が、h×pi×Hi n,a,t×qiを計算して、計算結果を記述したテーブルとpiとをリムーバブルメディアに出力し、逆ラプラス変換の数値解算出装置で、これらを利用して、式(38)で表わされる計算をより簡単に計算するようにしてもよい。
【0353】
(2) 離散化の手法
本発明の実施形態で説明した式(22)〜(24)を離散化する手法は、様々な方法があり、式(25)〜(36)はその一例にすぎず、その他の方法で行なってもよい。
【0354】
(3) 64ビットレジスタ
本発明の第3の実施形態では、CPUが64ビット長のレジスタを備える場合を例にして説明したが、これは一例であって、Kビット長(Kは自然数)のレジスタ、たとえば16ビット長、32ビット長、128ビット長などのレジスタであってもよい。
【0355】
(4) LU分解
本発明の第1の実施形態では、第二種積分方程式を離散化して得られる連立方程式を構成する係数行列AをLU分解したが、これに限定するものではなく、下三角行列と上三角行列の積に分解できるものなら、どのような分解方法であってもよい。
【0356】
今回開示された実施の形態はすべての点で例示であって制限的なものではないと考えられるべきである。本発明の範囲は上記した説明ではなくて特許請求の範囲によって示され、特許請求の範囲と均等の意味および範囲内でのすべての変更が含まれることが意図される。
【図面の簡単な説明】
【0357】
【図1】第1の実施形態の逆ラプラス変換装置の構成を表わす図である。
【図2】Hテーブルを表わす図である。
【図3】第1の実施形態におけるHテーブルを作成する手順を表わすフローチャートである。
【図4】第1の実施形態におけるHテーブルを用いて逆ラプラス変換の数値解を算出する手順を表わすフローチャートである。
【図5】第1の実施形態の変形例の逆ラプラス変換装置の構成を表わす図である。
【図6】第1の実施形態の変形例におけるHテーブルを作成する手順を表わすフローチャートである。
【図7】第2の実施形態におけるHテーブルを作成する手順を表わすフローチャートである。
【図8】第3の実施形態のCPUの具体的な構成を表わす図である。
【図9】フラグレジスタの具体的な構成を表わす図である。
【図10】(a)、(b)、(c)は、図8の逆ラプラス変換装置において演算の対象となる数値の表現を説明するための図である。
【図11】フラグレジスタを備えるCPUの一例であるAMD64の加算処理の内容を表わす図である。
【図12】Alphaの加算処理の内容を表わす図である。
【図13】フラグレジスタを備えるCPUの一例であるAMD64の減算処理の内容を表わす図である。
【図14】Alphaの減算処理の内容を表わす図である。
【図15】本発明の実施形態で用いる乗算処理の内容を説明するための図である。
【図16】本発明の実施形態における乗算処理の内容を表わす図である。
【図17】本発明の実施形態で用いる乗算処理の内容を説明するための図である。
【図18】フラグレジスタを備えるCPUの一例であるAMD64のmul_addの内容を表わす図である。
【図19】Alphaのmul_addの内容を表わす図である。
【図20】ニュートン法を用いた除算処理の内容を説明するための図である。
【図21】図20中のサブルーチンcmpの内容を表わす図である。
【図22】(a)は、正則化パラメータa=10-1、10-4、10-8、10-12のときの、計算結果fa,s(n)(t)を示す図である。(b)は、正則化パラメータa=10-100、10-400のときの、計算結果fa,s(n)(t)を示す図である。
【図23】(a)は、s=0.1での計算結果fa,s(n)(t)を示す図であり、(b)は、s=0.01での計算結果fa,s(n)(t)を示す図である。
【図24】Hテーブルを作成する逆ラプラス変換用テーブル作成装置の構成を表わす図である。
【図25】逆ラプラス変換の数値解を算出する逆ラプラス変換の数値解算出装置の構成を表わす図である。
【符号の説明】
【0358】
1,81 逆ラプラス変換装置、61 逆ラプラス変換用テーブル作成装置、71 逆ラプラス変換数値解算出装置、2,62,72 入力部、3,63 CPU、4,84 テーブル作成部、5 逆変換部、6,64,74 メインメモリ、7,65,75 プログラム記憶部、8,66,76 変数記憶部、9 三角行列記憶部、12 Hテーブル記憶部、13 出力部、14 表示装置、67,77,99 リムーバブルメディア、53 制御部、54 演算器、56 汎用レジスタ群、57 フラグレジスタ、89 逆行列記憶部。
【特許請求の範囲】
【請求項1】
与えられた正則化パラメータと軟化子パラメータのもとで、ラプラス変換像に対する原像の数値解を算出するために、コンピュータに、
原点で零であり、かつ絶対連続な函数からなる重み付き再生核ヒルベルト空間上で、重み付き二乗可積分空間を観測空間としたチホノフ正則化法により導かれる第二種積分方程式を離散化して得られる連立方程式の解を求め、前記連立方程式の解に基づく前記第二種積分方程式の数値解を含む情報を記述したテーブルを作成するステップと、
前記作成したテーブルを記憶部に保存するステップと、
前記記憶部のテーブルを参照して、前記第二種積分方程式の数値解と軟化子函数を乗じたラプラス変換像との重み付き二乗可積分空間での内積を数値計算で求めることにより、原像の数値解を算出するステップと、
前記原像の数値解を出力するステップとを実行させる逆ラプラス変換プログラム。
【請求項2】
正則化パラメータをaとし、
軟化子パラメータsとし、
任意の関数gのラプラス変換像をLgとし、
数値解を算出する原像に対するラプラス変換像をF(p)とし、
軟化子函数をRs(p)としたときに、
原点で零であり、かつ絶対連続な函数fからなる重み付き再生核ヒルベルト空間Hkのノルムは、式(A1)で表わされ、
【数1】
前記重み付き再生核ヒルベルト空間Hkでの再生核K(t,t′)が式(A2)で表わされ、
【数2】
前記重み付き二乗可積分空間は、L2(R+,u(p)dp)で表わされ、
ρ(t)とu(p)の間には、式(A3)の条件が成立し、
【数3】
前記第二種積分方程式は、式(A4)で表わされ、
【数4】
前記内積は、式(A5)で表わされる、
【数5】
請求項1記載の逆ラプラス変換プログラム。
【請求項3】
ρ(t)は、式(A6)で表わされ、u(p)は、式(A7)で表わされ、Rs(p)は、式(A8)で表わされる、
【数6】
請求項2記載の逆ラプラス変換プログラム。
【請求項4】
前記コンピュータは、
Kビット長(Kは自然数)の汎用レジスタと、
直前の演算により繰上りまたは繰下がりが生じたか否かを示すフラグを格納するフラグレジスタと、
前記フラグを参照して、直前の演算により繰上りまたは繰下がりが生じたか否かに応じた演算を行なう演算器とを備え、
前記テーブルを作成するステップおよび前記原像の数値解を算出するステップは、前記数値計算の対象となる変数について、その仮数部をK×Nビット(Nは2以上の自然数)の多倍長で表わし、前記変数の仮数部をKビットごとに区切って、区切られた下位側から順次、前記汎用レジスタを用いて前記演算器に演算を行なわせるステップを含む、請求項3記載の逆ラプラス変換プログラム。
【請求項5】
前記連立方程式は、Ax=bの形で表わされ、Aは各要素がtに依存しない係数行列であり、xは連立方程式の解を要素とするベクトルであり、bは各要素がtに依存するベクトルであり、
前記テーブルを作成するステップは、
前記行列Aの逆行列A-1を算出して、前記逆行列A-1の要素を前記記憶部に記憶するステップと、
計算範囲内のすべてのtについて、前記記憶部に記憶されている前記逆行列A-1の要素の値に基づいて、前記連立方程式の解を算出するステップとを含む、請求項3記載の逆ラプラス変換プログラム。
【請求項6】
前記連立方程式は、Ax=bの形で表わされ、Aは各要素がtに依存しない係数行列であり、xは連立方程式の解を要素とするベクトルであり、bは各要素がtに依存するベクトルであり、
前記テーブルを作成するステップは、
前記行列Aを上三角行列と下三角行列の積に分解して、前記上三角行列と前記下三角行列の要素を前記記憶部に記憶するステップと、
計算範囲内のすべてのtについて、前記記憶部に記憶されている前記上三角行列と前記下三角行列の要素の値に基づいて、前記連立方程式の解を算出するステップとを含む、請求項3記載の逆ラプラス変換プログラム。
【請求項7】
前記テーブルを作成するステップは、
前記正則化パラメータaおよび前記軟化子パラメータsと、その他のパラメータn、U、Lとから、式(A9)〜(A14)に従って、それぞれ変数h、xj、pj、qj、aij、H(pj, t)を計算するステップと、
【数7】
前記連立方程式は、式(A15)で表わされ、
式(A16)で表わされる行列Aをコレスキー分解し、式(A15)の解を算出するステップと、
【数8】
式(A17)に従って、式(A15)の解から前記第二種積分方程式の数値解{Hin,a,t:i=0,1,2,・・・,n}を算出し、前記第二種積分方程式の数値解を含む情報を記述したテーブルを作成するステップとを含み、
【数9】
前記原像の数値解を算出するステップは、前記テーブルを参照して、式(A18)に従って、ラプラス変換像F(pj)から、前記原像の数値解fa,s(n)(t)を算出するステップを含む、
【数10】
請求項3記載の逆ラプラス変換プログラム。
【請求項8】
与えられた正則化パラメータと軟化子パラメータのもとで、ラプラス変換像に対する原像の数値解を算出するためのテーブルを作成するために、コンピュータに、
原点で零であり、かつ絶対連続な函数からなる重み付き再生核ヒルベルト空間上で、重み付き二乗可積分空間を観測空間としたチホノフ正則化法により導かれる第二種積分方程式を離散化して得られる連立方程式の解を求め、前記連立方程式の解に基づく前記第二種積分方程式の数値解を含む情報を記述したテーブルを作成するステップと、
前記作成したテーブルを記憶部に保存するステップとを実行させるテーブル作成プログラム。
【請求項9】
正則化パラメータをaとし、
軟化子パラメータsとし、
任意の関数gのラプラス変換像をLgとし、
原点で零であり、かつ絶対連続な函数fからなる重み付き再生核ヒルベルト空間Hkのノルムは、式(A19)で表わされ、
【数11】
前記重み付き再生核ヒルベルト空間Hkでの再生核K(t,t′)が式(A20)で表わされ、
【数12】
前記重み付き二乗可積分空間は、L2(R+,u(p)dp)で表わされ、
ρ(t)とu(p)の間には、式(A21)の条件が成立し、
【数13】
前記第二種積分方程式は、式(A22)で表わされる、
【数14】
請求項8記載のテーブル作成プログラム。
【請求項10】
ρ(t)は、式(A23)で表わされ、u(p)は、式(A24)で表わされる、
【数15】
請求項9記載のテーブル作成プログラム。
【請求項11】
前記コンピュータは、
Kビット長(Kは自然数)の汎用レジスタと、
直前の演算により繰上りまたは繰下がりが生じたか否かを示すフラグを格納するフラグレジスタと、
前記フラグを参照して、直前の演算により繰上りまたは繰下がりが生じたか否かに応じた演算を行なう演算器とを備え、
前記テーブルを作成するステップは、前記数値計算の対象となる変数について、その仮数部をK×Nビット(Nは2以上の自然数)の多倍長で表わし、前記変数の仮数部をKビットごとに区切って、区切られた下位側から順次、前記汎用レジスタを用いて前記演算器に演算を行なわせるステップを含む、請求項10記載のテーブル作成プログラム。
【請求項12】
前記連立方程式は、Ax=bの形で表わされ、Aは各要素がtに依存しない係数行列であり、xは連立方程式の解を要素とするベクトルであり、bは各要素がtに依存するベクトルであり、
前記テーブルを作成するステップは、
前記行列Aの逆行列A-1を算出して、前記逆行列A-1の要素を前記記憶部に記憶するステップと、
計算範囲内のすべてのtについて、前記記憶部に記憶されている前記逆行列A-1の要素の値に基づいて、前記連立方程式の解を算出するステップとを含む、請求項10記載のテーブル作成プログラム。
【請求項13】
前記連立方程式は、Ax=bの形で表わされ、Aは各要素がtに依存しない係数行列であり、xは連立方程式の解を要素とするベクトルであり、bは各要素がtに依存するベクトルであり、
前記テーブルを作成するステップは、
前記行列Aを上三角行列と下三角行列の積に分解して、前記上三角行列と前記下三角行列の要素を前記記憶部に記憶するステップと、
計算範囲内のすべてのtについて、前記記憶部に記憶されている前記上三角行列と前記下三角行列の要素の値に基づいて、前記連立方程式の解を算出するステップとを含む、請求項10記載のテーブル作成プログラム。
【請求項14】
前記テーブルを作成するステップは、
前記正則化パラメータaおよび前記軟化子パラメータsと、その他のパラメータn、U、Lとから、式(A25)〜(A30)に従って、それぞれ変数h、xj、pj、qj、aij、H(pj, t)を計算するステップと、
【数16】
前記連立方程式は、式(A31)で表わされ、
式(A32)で表わされる行列Aをコレスキー分解し、式(A31)の解を算出するステップと、
【数17】
式(A33)に従って、式(A31)の解から前記第二種積分方程式の数値解{Hin,a,t:i=0,1,2,・・・,n}を算出し、前記第二種積分方程式の数値解を含む情報を記述したテーブルを作成するステップとを含む、
【数18】
請求項10記載のテーブル作成プログラム。
【請求項15】
請求項8に記載の逆ラプラス変換のためのテーブル作成プログラムによって作成されたテーブルを用いて、ラプラス変換像に対する原像の数値解を作成するために、コンピュータに、
請求項8に記載の逆ラプラス変換のためのテーブル作成プログラムで作成されたテーブルが入力され、入力されたテーブルを記憶部に保存するステップと、
前記記憶部のテーブルを参照して、前記第二種積分方程式の数値解と軟化子函数を乗じたラプラス変換像との重み付き二乗可積分空間での内積を数値計算で求めることにより、原像の数値解を算出するステップと、
前記原像の数値解を出力するステップとを実行させる逆ラプラス変換の数値解算出プログラム。
【請求項16】
請求項9に記載の逆ラプラス変換のためのテーブル作成プログラムによって作成されたテーブルを用いて、ラプラス変換像に対する原像の数値解を作成するために、コンピュータに、
請求項9に記載の逆ラプラス変換のためのテーブル作成プログラムで作成されたテーブルが入力され、入力されたテーブルを記憶部に保存するステップと、
前記記憶部のテーブルを参照して、前記第二種積分方程式の数値解と軟化子函数Rs(p)を乗じたラプラス変換像F(p)との重み付き二乗可積分空間での内積を数値計算で求めることにより、原像の数値解を算出するステップと、前記内積は、式(A34)で表わされ、
【数19】
前記原像の数値解を出力するステップとを実行させる逆ラプラス変換の数値解算出プログラム。
【請求項17】
請求項10、12、13に記載の逆ラプラス変換のためのテーブル作成プログラムによって作成されたテーブルを用いて、ラプラス変換像に対する原像の数値解を作成するために、コンピュータに、
請求項10、12、13に記載の逆ラプラス変換のためのテーブル作成プログラムで作成されたテーブルが入力され、入力されたテーブルを記憶部に保存するステップと、
前記記憶部のテーブルを参照して、前記第二種積分方程式の数値解と軟化子函数Rs(p)を乗じたラプラス変換像F(p)との重み付き二乗可積分空間での内積を数値計算で求めることにより、原像の数値解を算出するステップと、前記軟化子函数Rs(p)は、式(A35)で表わされ、
【数20】
前記内積は、式(A36)で表わされ、
【数21】
前記原像の数値解を出力するステップとを実行させる逆ラプラス変換の数値解算出プログラム。
【請求項18】
請求項11記載の逆ラプラス変換のためのテーブル作成プログラムによって作成されたテーブルを用いて、逆ラプラス変換の数値解を作成するために、
Kビット長(Kは自然数)の汎用レジスタと、直前の演算により繰上りまたは繰下がりが生じたか否かを示すフラグを格納するフラグレジスタと、前記フラグを参照して、直前の演算により繰上りまたは繰下がりが生じたか否かに応じた演算を行なう演算器とを備えたコンピュータに、
請求項11記載の逆ラプラス変換のためのテーブル作成プログラムで作成されたテーブルが入力され、入力されたテーブルを記憶部に保存するステップと、
前記記憶部のテーブルを参照して、前記第二種積分方程式の数値解と軟化子函数Rs(p)を乗じたラプラス変換像F(p)との重み付き二乗可積分空間での内積を数値計算で求めることにより、原像の数値解を算出するステップと、前記軟化子函数Rs(p)は、式(A37)で表わされ、
【数22】
前記内積は、式(A38)で表わされ、
【数23】
前記原像の数値解を出力するステップとを実行させ、
前記原像の数値解を算出するステップは、前記数値計算の対象となる変数について、その仮数部をK×Nビット(Nは2以上の自然数)の多倍長で表わし、前記変数の仮数部をKビットごとに区切って、区切られた下位側から順次、前記汎用レジスタを用いて前記演算器に演算を行なわせるステップを含む、逆ラプラス変換の数値解算出プログラム。
【請求項19】
請求項14記載の逆ラプラス変換のためのテーブル作成プログラムによって作成されたテーブルを用いて、逆ラプラス変換の数値解を作成するために、コンピュータに、
請求項14記載の逆ラプラス変換のためのテーブル作成プログラムで作成されたテーブルが入力され、入力されたテーブルを記憶部に保存するステップと、
前記記憶部のテーブルを参照して、式(A39)に従って、前記第二種積分方程式の数値解とラプラス変換像F(pj)から前記原像の数値解fa,s(n)(t)を算出するステップと、
【数24】
前記原像の数値解を出力するステップとを実行させる逆ラプラス変換の数値解算出プログラム。
【請求項20】
与えられた正則化パラメータと軟化子パラメータのもとで、ラプラス変換像に対する原像の数値解を算出する逆ラプラス変換装置であって、
原点で零であり、かつ絶対連続な函数からなる重み付き再生核ヒルベルト空間上で、重み付き二乗可積分空間を観測空間としたチホノフ正則化法により導かれる第二種積分方程式を離散化して得られる連立方程式の解を求め、前記連立方程式の解に基づく前記第二種積分方程式の数値解を含む情報を記述したテーブルを作成するテーブル作成部と、
前記作成したテーブルを記憶する記憶部と、
前記記憶部のテーブルを参照して、前記第二種積分方程式の数値解と軟化子函数を乗じたラプラス変換像との重み付き二乗可積分空間での内積を数値計算で求めることにより、原像の数値解を算出する逆変換部と、
前記原像の数値解を出力する出力部とを備えた逆ラプラス変換装置。
【請求項21】
正則化パラメータをaとし、
軟化子パラメータsとし、
任意の関数gのラプラス変換像をLgとし、
数値解を算出する原像に対するラプラス変換像をF(p)とし、
軟化子函数をRs(p)としたときに、
原点で零であり、かつ絶対連続な函数fからなる重み付き再生核ヒルベルト空間Hkのノルムは、式(A40)で表わされ、
【数25】
前記重み付き再生核ヒルベルト空間Hkでの再生核K(t,t′)が式(A41)で表わされ、
【数26】
前記重み付き二乗可積分空間は、L2(R+,u(p)dp)で表わされ、
ρ(t)とu(p)の間には、式(A42)の条件が成立し、
【数27】
前記第二種積分方程式は、式(A43)で表わされ、
【数28】
前記内積は、式(A44)で表わされる、
【数29】
請求項20記載の逆ラプラス変換装置。
【請求項22】
ρ(t)は、式(A45)で表わされ、u(p)は、式(A46)で表わされ、Rs(p)は、式(A47)で表わされる、
【数30】
請求項21記載の逆ラプラス変換装置。
【請求項23】
前記テーブル作成部および前記逆変換部とは、共通して、
Kビット長(Kは自然数)の汎用レジスタと、
直前の演算により繰上りまたは繰下がりが生じたか否かを示すフラグを格納するフラグレジスタと、
前記フラグを参照して、直前の演算により繰上りまたは繰下がりが生じたか否かに応じた演算を行なう演算器とを備え、
前記演算器は、前記数値計算の対象となる変数について、その仮数部をK×Nビット(Nは2以上の自然数)の多倍長で表わし、前記変数の仮数部をKビットごとに区切って、区切られた下位側から順次、前記汎用レジスタを用いて演算を行なう、請求項22記載の逆ラプラス変換装置。
【請求項24】
前記連立方程式は、Ax=bの形で表わされ、Aは各要素がtに依存しない係数行列であり、xは連立方程式の解を要素とするベクトルであり、bは各要素がtに依存するベクトルであり、
前記テーブル作成部は、前記行列Aの逆行列A-1を算出して、前記逆行列A-1の要素を前記記憶部に記憶し、
計算範囲内のすべてのtについて、前記記憶部に記憶されている前記逆行列A-1の要素の値に基づいて、前記連立方程式の解を算出する、請求項22記載の逆ラプラス変換装置。
【請求項25】
前記連立方程式は、Ax=bの形で表わされ、Aは各要素がtに依存しない係数行列であり、xは連立方程式の解を要素とするベクトルであり、bは各要素がtに依存するベクトルであり、
前記テーブル作成部は、前記行列Aを上三角行列と下三角行列の積に分解して、前記上三角行列と前記下三角行列の要素を前記記憶部に記憶し、
計算範囲内のすべてのtについて、前記記憶部に記憶されている前記上三角行列と前記下三角行列の要素の値に基づいて、前記連立方程式の解を算出する、請求項22記載の逆ラプラス変換装置。
【請求項26】
前記テーブル作成部は、
前記正則化パラメータaおよび前記軟化子パラメータsと、その他のパラメータn、U、Lとから、式(A48)〜(A53)に従って、それぞれ変数h、xj、pj、qj、aij、H(pj, t)を計算し、
【数31】
前記連立方程式は、式(A54)で表わされ、前記テーブル作成部は、式(A55)で表わされる行列Aをコレスキー分解し、式(A54)の解を算出し、
【数32】
前記テーブル作成部は、さらに式(A56)に従って、式(A54)の解から前記第二種積分方程式の数値解{Hin,a,t:i=0,1,2,・・・,n}を算出し、前記第二種積分方程式の数値解を含む情報を記述したテーブルを作成し、
【数33】
前記逆変換部は、前記テーブルを参照して、式(A57)に従って、ラプラス変換像F(pj)から前記原像の数値解fa,s(n)(t)を算出する、
【数34】
請求項22記載の逆ラプラス変換装置。
【請求項1】
与えられた正則化パラメータと軟化子パラメータのもとで、ラプラス変換像に対する原像の数値解を算出するために、コンピュータに、
原点で零であり、かつ絶対連続な函数からなる重み付き再生核ヒルベルト空間上で、重み付き二乗可積分空間を観測空間としたチホノフ正則化法により導かれる第二種積分方程式を離散化して得られる連立方程式の解を求め、前記連立方程式の解に基づく前記第二種積分方程式の数値解を含む情報を記述したテーブルを作成するステップと、
前記作成したテーブルを記憶部に保存するステップと、
前記記憶部のテーブルを参照して、前記第二種積分方程式の数値解と軟化子函数を乗じたラプラス変換像との重み付き二乗可積分空間での内積を数値計算で求めることにより、原像の数値解を算出するステップと、
前記原像の数値解を出力するステップとを実行させる逆ラプラス変換プログラム。
【請求項2】
正則化パラメータをaとし、
軟化子パラメータsとし、
任意の関数gのラプラス変換像をLgとし、
数値解を算出する原像に対するラプラス変換像をF(p)とし、
軟化子函数をRs(p)としたときに、
原点で零であり、かつ絶対連続な函数fからなる重み付き再生核ヒルベルト空間Hkのノルムは、式(A1)で表わされ、
【数1】
前記重み付き再生核ヒルベルト空間Hkでの再生核K(t,t′)が式(A2)で表わされ、
【数2】
前記重み付き二乗可積分空間は、L2(R+,u(p)dp)で表わされ、
ρ(t)とu(p)の間には、式(A3)の条件が成立し、
【数3】
前記第二種積分方程式は、式(A4)で表わされ、
【数4】
前記内積は、式(A5)で表わされる、
【数5】
請求項1記載の逆ラプラス変換プログラム。
【請求項3】
ρ(t)は、式(A6)で表わされ、u(p)は、式(A7)で表わされ、Rs(p)は、式(A8)で表わされる、
【数6】
請求項2記載の逆ラプラス変換プログラム。
【請求項4】
前記コンピュータは、
Kビット長(Kは自然数)の汎用レジスタと、
直前の演算により繰上りまたは繰下がりが生じたか否かを示すフラグを格納するフラグレジスタと、
前記フラグを参照して、直前の演算により繰上りまたは繰下がりが生じたか否かに応じた演算を行なう演算器とを備え、
前記テーブルを作成するステップおよび前記原像の数値解を算出するステップは、前記数値計算の対象となる変数について、その仮数部をK×Nビット(Nは2以上の自然数)の多倍長で表わし、前記変数の仮数部をKビットごとに区切って、区切られた下位側から順次、前記汎用レジスタを用いて前記演算器に演算を行なわせるステップを含む、請求項3記載の逆ラプラス変換プログラム。
【請求項5】
前記連立方程式は、Ax=bの形で表わされ、Aは各要素がtに依存しない係数行列であり、xは連立方程式の解を要素とするベクトルであり、bは各要素がtに依存するベクトルであり、
前記テーブルを作成するステップは、
前記行列Aの逆行列A-1を算出して、前記逆行列A-1の要素を前記記憶部に記憶するステップと、
計算範囲内のすべてのtについて、前記記憶部に記憶されている前記逆行列A-1の要素の値に基づいて、前記連立方程式の解を算出するステップとを含む、請求項3記載の逆ラプラス変換プログラム。
【請求項6】
前記連立方程式は、Ax=bの形で表わされ、Aは各要素がtに依存しない係数行列であり、xは連立方程式の解を要素とするベクトルであり、bは各要素がtに依存するベクトルであり、
前記テーブルを作成するステップは、
前記行列Aを上三角行列と下三角行列の積に分解して、前記上三角行列と前記下三角行列の要素を前記記憶部に記憶するステップと、
計算範囲内のすべてのtについて、前記記憶部に記憶されている前記上三角行列と前記下三角行列の要素の値に基づいて、前記連立方程式の解を算出するステップとを含む、請求項3記載の逆ラプラス変換プログラム。
【請求項7】
前記テーブルを作成するステップは、
前記正則化パラメータaおよび前記軟化子パラメータsと、その他のパラメータn、U、Lとから、式(A9)〜(A14)に従って、それぞれ変数h、xj、pj、qj、aij、H(pj, t)を計算するステップと、
【数7】
前記連立方程式は、式(A15)で表わされ、
式(A16)で表わされる行列Aをコレスキー分解し、式(A15)の解を算出するステップと、
【数8】
式(A17)に従って、式(A15)の解から前記第二種積分方程式の数値解{Hin,a,t:i=0,1,2,・・・,n}を算出し、前記第二種積分方程式の数値解を含む情報を記述したテーブルを作成するステップとを含み、
【数9】
前記原像の数値解を算出するステップは、前記テーブルを参照して、式(A18)に従って、ラプラス変換像F(pj)から、前記原像の数値解fa,s(n)(t)を算出するステップを含む、
【数10】
請求項3記載の逆ラプラス変換プログラム。
【請求項8】
与えられた正則化パラメータと軟化子パラメータのもとで、ラプラス変換像に対する原像の数値解を算出するためのテーブルを作成するために、コンピュータに、
原点で零であり、かつ絶対連続な函数からなる重み付き再生核ヒルベルト空間上で、重み付き二乗可積分空間を観測空間としたチホノフ正則化法により導かれる第二種積分方程式を離散化して得られる連立方程式の解を求め、前記連立方程式の解に基づく前記第二種積分方程式の数値解を含む情報を記述したテーブルを作成するステップと、
前記作成したテーブルを記憶部に保存するステップとを実行させるテーブル作成プログラム。
【請求項9】
正則化パラメータをaとし、
軟化子パラメータsとし、
任意の関数gのラプラス変換像をLgとし、
原点で零であり、かつ絶対連続な函数fからなる重み付き再生核ヒルベルト空間Hkのノルムは、式(A19)で表わされ、
【数11】
前記重み付き再生核ヒルベルト空間Hkでの再生核K(t,t′)が式(A20)で表わされ、
【数12】
前記重み付き二乗可積分空間は、L2(R+,u(p)dp)で表わされ、
ρ(t)とu(p)の間には、式(A21)の条件が成立し、
【数13】
前記第二種積分方程式は、式(A22)で表わされる、
【数14】
請求項8記載のテーブル作成プログラム。
【請求項10】
ρ(t)は、式(A23)で表わされ、u(p)は、式(A24)で表わされる、
【数15】
請求項9記載のテーブル作成プログラム。
【請求項11】
前記コンピュータは、
Kビット長(Kは自然数)の汎用レジスタと、
直前の演算により繰上りまたは繰下がりが生じたか否かを示すフラグを格納するフラグレジスタと、
前記フラグを参照して、直前の演算により繰上りまたは繰下がりが生じたか否かに応じた演算を行なう演算器とを備え、
前記テーブルを作成するステップは、前記数値計算の対象となる変数について、その仮数部をK×Nビット(Nは2以上の自然数)の多倍長で表わし、前記変数の仮数部をKビットごとに区切って、区切られた下位側から順次、前記汎用レジスタを用いて前記演算器に演算を行なわせるステップを含む、請求項10記載のテーブル作成プログラム。
【請求項12】
前記連立方程式は、Ax=bの形で表わされ、Aは各要素がtに依存しない係数行列であり、xは連立方程式の解を要素とするベクトルであり、bは各要素がtに依存するベクトルであり、
前記テーブルを作成するステップは、
前記行列Aの逆行列A-1を算出して、前記逆行列A-1の要素を前記記憶部に記憶するステップと、
計算範囲内のすべてのtについて、前記記憶部に記憶されている前記逆行列A-1の要素の値に基づいて、前記連立方程式の解を算出するステップとを含む、請求項10記載のテーブル作成プログラム。
【請求項13】
前記連立方程式は、Ax=bの形で表わされ、Aは各要素がtに依存しない係数行列であり、xは連立方程式の解を要素とするベクトルであり、bは各要素がtに依存するベクトルであり、
前記テーブルを作成するステップは、
前記行列Aを上三角行列と下三角行列の積に分解して、前記上三角行列と前記下三角行列の要素を前記記憶部に記憶するステップと、
計算範囲内のすべてのtについて、前記記憶部に記憶されている前記上三角行列と前記下三角行列の要素の値に基づいて、前記連立方程式の解を算出するステップとを含む、請求項10記載のテーブル作成プログラム。
【請求項14】
前記テーブルを作成するステップは、
前記正則化パラメータaおよび前記軟化子パラメータsと、その他のパラメータn、U、Lとから、式(A25)〜(A30)に従って、それぞれ変数h、xj、pj、qj、aij、H(pj, t)を計算するステップと、
【数16】
前記連立方程式は、式(A31)で表わされ、
式(A32)で表わされる行列Aをコレスキー分解し、式(A31)の解を算出するステップと、
【数17】
式(A33)に従って、式(A31)の解から前記第二種積分方程式の数値解{Hin,a,t:i=0,1,2,・・・,n}を算出し、前記第二種積分方程式の数値解を含む情報を記述したテーブルを作成するステップとを含む、
【数18】
請求項10記載のテーブル作成プログラム。
【請求項15】
請求項8に記載の逆ラプラス変換のためのテーブル作成プログラムによって作成されたテーブルを用いて、ラプラス変換像に対する原像の数値解を作成するために、コンピュータに、
請求項8に記載の逆ラプラス変換のためのテーブル作成プログラムで作成されたテーブルが入力され、入力されたテーブルを記憶部に保存するステップと、
前記記憶部のテーブルを参照して、前記第二種積分方程式の数値解と軟化子函数を乗じたラプラス変換像との重み付き二乗可積分空間での内積を数値計算で求めることにより、原像の数値解を算出するステップと、
前記原像の数値解を出力するステップとを実行させる逆ラプラス変換の数値解算出プログラム。
【請求項16】
請求項9に記載の逆ラプラス変換のためのテーブル作成プログラムによって作成されたテーブルを用いて、ラプラス変換像に対する原像の数値解を作成するために、コンピュータに、
請求項9に記載の逆ラプラス変換のためのテーブル作成プログラムで作成されたテーブルが入力され、入力されたテーブルを記憶部に保存するステップと、
前記記憶部のテーブルを参照して、前記第二種積分方程式の数値解と軟化子函数Rs(p)を乗じたラプラス変換像F(p)との重み付き二乗可積分空間での内積を数値計算で求めることにより、原像の数値解を算出するステップと、前記内積は、式(A34)で表わされ、
【数19】
前記原像の数値解を出力するステップとを実行させる逆ラプラス変換の数値解算出プログラム。
【請求項17】
請求項10、12、13に記載の逆ラプラス変換のためのテーブル作成プログラムによって作成されたテーブルを用いて、ラプラス変換像に対する原像の数値解を作成するために、コンピュータに、
請求項10、12、13に記載の逆ラプラス変換のためのテーブル作成プログラムで作成されたテーブルが入力され、入力されたテーブルを記憶部に保存するステップと、
前記記憶部のテーブルを参照して、前記第二種積分方程式の数値解と軟化子函数Rs(p)を乗じたラプラス変換像F(p)との重み付き二乗可積分空間での内積を数値計算で求めることにより、原像の数値解を算出するステップと、前記軟化子函数Rs(p)は、式(A35)で表わされ、
【数20】
前記内積は、式(A36)で表わされ、
【数21】
前記原像の数値解を出力するステップとを実行させる逆ラプラス変換の数値解算出プログラム。
【請求項18】
請求項11記載の逆ラプラス変換のためのテーブル作成プログラムによって作成されたテーブルを用いて、逆ラプラス変換の数値解を作成するために、
Kビット長(Kは自然数)の汎用レジスタと、直前の演算により繰上りまたは繰下がりが生じたか否かを示すフラグを格納するフラグレジスタと、前記フラグを参照して、直前の演算により繰上りまたは繰下がりが生じたか否かに応じた演算を行なう演算器とを備えたコンピュータに、
請求項11記載の逆ラプラス変換のためのテーブル作成プログラムで作成されたテーブルが入力され、入力されたテーブルを記憶部に保存するステップと、
前記記憶部のテーブルを参照して、前記第二種積分方程式の数値解と軟化子函数Rs(p)を乗じたラプラス変換像F(p)との重み付き二乗可積分空間での内積を数値計算で求めることにより、原像の数値解を算出するステップと、前記軟化子函数Rs(p)は、式(A37)で表わされ、
【数22】
前記内積は、式(A38)で表わされ、
【数23】
前記原像の数値解を出力するステップとを実行させ、
前記原像の数値解を算出するステップは、前記数値計算の対象となる変数について、その仮数部をK×Nビット(Nは2以上の自然数)の多倍長で表わし、前記変数の仮数部をKビットごとに区切って、区切られた下位側から順次、前記汎用レジスタを用いて前記演算器に演算を行なわせるステップを含む、逆ラプラス変換の数値解算出プログラム。
【請求項19】
請求項14記載の逆ラプラス変換のためのテーブル作成プログラムによって作成されたテーブルを用いて、逆ラプラス変換の数値解を作成するために、コンピュータに、
請求項14記載の逆ラプラス変換のためのテーブル作成プログラムで作成されたテーブルが入力され、入力されたテーブルを記憶部に保存するステップと、
前記記憶部のテーブルを参照して、式(A39)に従って、前記第二種積分方程式の数値解とラプラス変換像F(pj)から前記原像の数値解fa,s(n)(t)を算出するステップと、
【数24】
前記原像の数値解を出力するステップとを実行させる逆ラプラス変換の数値解算出プログラム。
【請求項20】
与えられた正則化パラメータと軟化子パラメータのもとで、ラプラス変換像に対する原像の数値解を算出する逆ラプラス変換装置であって、
原点で零であり、かつ絶対連続な函数からなる重み付き再生核ヒルベルト空間上で、重み付き二乗可積分空間を観測空間としたチホノフ正則化法により導かれる第二種積分方程式を離散化して得られる連立方程式の解を求め、前記連立方程式の解に基づく前記第二種積分方程式の数値解を含む情報を記述したテーブルを作成するテーブル作成部と、
前記作成したテーブルを記憶する記憶部と、
前記記憶部のテーブルを参照して、前記第二種積分方程式の数値解と軟化子函数を乗じたラプラス変換像との重み付き二乗可積分空間での内積を数値計算で求めることにより、原像の数値解を算出する逆変換部と、
前記原像の数値解を出力する出力部とを備えた逆ラプラス変換装置。
【請求項21】
正則化パラメータをaとし、
軟化子パラメータsとし、
任意の関数gのラプラス変換像をLgとし、
数値解を算出する原像に対するラプラス変換像をF(p)とし、
軟化子函数をRs(p)としたときに、
原点で零であり、かつ絶対連続な函数fからなる重み付き再生核ヒルベルト空間Hkのノルムは、式(A40)で表わされ、
【数25】
前記重み付き再生核ヒルベルト空間Hkでの再生核K(t,t′)が式(A41)で表わされ、
【数26】
前記重み付き二乗可積分空間は、L2(R+,u(p)dp)で表わされ、
ρ(t)とu(p)の間には、式(A42)の条件が成立し、
【数27】
前記第二種積分方程式は、式(A43)で表わされ、
【数28】
前記内積は、式(A44)で表わされる、
【数29】
請求項20記載の逆ラプラス変換装置。
【請求項22】
ρ(t)は、式(A45)で表わされ、u(p)は、式(A46)で表わされ、Rs(p)は、式(A47)で表わされる、
【数30】
請求項21記載の逆ラプラス変換装置。
【請求項23】
前記テーブル作成部および前記逆変換部とは、共通して、
Kビット長(Kは自然数)の汎用レジスタと、
直前の演算により繰上りまたは繰下がりが生じたか否かを示すフラグを格納するフラグレジスタと、
前記フラグを参照して、直前の演算により繰上りまたは繰下がりが生じたか否かに応じた演算を行なう演算器とを備え、
前記演算器は、前記数値計算の対象となる変数について、その仮数部をK×Nビット(Nは2以上の自然数)の多倍長で表わし、前記変数の仮数部をKビットごとに区切って、区切られた下位側から順次、前記汎用レジスタを用いて演算を行なう、請求項22記載の逆ラプラス変換装置。
【請求項24】
前記連立方程式は、Ax=bの形で表わされ、Aは各要素がtに依存しない係数行列であり、xは連立方程式の解を要素とするベクトルであり、bは各要素がtに依存するベクトルであり、
前記テーブル作成部は、前記行列Aの逆行列A-1を算出して、前記逆行列A-1の要素を前記記憶部に記憶し、
計算範囲内のすべてのtについて、前記記憶部に記憶されている前記逆行列A-1の要素の値に基づいて、前記連立方程式の解を算出する、請求項22記載の逆ラプラス変換装置。
【請求項25】
前記連立方程式は、Ax=bの形で表わされ、Aは各要素がtに依存しない係数行列であり、xは連立方程式の解を要素とするベクトルであり、bは各要素がtに依存するベクトルであり、
前記テーブル作成部は、前記行列Aを上三角行列と下三角行列の積に分解して、前記上三角行列と前記下三角行列の要素を前記記憶部に記憶し、
計算範囲内のすべてのtについて、前記記憶部に記憶されている前記上三角行列と前記下三角行列の要素の値に基づいて、前記連立方程式の解を算出する、請求項22記載の逆ラプラス変換装置。
【請求項26】
前記テーブル作成部は、
前記正則化パラメータaおよび前記軟化子パラメータsと、その他のパラメータn、U、Lとから、式(A48)〜(A53)に従って、それぞれ変数h、xj、pj、qj、aij、H(pj, t)を計算し、
【数31】
前記連立方程式は、式(A54)で表わされ、前記テーブル作成部は、式(A55)で表わされる行列Aをコレスキー分解し、式(A54)の解を算出し、
【数32】
前記テーブル作成部は、さらに式(A56)に従って、式(A54)の解から前記第二種積分方程式の数値解{Hin,a,t:i=0,1,2,・・・,n}を算出し、前記第二種積分方程式の数値解を含む情報を記述したテーブルを作成し、
【数33】
前記逆変換部は、前記テーブルを参照して、式(A57)に従って、ラプラス変換像F(pj)から前記原像の数値解fa,s(n)(t)を算出する、
【数34】
請求項22記載の逆ラプラス変換装置。
【図1】
【図2】
【図3】
【図4】
【図5】
【図6】
【図7】
【図8】
【図9】
【図10】
【図11】
【図12】
【図13】
【図14】
【図15】
【図16】
【図17】
【図18】
【図19】
【図20】
【図21】
【図22】
【図23】
【図24】
【図25】
【図2】
【図3】
【図4】
【図5】
【図6】
【図7】
【図8】
【図9】
【図10】
【図11】
【図12】
【図13】
【図14】
【図15】
【図16】
【図17】
【図18】
【図19】
【図20】
【図21】
【図22】
【図23】
【図24】
【図25】
【公開番号】特開2009−64424(P2009−64424A)
【公開日】平成21年3月26日(2009.3.26)
【国際特許分類】
【出願番号】特願2008−203384(P2008−203384)
【出願日】平成20年8月6日(2008.8.6)
【出願人】(504132272)国立大学法人京都大学 (1,269)
【出願人】(504145364)国立大学法人群馬大学 (352)
【Fターム(参考)】
【公開日】平成21年3月26日(2009.3.26)
【国際特許分類】
【出願日】平成20年8月6日(2008.8.6)
【出願人】(504132272)国立大学法人京都大学 (1,269)
【出願人】(504145364)国立大学法人群馬大学 (352)
【Fターム(参考)】
[ Back to top ]