知財求人 - 知財ポータルサイト「IP Force」
▶ 大貫 進一郎の特許一覧 ▶ 地方独立行政法人東京都立産業技術研究センターの特許一覧
(19)【発行国】日本国特許庁(JP)
(12)【公報種別】特許公報(B2)
(11)【特許番号】
(24)【登録日】2022-09-29
(45)【発行日】2022-10-07
(54)【発明の名称】演算システム、演算装置及びプログラム
(51)【国際特許分類】
G06F 30/23 20200101AFI20220930BHJP
G06F 30/10 20200101ALI20220930BHJP
G16Z 99/00 20190101ALI20220930BHJP
【ＦＩ】
G06F30/23
G06F30/10
G16Z99/00
【請求項の数】
4
(21)【出願番号】P 2018012042
(22)【出願日】2018-01-26
(65)【公開番号】P2019128918
(43)【公開日】2019-08-01
【審査請求日】2021-01-15
(73)【特許権者】
【識別番号】521022808
【氏名又は名称】大貫 進一郎
(73)【特許権者】
【識別番号】506209422
【氏名又は名称】地方独立行政法人東京都立産業技術研究センター
(74)【代理人】
【識別番号】100161207
【弁理士】
【氏名又は名称】西澤 和純
(74)【代理人】
【識別番号】100175824
【弁理士】
【氏名又は名称】小林 淳一
(74)【代理人】
【識別番号】100126882
【弁理士】
【氏名又は名称】五十嵐 光永
(72)【発明者】
【氏名】大貫 進一郎
(72)【発明者】
【氏名】山口 隆志
(72)【発明者】
【氏名】呉 迪
(72)【発明者】
【氏名】大西 崚平
【審査官】合田 幸裕
(56)【参考文献】
【文献】特開２０１０－１８６３７２（ＪＰ，Ａ）
【文献】国際公開第２０１６／０２７４５２（ＷＯ，Ａ１）
【文献】岸本誠也, 大貫進一郎，並列計算による BIEM-FILT の高速化，電子情報通信学会２０１４年総合大会講演論文集（エレクトロニクス講演論文集１），日本，一般社団法人電子情報通信学会，2014年03月04日，第26頁
(58)【調査した分野】(Int.Cl.，ＤＢ名)
Ｇ０６Ｆ ３０／２３
Ｇ０６Ｆ ３０／１０
Ｇ１６Ｚ ９９／００
ＩＥＥＥ Ｘｐｌｏｒｅ
ＪＳＴＰｌｕｓ（ＪＤｒｅａｍＩＩＩ）
(57)【特許請求の範囲】
【請求項1】
時間変化する電磁界強度の波形についての複素周波数領域の関数に対する演算の時間軸における起点である初期時刻と、当該初期時刻からの演算時間幅とを、演算条件として取得する演算条件取得部と、前記演算条件取得部によって取得される前記初期時刻における前記波形についての演算結果である初期値を、前記関数に対する逆ラプラス変換法による離散的な時間応答解析によって算出するとともに、算出された前記初期値に基づく時間軸に沿った逐次演算によって、前記演算条件取得部によって取得される前記演算時間幅ぶんの前記波形についての演算結果を算出する演算部と、
前記演算時間幅毎に前記演算部によって算出された前記演算結果を演算管理装置に対して 出力する出力部と、
を備える演算装置と、
それぞれの前記演算装置に問い合わせることにより取得した前記演算装置ごとの前記演算部の演算性能に基づいて、前記初期時刻からの演算時間幅を前記演算装置ごとに算出する演算時間幅算出部と、
複数の前記演算装置に対して、前記演算装置ごとに互いに異ならせた前記初期時刻と、当該初期時刻からの前記演算時間幅とを、演算条件として供給する演算条件供給部と、
供給された前記初期時刻に基づいて、複数の前記演算装置がそれぞれ算出する前記初期値に応じた演算結果を取得する演算結果取得部と、
前記演算結果取得部が取得する複数の前記初期値に応じた演算結果を、前記初期時刻の順に統合する統合部と、
を備える演算管理装置と、
を含む演算システム。
【請求項2】
前記逆ラプラス変換法とは、ＦＩＬＴ法である請求項１に記載の演算システム。
【請求項3】
時間変化する電磁界強度の波形についての複素周波数領域の関数に対する演算の時間軸における起点である初期時刻と、当該初期時刻からの演算時間幅とを、演算条件として取得する演算条件取得部と、
前記演算条件取得部によって取得される前記初期時刻における前記波形についての演算結果である初期値を、前記関数に対する逆ラプラス変換法による離散的な時間応答解析によって算出するとともに、算出された前記初期値に基づく時間軸に沿った逐次演算によって、前記演算条件取得部によって取得される前記演算時間幅ぶんの前記波形についての演算結果を算出する演算部と、
前記演算時間幅毎に前記演算部によって算出された前記演算結果を演算管理装置に対して出力する出力部と、
を備える演算装置。
【請求項4】
演算装置が備えるコンピュータに、時間変化する電磁界強度の波形についての複素周波数領域の関数に対する演算の時間軸における起点である初期時刻と、当該初期時刻からの演算時間幅とを、演算条件として取得する演算条件取得ステップと、
前記演算条件取得ステップにおいて取得される前記初期時刻における前記波形についての演算結果である初期値を、前記関数に対する逆ラプラス変換法による離散的な時間応答解析によって算出するとともに、算出された前記初期値に基づく時間軸に沿った逐次演算によって、前記演算条件取得ステップにおいて取得される前記演算時間幅ぶんの前記波形についての演算結果を算出する演算ステップと、
前記演算ステップにおいて前記演算時間幅毎に算出された前記演算結果を演算管理装置に対して出力する出力ステップと、
を実行させるためのプログラム。
【発明の詳細な説明】
【技術分野】
【0001】
本発明は、演算システム、演算装置及びプログラムに関する。
【背景技術】
【0002】
近年、ＦＤＴＤ法（Ｆｉｎｉｔｅ－Ｄｉｆｆｅｒｅｎｃｅ Ｔｉｍｅ－Ｄｏｍａｉｎ ｍｅｔｈｏｄ； ＦＤＴＤ ｍｅｔｈｏｄ）などの、Ｍａｘｗｅｌｌ方程式（以下、マクスウェルの方程式とも称する。）を時間・空間で差分化して解く有限差分時間領域法を用いた電磁界解析手法がある（例えば、特許文献１を参照）。この従来技術においては、複数台の演算装置によって並列演算することにより、大規模な電磁界解析を行うことが開示されている。
【先行技術文献】
【特許文献】
【0003】
【文献】特開２００９－０５３０７５号公報
【発明の概要】
【発明が解決しようとする課題】
【0004】
しかしながら、上述の従来技術によると、並列演算の結果が他の演算に影響することがあり、このような場合には、並列化された他の演算の結果を待つ必要が生じるため、複数台の演算装置による同時並列演算ができないという問題があった。
【0005】
本発明は、上記問題を解決すべくなされたもので、その目的は、電磁界解析の演算を同時並列化することができる演算装置、演算管理装置及びプログラムを提供することにある。
【課題を解決するための手段】
【0006】
本発明の一実施形態は、時間変化する電磁界強度の波形についての複素周波数領域の関数に対する演算の時間軸における起点である初期時刻と、当該初期時刻からの演算時間幅とを、演算条件として取得する演算条件取得部と、前記演算条件取得部によって取得される前記初期時刻における前記波形についての演算結果である初期値を、前記関数に対する逆ラプラス変換法による離散的な時間応答解析によって算出するとともに、算出された前記初期値に基づく時間軸に沿った逐次演算によって、前記演算条件取得部によって取得される前記演算時間幅ぶんの前記波形についての演算結果を算出する演算部と、前記演算時間幅毎に前記演算部によって算出された前記演算結果を演算管理装置に対して出力する出力部と、を備える演算装置と、それぞれの前記演算装置に問い合わせることにより取得した前記演算装置ごとの前記演算部の演算性能に基づいて、前記初期時刻からの演算時間幅を前記演算装置ごとに算出する演算時間幅算出部と、複数の前記演算装置に対して、前記演算装置ごとに互いに異ならせた前記初期時刻と、当該初期時刻からの前記演算時間幅とを、演算条件として供給する演算条件供給部と、供給された前記初期時刻に基づいて、複数の前記演算装置がそれぞれ算出する前記初期値に応じた演算結果を取得する演算結果取得部と、前記演算結果取得部が取得する複数の前記初期値に応じた演算結果を、前記初期時刻の順に統合する統合部と、を備える演算管理装置と、を含む演算システムである。
【0007】
また、本発明の一実施形態は、上述の演算システムにおいて、前記逆ラプラス変換法とは、ＦＩＬＴ法である。
【0008】
また、本発明の一実施形態は、時間変化する電磁界強度の波形についての複素周波数領域の関数に対する演算の時間軸における起点である初期時刻と、当該初期時刻からの演算時間幅とを、演算条件として取得する演算条件取得部と、前記演算条件取得部によって取得される前記初期時刻における前記波形についての演算結果である初期値を、前記関数に対する逆ラプラス変換法による離散的な時間応答解析によって算出するとともに、算出された前記初期値に基づく時間軸に沿った逐次演算によって、前記演算条件取得部によって取得される前記演算時間幅ぶんの前記波形についての演算結果を算出する演算部と、前記演算時間幅毎に前記演算部によって算出された前記演算結果を演算管理装置に対して出力する出力部と、を備える演算装置である。
【0012】
また、本発明の一実施形態は、演算装置が備えるコンピュータに、時間変化する電磁界強度の波形についての複素周波数領域の関数に対する演算の時間軸における起点である初期時刻と、当該初期時刻からの演算時間幅とを、演算条件として取得する演算条件取得ステップと、前記演算条件取得ステップにおいて取得される前記初期時刻における前記波形についての演算結果である初期値を、前記関数に対する逆ラプラス変換法による離散的な時間応答解析によって算出するとともに、算出された前記初期値に基づく時間軸に沿った逐次演算によって、前記演算条件取得ステップにおいて取得される前記演算時間幅ぶんの前記波形についての演算結果を算出する演算ステップと、前記演算ステップにおいて前記演算時間幅毎に算出された前記演算結果を演算管理装置に対して出力する出力ステップと、を実行させるためのプログラムである。
【発明の効果】
【0013】
この発明によれば、電磁界解析の演算を同時並列化することができる演算システム、演算装置及びプログラムを提供することができる。
【図面の簡単な説明】
【0014】
【図1】本実施形態の演算装置の構成の一例を示す図である。
【図2】本実施形態の演算装置の演算結果の一例を示す図である。
【図3】本実施形態の演算部による演算に用いられる微小直方体の一例を示す図である。
【図4】電界及び磁界の時間配置の割り当ての一例を示す図である。
【図5】本実施形態の解析対象物の一例を示す図である。
【図6】本実施形態の演算システムの構成の一例を示す図である。
【図7】本実施形態の演算システムの動作の一例を示す図である。
【図8】本実施形態の演算装置による並列演算の演算結果一例を示す図である。
【図9】本実施形態の複数の演算装置によって算出される応答波形の具体例を示す図である。
【図10】従来手法により１台の演算装置によって算出される応答波形の具体例を示す図である。
【発明を実施するための形態】
【0015】
【0016】
［演算装置の構成］
演算装置１０は、例えばパーソナルコンピュータなどの情報処理装置であって、供給される情報に基づいて演算を行い、演算結果を出力する。
演算装置１０による演算の一例について説明する。この一例における演算装置１０は、解析領域内のある位置における応答波形ｆｗを算出する。例えば、演算装置１０は、電磁波を放射するアンテナを波源とし、波源から放射される電磁波が入射する物体（例えば、球体や円筒体）を散乱体として、この波源及び散乱体を含む解析領域が与えられた場合に、解析領域内のある位置における電界強度の時間変化である応答波形ｆｗを算出する。
なお、以下の説明において、応答波形ｆｗを電磁界強度の時間変化を示すグラフとして表現するが、応答波形ｆｗは、電磁界強度の時間変化を示す情報であればその表現形式は問わない。例えば、応答波形ｆｗは、電磁界の空間分布や電磁界を、２次元平面内や３次元空間内の分布図などの形式によって示す情報であってもよい。
また、以下の説明において、解析領域に散乱体が含まれる場合を一例にして説明するが、これに限られない。例えば、解析領域には、電波を吸収する吸収体や、電波を放射する放射体などが含まれていてもよい。
演算装置１０は、例えばパーソナルコンピュータなどの情報処理装置であって、供給される情報に基づいて演算を行い、演算結果を出力する。
演算装置１０による演算の一例について説明する。この一例における演算装置１０は、解析領域内のある位置における応答波形ｆｗを算出する。例えば、演算装置１０は、電磁波を放射するアンテナを波源とし、波源から放射される電磁波が入射する物体（例えば、球体や円筒体）を散乱体として、この波源及び散乱体を含む解析領域が与えられた場合に、解析領域内のある位置における電界強度の時間変化である応答波形ｆｗを算出する。
なお、以下の説明において、応答波形ｆｗを電磁界強度の時間変化を示すグラフとして表現するが、応答波形ｆｗは、電磁界強度の時間変化を示す情報であればその表現形式は問わない。例えば、応答波形ｆｗは、電磁界の空間分布や電磁界を、２次元平面内や３次元空間内の分布図などの形式によって示す情報であってもよい。
また、以下の説明において、解析領域に散乱体が含まれる場合を一例にして説明するが、これに限られない。例えば、解析領域には、電波を吸収する吸収体や、電波を放射する放射体などが含まれていてもよい。
【0017】
より具体的には、演算装置１０には、演算対象情報ＭＤＬと、初期時刻ｔ０と、演算時間幅ｔｗとが、上位装置（不図示）から演算条件ＯＣとして供給される。演算装置１０は、供給される演算条件ＯＣに基づいて応答波形ｆｗを求める演算を行い、求めた応答波形ｆｗを演算結果ＲＳとして上位装置に出力する。この一例の場合、演算対象情報ＭＤＬには、波源から放射される電磁波の電磁界分布、散乱体の形状、相対位置、誘電率、透磁率など応答波形の解析に必要な諸情報が含まれる。演算装置１０による演算結果ＲＳの一例について、図２を参照して説明する。
【0018】
図２は、本実施形態の演算装置１０の演算結果ＲＳの一例を示す図である。演算装置１０は、応答波形ｆｗを演算結果ＲＳとして出力する。この一例において、応答波形ｆｗとは、解析領域内のある位置における電界強度の時間変化である。この場合、応答波形ｆｗは、横軸が時間（単位：ｆｓ（フェムト秒））、縦軸が電界強度（単位：Ｖ／ｍ）によって示される。
演算装置１０は、応答波形ｆｗのうち、ある演算時間幅ｔｗ（例えば、初期時刻ｔ０から時刻ｔ１まで）についての波形を算出する。ここで、初期時刻ｔ０とは、演算時間幅の開始時刻である。初期時刻ｔ０における応答波形ｆｗの値（例えば、電界強度）を初期値ｆ０ともいう。
つまり、演算装置１０は、初期時刻ｔ０から演算時間幅ｔｗぶんの応答波形ｆｗを算出し、算出した応答波形ｆｗを演算結果ＲＳとして出力する。
演算装置１０は、応答波形ｆｗのうち、ある演算時間幅ｔｗ（例えば、初期時刻ｔ０から時刻ｔ１まで）についての波形を算出する。ここで、初期時刻ｔ０とは、演算時間幅の開始時刻である。初期時刻ｔ０における応答波形ｆｗの値（例えば、電界強度）を初期値ｆ０ともいう。
つまり、演算装置１０は、初期時刻ｔ０から演算時間幅ｔｗぶんの応答波形ｆｗを算出し、算出した応答波形ｆｗを演算結果ＲＳとして出力する。
【0019】
図１に戻り演算装置１０の構成について説明を続ける。演算装置１０は、演算条件取得部１１０と、演算部１２０と、出力部１３０とを備える。
演算条件取得部１１０は、上位装置（不図示）から演算条件ＯＣとして供給される演算対象情報ＭＤＬと、初期時刻ｔ０と、演算時間幅ｔｗとを取得する。
演算部１２０は、演算条件取得部１１０が取得する演算条件ＯＣに基づいて、初期値ｆ０と応答波形ｆｗとを算出する。
出力部１３０は、算出した初期値ｆ０と応答波形ｆｗとを、上位装置（不図示）に対して出力する。
演算条件取得部１１０は、上位装置（不図示）から演算条件ＯＣとして供給される演算対象情報ＭＤＬと、初期時刻ｔ０と、演算時間幅ｔｗとを取得する。
演算部１２０は、演算条件取得部１１０が取得する演算条件ＯＣに基づいて、初期値ｆ０と応答波形ｆｗとを算出する。
出力部１３０は、算出した初期値ｆ０と応答波形ｆｗとを、上位装置（不図示）に対して出力する。
【0020】
［演算の具体例］
演算部１２０が行う演算の具体例について説明する。本実施形態の演算部１２０は、初期時刻ｔ０における初期値ｆ０の算出と、初期時刻ｔ０以降における応答波形ｆｗの算出とを、それぞれ異なる算出方法によって行う。まず、初期時刻ｔ０以降における応答波形ｆｗの算出について説明する。
演算部１２０が行う演算の具体例について説明する。本実施形態の演算部１２０は、初期時刻ｔ０における初期値ｆ０の算出と、初期時刻ｔ０以降における応答波形ｆｗの算出とを、それぞれ異なる算出方法によって行う。まず、初期時刻ｔ０以降における応答波形ｆｗの算出について説明する。
【0021】
［演算例（１）ＦＤＴＤ法］
演算部１２０は、ＦＤＴＤ法によって応答波形ｆｗを算出する。
演算部１２０は、ＦＤＴＤ法によって応答波形ｆｗを算出する。
【0022】
図３は、本実施形態の演算部１２０による演算に用いられる微小直方体の一例を示す図である。演算部１２０は、与えられた解析領域全体を、同図に示す微小直方体（例えば、Ｙｅｅ格子）によって分割する。ここで、同図のＥ（斜体）は電界を、Ｈ（斜体）は磁界を示す。電界Ｅ及び磁界Ｈの添え字ｘｙｚ（いずれも斜体）は、それぞれｘ軸方向成分、ｙ軸方向成分、ｚ軸方向成分を示す。演算部１２０は、分割された微小直方体のそれぞれについてＭａｘｗｅｌｌ方程式（式（１）及び式（２））を適用して、時間及び３次元空間のそれぞれについて定式化した結果を用いて、演算を行う。具体的には次の通りである。
【0023】
【数1】
【0024】
【数2】
【0025】
まず、時間についての差分を求める。式（１）及び式（２）に示すファラデーの法則及びアンペアの法則の時間微分項についてまとめると式（３）及び式（４）のようになる。
【0026】
【数3】
【0027】
【数4】
【0028】
図４は、電界Ｅ及び磁界Ｈの時間配置の割り当ての一例を示す図である。
次に、図４に示すような電界Ｅ及び磁界Ｈの時間配置を、時間軸ｔに対して割り当て、中心差分を行うと、式（３）及び式（４）は、式（５）及び式（６）となる。
次に、図４に示すような電界Ｅ及び磁界Ｈの時間配置を、時間軸ｔに対して割り当て、中心差分を行うと、式（３）及び式（４）は、式（５）及び式（６）となる。
【0029】
【数5】
【0030】
【数6】
【0031】
ここで式（６）のＥｎ－１／２は定義されていないので式（７）のように近似する。
【0032】
【数7】
【0033】
式（６）のＥｎ－１／２を式（７）によって近似し、式（６）の右辺第一項に代入することにより、式（８）が得られる。
【0034】
【数8】
【0035】
式（５）及び式（８）についてまとめると、式（９）～式（１１）になる。
【0036】
【数9】
【0037】
【数10】
【0038】
【数11】
【0039】
以上により、時間についての定式化が完了する。
次に３次元空間についての定式化を行う．式（９）と式（１０）を（電界Ｅｘ、電界Ｅｙ、電界Ｅｚ、磁界Ｈｘ、磁界Ｈｙ、磁界Ｈｚ）の６成分に分け図４に示した電磁界の空間配置を適用して定式化すると以下のようになる。なお、以下の式（１２）～式（１９）において電界Ｅ及び磁界Ｈの記号には斜体表記を用いる。
次に３次元空間についての定式化を行う．式（９）と式（１０）を（電界Ｅｘ、電界Ｅｙ、電界Ｅｚ、磁界Ｈｘ、磁界Ｈｙ、磁界Ｈｚ）の６成分に分け図４に示した電磁界の空間配置を適用して定式化すると以下のようになる。なお、以下の式（１２）～式（１９）において電界Ｅ及び磁界Ｈの記号には斜体表記を用いる。
【0040】
【数12】
【0041】
【数13】
【0042】
【数14】
【0043】
【数15】
【0044】
【数16】
【0045】
【数17】
【0046】
【数18】
【0047】
【数19】
【0048】
【数20】
【0049】
【数21】
【0050】
演算部１２０は、以上のようにして時間及び３次元空間について定式化された演算式を、与えられた初期値ｆ０に対して適用することにより、初期時刻ｔ０以降における応答波形ｆｗを算出する。
【0051】
［演算例（２）数値逆ラプラス変換法（ＦＩＬＴ法）］
演算部１２０は、数値逆ラプラス変換法（ＦＩＬＴ法；Ｆａｓｔ Ｉｎｖｅｒｓｅ Ｌａｐｌａｃｅ Ｔｒａｎｓｆｏｒｍ ｍｅｔｈｏｄ）によって応答波形ｆｗを算出する。なお、数値逆ラプラス変換法は、ＮＩＬＴ法（Ｎｕｍｅｒｉｃａｌ Ｉｎｖｅｒｓｅ Ｌａｐｌａｃｅ Ｔｒａｎｓｆｏｒｍ ｍｅｔｈｏｄ)とも称される。
複素周波数領域の関数Ｆ（ｓ）は、数値逆ラプラス変換法（ＦＩＬＴ法）により時間領域に変換できる。この数値逆ラプラス変換法（ＦＩＬＴ法）は、逆ラプラス変換法による離散的な時間応答解析の一例である。一般的に、逆ラプラス変換はＢｒｏｍｗｉｃｈ積分より以下の式で定義される。
演算部１２０は、数値逆ラプラス変換法（ＦＩＬＴ法；Ｆａｓｔ Ｉｎｖｅｒｓｅ Ｌａｐｌａｃｅ Ｔｒａｎｓｆｏｒｍ ｍｅｔｈｏｄ）によって応答波形ｆｗを算出する。なお、数値逆ラプラス変換法は、ＮＩＬＴ法（Ｎｕｍｅｒｉｃａｌ Ｉｎｖｅｒｓｅ Ｌａｐｌａｃｅ Ｔｒａｎｓｆｏｒｍ ｍｅｔｈｏｄ)とも称される。
複素周波数領域の関数Ｆ（ｓ）は、数値逆ラプラス変換法（ＦＩＬＴ法）により時間領域に変換できる。この数値逆ラプラス変換法（ＦＩＬＴ法）は、逆ラプラス変換法による離散的な時間応答解析の一例である。一般的に、逆ラプラス変換はＢｒｏｍｗｉｃｈ積分より以下の式で定義される。
【0052】
【数22】
【0053】
ＦＩＬＴ法では式（２２）の指数関数について、以下の式（２３）による近似を行う。
【0054】
【数23】
【0055】
【数24】
【0056】
ただし、
【0057】
【数25】
【0058】
である。
式（２４）の無限級数を数値計算のため有限の項数で打ち切り、式（２６）を得る。
式（２４）の無限級数を数値計算のため有限の項数で打ち切り、式（２６）を得る。
【0059】
【数26】
【0060】
ここで、式（２６）におけるＭ（斜体）は、ＦＩＬＴ法における打ち切り項数である。式（２６）は、交代級数であるためＥｕｌｅｒ変換を適用することにより、式（２７）に変換される。
【0061】
【数27】
【0062】
ここで、
【0063】
【数28】
【0064】
であり、ｐ（斜体）は、Ｅｕｌｅｒ変換の項数を示す。式（２７）は、ある観測時刻において他の観測時刻とは互いに独立に計算が可能であり、各観測時刻での電磁界を容易に並列計算できる。
また、ＦＩＬＴ法における誤差は式（２３）の近似より以下の式（２９）で表せる。
また、ＦＩＬＴ法における誤差は式（２３）の近似より以下の式（２９）で表せる。
【0065】
【数29】
【0066】
式（２９）より、近似のパラメータに対し、計算誤差は指数関数的に減少する。また、主に第一項が計算誤差となり、その誤差は、１０－αとして見積れる。
【0067】
［演算例（３）ＦＤＣＦＤ－ＦＩＬＴ法］
一般的なＦＤＦＤ法（Ｆｉｎｉｔｅ－Ｄｉｆｆｅｒｅｎｃｅ Ｆｒｅｑｕｅｎｃｙ－Ｄｏｍａｉｎ ｍｅｔｈｏｄ）では、上述したＦＤＴＤ法と同様の空間メッシュ分割を用いつつ、周波数領域において行列方程式を解く。本実施形態の演算部１２０は、ＦＤＦＤ法を複素周波数領域に拡張したＦＤＣＦＤ法（Ｆｉｎｉｔｅ－Ｄｉｆｆｅｒｅｎｃｅ Ｃｏｍｐｌｅｘ Ｆｒｅｑｕｅｎｃｙ－Ｄｏｍａｉｎ ｍｅｔｈｏｄ）と上述のＦＩＬＴ法とを用いる。この複素周波数領域でのＭａｘｗｅｌｌ方程式は、式（３０）及び式（３１）によって示される。
一般的なＦＤＦＤ法（Ｆｉｎｉｔｅ－Ｄｉｆｆｅｒｅｎｃｅ Ｆｒｅｑｕｅｎｃｙ－Ｄｏｍａｉｎ ｍｅｔｈｏｄ）では、上述したＦＤＴＤ法と同様の空間メッシュ分割を用いつつ、周波数領域において行列方程式を解く。本実施形態の演算部１２０は、ＦＤＦＤ法を複素周波数領域に拡張したＦＤＣＦＤ法（Ｆｉｎｉｔｅ－Ｄｉｆｆｅｒｅｎｃｅ Ｃｏｍｐｌｅｘ Ｆｒｅｑｕｅｎｃｙ－Ｄｏｍａｉｎ ｍｅｔｈｏｄ）と上述のＦＩＬＴ法とを用いる。この複素周波数領域でのＭａｘｗｅｌｌ方程式は、式（３０）及び式（３１）によって示される。
【0068】
【数30】
【0069】
【数31】
【0070】
このＦＤＣＦＤ法のアルゴリズムの定式化について、２次元の場合について説明する。
計算空間を離散化するために、式（３０）及び式（３１）は次の式（３２）～式（３５）に示す形式に書き換えることができる。
計算空間を離散化するために、式（３０）及び式（３１）は次の式（３２）～式（３５）に示す形式に書き換えることができる。
【0071】
【数32】
【0072】
【数33】
【0073】
【数34】
【0074】
【数35】
【0075】
計算空間全体に対して同様の手順を適用するために、式（３６）に示す行列形式の連立方程式を用いる。
【0076】
【数36】
【0077】
ここで、式（３７）に示すように、金属の周波数分散を表すために、複素周波数s領域の複素誘電率をＬｏｒｅｎｔｚ－Ｄｒｕｄｅによって示す。
【0078】
【数37】
【0079】
ここで、ωpはプラズマ周波数、ｋは、共振周波数ωj、媒質固有の定数ｆｊ、衝突周波数Γｊである振動子の数である（式中の記号はいずれも斜体）。上述の連立方程式を解くことにより、複素周波数領域の電磁界を求めることができる。複素周波数領域における式（３６）の解ｘ（斜体）は、上述したＦＩＬＴ法を適用することによって時間領域に変換することができる。ここで、本実施形態のアルゴリズムにおいて、Ｂｒｏｍｗｉｃｈ積分の指数関数は次の式（３８）によって近似することができる。
【0080】
【数38】
【0081】
ここで、αは、近似パラメータである。Ｂｒｏｍｗｉｃｈ積分に式（３８）を代入すると、式（３９）が得られる。
【0082】
【数39】
【0083】
ここで、
【0084】
【数40】
【0085】
【数41】
【0086】
である。式（３９）の無限級数を、数値計算のため有限の項数で打ち切り、Ｅｕｌｅｒ変換を適用することにより、逆ラプラス変換を示す式（４２）を得る。
【0087】
【数42】
【0088】
【数43】
【0089】
【数44】
【0090】
［演算例（４）複素周波数領域の積分方程式法］
本演算例において、散乱体はＰＥＣ（Ｐｅｒｆｅｃｔ Ｅｌｅｃｔｒｉｃ Ｃｏｎｄｕｃｔｏｒ）球であり、導体の厚さはゼロであり、開口は球面座標の範囲内に存在する（すなわち、球表面の一部に開口がある）と仮定する。このＰＥＣ球に対する入射波を式（４５）に示す。
本演算例において、散乱体はＰＥＣ（Ｐｅｒｆｅｃｔ Ｅｌｅｃｔｒｉｃ Ｃｏｎｄｕｃｔｏｒ）球であり、導体の厚さはゼロであり、開口は球面座標の範囲内に存在する（すなわち、球表面の一部に開口がある）と仮定する。このＰＥＣ球に対する入射波を式（４５）に示す。
【0091】
【数45】
【0092】
ここで、
【0093】
【数46】
【0094】
である。また、正規化された複素周波数Ｓ＝ｓ（ε０μ０）１／２ａ、正規化された距離Ｒ＝ｒ／ａ、ａは散乱体を取り囲む球の半径、θｉｎ及びφｉｎは入射角、Ｅ０＾（イーゼロ・ハット）は入射波のスペクトルである。
この一例の場合、複素周波数領域における電界積分方程式は、式（４７）で与えられる。
この一例の場合、複素周波数領域における電界積分方程式は、式（４７）で与えられる。
【0095】
【数47】
【0096】
ここで、Ｊ（Ｒ）は未知の表面電流分布、ｔ＾（ティー・ハット）は散乱体の表面上の単位接線ベクトルであり、
【0097】
【数48】
【0098】
【数49】
【0099】
である。散乱体の表面は三角形領域で分割され、式（５１）に示すＲＷＧ（Ｒａｏ‐Ｗｉｌｔｏｎ‐Ｇｌｉｓｓｏｎ）基底関数などを使用して拡張される。すなわち、電流密度Ｊは、ＲＷＧ基底関数などを使用して式（５０）のように展開される。
【0100】
【数50】
【0101】
ここで、Ｉｎは未知の展開係数、Ｎは未知数の数であり、
【0102】
【数51】
【0103】
である。ここでＡｎ
＋、Ａｎ
－は三角形Ｔｎ
＋及び三角形Ｔｎ
－の面積、ρｎ
＋及びρｎ
－は三角形Ｔｎ
＋及び三角形Ｔｎ
－の頂点の位置ベクトル、ｌｎは辺の長さである。
式（４７）は、Ｉｎを基に式（５０）に置き換えることにより離散化され（すなわち、式（５０）のＪ（Ｒ）を式（５１）の基底関数ｆｎ（Ｒ）と式（５０）の展開係数Ｉｎを基に置き換えることにより離散化され）、一組のテスト関数が乗算される。既知のモーメント法の手順と同様にして内積を求めることにより、表面電流を得る。表面電流分布に基づき、観測角（θ、φ）における散乱電界Ｅθ （Ｓ）は、式（５２）によって示される。
式（４７）は、Ｉｎを基に式（５０）に置き換えることにより離散化され（すなわち、式（５０）のＪ（Ｒ）を式（５１）の基底関数ｆｎ（Ｒ）と式（５０）の展開係数Ｉｎを基に置き換えることにより離散化され）、一組のテスト関数が乗算される。既知のモーメント法の手順と同様にして内積を求めることにより、表面電流を得る。表面電流分布に基づき、観測角（θ、φ）における散乱電界Ｅθ （Ｓ）は、式（５２）によって示される。
【0104】
【数52】
【0105】
ここで、
【0106】
【数53】
【0107】
【数54】
【0108】
この複素周波数領域の関数は、上述したＦＩＬＴ法を適用することにより、時間領域に数値的に変換される。
【0109】
［演算例（５）静的近似を用いた複素周波数領域の積分方程式］
図５は、本実施形態の解析対象物の一例を示す図である。本演算例において、解析対象物は同図に示す一様な金属で構成される任意形状物体とする。具体的には、同図に示すように、解析対象物は、任意形状の金属を分割した金属物体ｍｔで構成され、そのサイズは入射波長に対し０．１倍以下とする。
金属の誘電率は、Ｌｏｒｅｎｔｚ－Ｄｒｕｄｅモデルとし式（５５）により定義する。
図５は、本実施形態の解析対象物の一例を示す図である。本演算例において、解析対象物は同図に示す一様な金属で構成される任意形状物体とする。具体的には、同図に示すように、解析対象物は、任意形状の金属を分割した金属物体ｍｔで構成され、そのサイズは入射波長に対し０．１倍以下とする。
金属の誘電率は、Ｌｏｒｅｎｔｚ－Ｄｒｕｄｅモデルとし式（５５）により定義する。
【0110】
【数55】
【0111】
ここで、sは複素周波数、ωｐはプラズマ周波数、ωｊは共振周波数、ν０とν１とは衝突周波数、Ａ０とＡｉとは媒質固有の定数、Ｋは振動子の個数である（式中の記号はいずれも斜体）。解析対象物の誘電率が負であり、解析対象物のサイズが入射波長に比べ十分小さい場合には、解析対象物近傍の電磁界が増強されるプラズモン共鳴が観測できる。
定式化を行うにあたり、散乱電磁界および入射界を式（５６）で表わす。
定式化を行うにあたり、散乱電磁界および入射界を式（５６）で表わす。
【0112】
【数56】
【0113】
解析対象物が入射波長に比べ十分小さいとし、静的近似を用いることで式（５６）を式（５７）、式（５８）に展開する。
【0114】
【数57】
【0115】
【数58】
【0116】
ここで、β＝ｓ（ε０μ０）１／２ｄ、ｄは解析対象物のサイズを表す。式（５７）及び式（５８）において０次項は式（５９）の境界条件を満足する。
【0117】
【数59】
【0118】
ただし、Ｅ０
±＝ε－１／２ｅ０
±、ｆｉｎ＾（ｓ）（エフイン・ハット・エス）は入射波のスペクトル、ｎは単位法線ベクトルを示す。また、この静的近似により、解析対象物の境界面における電界の法線方向成分は未知表面電荷密度σ＾（シグマ・ハット）を用いて以下の式（６０）で表わせる。
【0119】
【数60】
【0120】
ただし、Ωは解析対象物表面を示す。
次に、境界面における電束密度の法線方向成分連続の条件より、式（６１）に示す境界型積分方程式が得られる。
次に、境界面における電束密度の法線方向成分連続の条件より、式（６１）に示す境界型積分方程式が得られる。
【0121】
【数61】
【0122】
【数62】
【0123】
である。
次に、未知電荷密度を基底関数ｊ＾ｉ（ジェイ・ハット・アイ）により近似展開する。
次に、未知電荷密度を基底関数ｊ＾ｉ（ジェイ・ハット・アイ）により近似展開する。
【0124】
【数63】
【0125】
ａｉは未知展開係数、Ｎは未知数の数を示す。式（６１）を式（６２）と試行関数ｔ＾ｉ（ティー・ハット・アイ）により離散化し、式（６４）の連立一次方程式が得られる。
【0126】
【数64】
【0127】
【数65】
【0128】
【数66】
【0129】
式（６４）の連立一次方程式の求解より、未知電荷密度を得る。
連立一次方程式の求解に反復法を用いた場合、ＢＩＥＭ（Ｂｏｕｎｄａｒｙ Ｉｎｔｅｇｒａｌ Ｅｑｕａｔｉｏｎ Ｍｅｔｈｏｄ；境界型積分方程式法）では計算時間と計算機メモリが未知数Ｎ（斜体）の２乗のオーダＯ（Ｎ＾２）に比例する（ここで、Ｎ＾２はＮの２乗を示す）。ＭＭ（Ｆａｓｔ Ｍｕｌｔｉｐｏｌｅ Ｍｅｔｈｏｄ；高速多重極法）による高速化を行うためには、解析対象物をボックス分割する（図５を参照）。ここで、自由空間中のグリーン関数を以下の式（６７）で表現する。
連立一次方程式の求解に反復法を用いた場合、ＢＩＥＭ（Ｂｏｕｎｄａｒｙ Ｉｎｔｅｇｒａｌ Ｅｑｕａｔｉｏｎ Ｍｅｔｈｏｄ；境界型積分方程式法）では計算時間と計算機メモリが未知数Ｎ（斜体）の２乗のオーダＯ（Ｎ＾２）に比例する（ここで、Ｎ＾２はＮの２乗を示す）。ＭＭ（Ｆａｓｔ Ｍｕｌｔｉｐｏｌｅ Ｍｅｔｈｏｄ；高速多重極法）による高速化を行うためには、解析対象物をボックス分割する（図５を参照）。ここで、自由空間中のグリーン関数を以下の式（６７）で表現する。
【0130】
【数67】
【0131】
ただし、ｒｍとｒｍ’とはグループｍとｍ’との中心までの距離ベクトルであり、
【0132】
【数68】
【0133】
【数69】
【0134】
である。ＬはＦＭＭにおける打ち切り項数を示す。要素間の相互作用は、互いに近接したボックス内に属する場合は式（６４）の行列要素を直接計算する。それ以外の場合は、式（６７）を用いて計算をグループ化して行うことが可能になる。これにより、ボックス分割を２回行うレベル２における計算量はＯ（Ｎ１．５）、さらにボックス分割を行う多重レベル拡張時ではＯ（Ｎ）まで計算量を低減できる。
また、固有モード解析のため、式（６１）は以下の固有方程式に変換される。
また、固有モード解析のため、式（６１）は以下の固有方程式に変換される。
【0135】
【数70】
【0136】
ここで、
【0137】
【数71】
【0138】
であり、ｋはモードナンバー、λｋは固有値であり、式（７０）を解くことで得る。また、共振時における誘電率εｋは式（７１）の関係より得られる。表面電荷の時間変化を表わすため、電荷密度を以下の式で展開する。
【0139】
【数72】
【0140】
ここで、ａｋ（ｔ）は展開係数を示す。ラプラス変換により式（７２）を複素周波数領域にて展開することで次式を得る。
【0141】
【数73】
【0142】
ここで、
【0143】
【数74】
【0144】
【数75】
【0145】
これより得られたσｋ（ｒ、ｓ）をＦＩＬＴ法により時間領域に変換する。
【0146】
［演算例（６）複素周波数領域の厳密解］
自由空間中に置かれたz軸に一様な二次元誘電体円柱の問題を解析する場合について説明する。入射波はz軸方向にのみ磁界成分を持つＨ波とし、円柱を構成する様々な媒質に対して、複素周波数領域における厳密解を求める。入射平面波Ｈｚ（ｉ）に対する散乱電磁界Ｈｚ（ｓ），Ｅφ（ｓ）、Ｅｒ（ｓ）は式（７６）～式（８０）によって表せる。
自由空間中に置かれたz軸に一様な二次元誘電体円柱の問題を解析する場合について説明する。入射波はz軸方向にのみ磁界成分を持つＨ波とし、円柱を構成する様々な媒質に対して、複素周波数領域における厳密解を求める。入射平面波Ｈｚ（ｉ）に対する散乱電磁界Ｈｚ（ｓ），Ｅφ（ｓ）、Ｅｒ（ｓ）は式（７６）～式（８０）によって表せる。
【0147】
【数76】
【0148】
【数77】
【0149】
【数78】
【0150】
【数79】
【0151】
【数80】
【0152】
ここで、Ｈ０（ｓ）：入射波形、Ａｎ：散乱係数、Ｉｎ（・）：第一種変形Ｂｅｓｓｅｌ関数、Ｋｎ（・）：第二種変形Ｂｅｓｓｅｌ関数、ｒ：観測距離、θ：観測角、φ：入射角、ε０：真空中の誘電率、μ０：真空中の透磁率、Ｚ０：真空中の波動インピーダンスである（いずれも斜体）。散乱係数Ａｎは、誘電体円柱を構成する媒質に対する境界条件を考慮することで、以下のように求まる。
【0153】
１）完全導体の場合
【0154】
【数81】
【0155】
ここで、a：円柱の半径である。
【0156】
２） 誘電体の場合
【0157】
【数82】
【0158】
【数83】
【0159】
ここで、ε１：円柱の誘電率、εｒ：円柱の比誘電率、μ１：円柱の透磁率、μｒ：円柱の比透磁率である。
【0160】
３） 完全導体に誘電体コーティングした場合
【0161】
【数84】
【0162】
【数85】
【0163】
【数86】
【0164】
ここで、ｂ：内部円柱の半径、ａ：外部円柱の半径、ε１：外部円柱の誘電率、εｒ１：外部円柱の比誘電率、μ１：外部円柱の透磁率、μｒ１：外部円柱の比透磁率である。
【0165】
４） ２層誘電体の場合
【0166】
【数87】
【0167】
【数88】
【0168】
【数89】
【0169】
【数90】
【0170】
ここで、ε２：内部円柱の誘電率、εｒ２：内部円柱の比誘電率、μ２：内部円柱の透磁率、μｒ２：内部円柱の比透磁率である。
【0171】
［演算システムの全体構成］
次に、図６を参照して本実施形態の演算システム１の全体構成について説明する。
図６は、本実施形態の演算システム１の構成の一例を示す図である。演算システムは、１、複数台の演算装置１０と、演算管理装置２０と、記憶装置３０とを備える。本実施形態の一例では、演算装置１０－１から演算装置１０－ｎ（ｎは自然数。）までのｎ台の演算装置１０を利用して、応答波形ｆｗを算出する場合について説明する。
次に、図６を参照して本実施形態の演算システム１の全体構成について説明する。
図６は、本実施形態の演算システム１の構成の一例を示す図である。演算システムは、１、複数台の演算装置１０と、演算管理装置２０と、記憶装置３０とを備える。本実施形態の一例では、演算装置１０－１から演算装置１０－ｎ（ｎは自然数。）までのｎ台の演算装置１０を利用して、応答波形ｆｗを算出する場合について説明する。
【0172】
演算管理装置２０は、演算装置１０に対して演算条件ＯＣを供給し、この演算条件ＯＣに基づいて演算装置１０が求めた演算結果ＲＳを取得する。演算管理装置２０は、複数台の演算装置１０に対して、互いに異なる初期時刻ｔ０を供給して演算させることにより、１つの応答波形ｆｗを複数台の演算装置１０に並列演算させる。
【0173】
この演算管理装置２０は、取得部２１０と、演算時間幅算出部２２０と、演算条件供給部２３０と、演算結果取得部２４０と、統合部２５０とを備える。
【0174】
取得部２１０は、不図示の上位装置から演算対象情報ＭＤＬと並列装置数ＮＰＣとを取得する。ここで、並列装置数ＮＰＣとは、応答波形ｆｗを算出する演算装置１０の並列数である。
【0175】
演算時間幅算出部２２０は、取得部２１０が取得した並列装置数ＮＰＣと、それぞれの演算装置１０の演算性能ＣＰとに基づいて、各演算装置１０に担当させる演算時間幅ｔｗを演算装置１０ごとに算出する。ここで、演算性能ＣＰは、予め記憶装置３０に記憶されていてもよいし、演算管理装置２０が演算装置１０に問い合わせることにより取得されてもよい。演算時間幅算出部２２０は、演算装置１０ごとに算出した演算時間幅ｔｗに基づいて、各演算装置１０の演算の初期時刻ｔ０を演算装置１０ごとに算出する。
【0176】
演算条件供給部２３０は、演算条件ＯＣを各演算装置１０に供給する。この演算条件ＯＣには、演算時間幅算出部２２０が算出した初期時刻ｔ０と、演算時間幅ｔｗとが含まれる。具体的には、演算条件供給部２３０は、演算装置１０－１に対して、初期時刻ｔ０－１と、演算時間幅ｔｗ－１とを含む演算条件ＯＣ１を供給する。演算条件供給部２３０は、演算装置１０－２に対して、初期時刻ｔ０－２と、演算時間幅ｔｗ－２とを含む演算条件ＯＣ２を供給する。同様に、演算条件供給部２３０は、演算装置１０－ｎに対して、初期時刻ｔ０－ｎと、演算時間幅ｔｗ－ｎとを含む演算条件ＯＣｎを供給する。
【0177】
各演算装置１０は、演算管理装置２０から供給される演算条件ＯＣに基づいて、応答波形ｆｗを演算結果ＲＳとして算出する。演算装置１０は、算出した演算結果ＲＳを演算管理装置２０に対して供給する。具体的には、演算装置１０－１は、演算結果ＲＳ１を演算管理装置２０に対して供給する。演算装置１０－２は、演算結果ＲＳ２を演算管理装置２０に対して供給する。同様にして、演算装置１０－ｎは、演算結果ＲＳｎを演算管理装置２０に対して供給する。
【0178】
演算結果取得部２４０は、各演算装置１０から供給される演算結果ＲＳを取得する。この演算結果ＲＳには、各演算装置１０に割り当てた演算時間幅ｔｗぶんの応答波形ｆｗの算出結果が含まれている。
【0179】
統合部２５０は、演算結果取得部２４０が、演算装置１０からそれぞれ取得した複数の演算結果ＲＳを統合する。
【0180】
【0181】
（ステップＳ１０）取得部２１０は、演算対象情報ＭＤＬを取得する。
（ステップＳ２０）取得部２１０は、並列装置数ＮＰＣを取得する。ここでは、一例として４台の演算装置１０によって並列演算する場合について説明する。この場合、並列装置数ＮＰＣは「４」である。
（ステップＳ３０）取得部２１０は、記憶装置３０から演算性能情報を取得する。この演算性能情報は、演算装置１０の演算性能ＣＰを演算装置１０ごとに示す。この一例では、４台の演算性能ＣＰが同等である場合について説明する。
（ステップＳ２０）取得部２１０は、並列装置数ＮＰＣを取得する。ここでは、一例として４台の演算装置１０によって並列演算する場合について説明する。この場合、並列装置数ＮＰＣは「４」である。
（ステップＳ３０）取得部２１０は、記憶装置３０から演算性能情報を取得する。この演算性能情報は、演算装置１０の演算性能ＣＰを演算装置１０ごとに示す。この一例では、４台の演算性能ＣＰが同等である場合について説明する。
【0182】
（ステップＳ４０）演算時間幅算出部２２０は、ステップＳ２０において取得された並列装置数ＮＰＣと、ステップＳ３０において取得された演算性能情報とに基づいて、演算装置１０ごとの演算時間幅ｔｗを算出する。演算時間幅算出部２２０が算出する演算時間幅ｔｗの具体例について、図８を参照して説明する。
【0183】
図８は、本実施形態の演算装置１０による並列演算の演算結果ＲＳ一例を示す図である。本実施形態の一例において、４台の演算装置１０（演算装置１０－１～演算装置１０－４）が並列演算する場合について説明する。演算管理装置２０は、開始時刻ｔｓから終了時刻ｔｅまでを演算対象として、演算時間幅ｔｗの分割を行う。この一例の場合、４台の演算装置１０が同等の演算性能ＣＰを有している。このため、演算時間幅算出部２２０は、開始時刻ｔｓから終了時刻ｔｅまでを均等に４分割することにより、演算時間幅ｔｗを分割する。具体的には、演算時間幅算出部２２０は、開始時刻ｔｓから終了時刻ｔｅまでを演算時間幅ｔｗ－１、演算時間幅ｔｗ－２、演算時間幅ｔｗ－３、演算時間幅ｔｗ－４に４等分する。ここで、演算時間幅ｔｗ－１とは、演算装置１０－１に割り当てられる演算時間幅ｔｗである。同様に、演算時間幅ｔｗ－２～－４とは、演算装置１０－２～－４に割り当てられる演算時間幅ｔｗである。
【0184】
すなわち、演算管理装置２０は、演算時間幅算出部２２０を備えている。この演算時間幅算出部２２０は、それぞれの演算装置１０が備える演算部１２０の演算性能ＣＰに基づいて、初期時刻ｔ０からの演算時間幅ｔｗを演算装置１０ごとに算出する。
【0185】
すなわち、演算時間幅算出部２２０は、演算部１２０の演算性能ＣＰに基づいて演算時間幅ｔｗを算出する。つまり、演算時間幅ｔｗが、演算部１２０の演算性能ＣＰに基づいて定められている。
【0186】
演算時間幅算出部２２０は、分割した演算時間幅ｔｗ－１～－４を時系列に並べ、それぞれの演算時間幅ｔｗの開始時刻を、初期時刻ｔ０として算出する。具体的には、演算時間幅算出部２２０は、演算時間幅ｔｗ－１の開始時刻を初期時刻ｔ０－１とし、演算時間幅ｔｗ－２の開始時刻を初期時刻ｔ０－２とし、演算時間幅ｔｗ－３の開始時刻を初期時刻ｔ０－３とし、演算時間幅ｔｗ－４の開始時刻を初期時刻ｔ０－４とする。
【0187】
（ステップＳ５０）図７に戻り、演算条件供給部２３０は、演算対象情報ＭＤＬと、演算時間幅算出部２２０が上述のようにして算出した演算時間幅ｔｗ及び初期時刻ｔ０とを演算条件ＯＣとして各演算装置１０に供給する。より具体的には、演算時間幅算出部２２０は、演算対象情報ＭＤＬと演算時間幅ｔｗ－１及び初期時刻ｔ０－１とを演算装置１０－１に供給する。同様にして演算条件供給部２３０は、演算対象情報ＭＤＬと演算時間幅ｔｗ－２及び初期時刻ｔ０－２とを演算装置１０－２に、演算対象情報ＭＤＬと演算時間幅ｔｗ－３及び初期時刻ｔ０－３とを演算装置１０－３に、演算対象情報ＭＤＬと演算時間幅ｔｗ－４及び初期時刻ｔ０－４とを演算装置１０－４に、それぞれ供給する。
【0188】
すなわち、演算条件供給部２３０は、複数の演算装置１０に対して、演算装置１０ごとに互いに異ならせた初期時刻ｔ０を演算条件ＯＣとして供給する。
【0189】
また、演算条件供給部２３０は、演算時間幅算出部２２０が算出する演算時間幅ｔｗを、演算条件ＯＣとして演算装置１０に供給する。
【0190】
各演算装置１０は、演算管理装置２０から供給された演算条件ＯＣに基づいて、応答波形ｆｗを演算結果ＲＳとして算出する。これら演算装置１０が算出する応答波形ｆｗについて、再び図８を参照して説明する。
【0191】
各演算装置１０の演算条件取得部１１０は、初期時刻ｔ０を取得する。ここで、初期時刻ｔ０とは、時間変化する応答波形ｆｗについての演算の時間軸における起点である。
【0192】
各演算装置１０は、初期時刻ｔ０における初期値ｆ０を算出する。また、各演算装置１０は、算出した初期値ｆ０に基づいて、演算時間幅ｔｗぶんの応答波形ｆｗを算出する。より具体的には、演算装置１０－１は、初期時刻ｔ０－１における初期値ｆ０－１を算出する。また、演算装置１０－１は、初期値ｆ０－１に基づいて、演算時間幅ｔｗ－１ぶんの応答波形ｆｗ１を算出する。演算装置１０－２～－４は、演算装置１０－１と同様にして応答波形ｆｗ－２～－４を算出する。
【0193】
ここで、演算装置１０の演算部１２０は、上述したＦＤＣＦＤ－ＦＩＬＴ法によって、初期値ｆ０を算出する。すなわち、演算装置１０は、逆ラプラス変換法による離散的な時間応答解析によって初期値ｆ０を算出する。ここで、初期値ｆ０とは、演算条件取得部１１０によって取得される初期時刻ｔ０における応答波形ｆｗ（波形）についての演算結果ＲＳである。また、ここでいう逆ラプラス変換法とは、ＦＩＬＴ法である。
また、演算装置１０の出力部１３０は、演算部１２０によって算出された初期値ｆ０を出力する。
また、演算装置１０の出力部１３０は、演算部１２０によって算出された初期値ｆ０を出力する。
【0194】
また、演算装置１０の演算部１２０は、上述したＦＤＴＤ法によって、初期時刻ｔ０以降の演算時間幅ｔｗぶんの応答波形ｆｗを算出する。すなわち、演算装置１０は、算出された初期値ｆ０に基づく時間軸に沿った逐次演算によって、応答波形ｆｗについての所定の演算時間幅ｔｗぶんの演算結果ＲＳを算出する。
また、演算装置１０の出力部１３０は、演算時間幅ｔｗ毎に演算部１２０によって算出された演算結果ＲＳを出力する。
また、演算装置１０の出力部１３０は、演算時間幅ｔｗ毎に演算部１２０によって算出された演算結果ＲＳを出力する。
【0195】
（ステップＳ６０）図７に戻り、演算結果取得部２４０は、各演算装置１０が算出した初期値ｆ０及び応答波形ｆｗ、すなわち演算結果ＲＳを取得する。具体的には、演算結果取得部２４０は、演算装置１０－１から初期値ｆ０－１と応答波形ｆｗ１とを取得する。これと同様にして、演算結果取得部２４０は、演算装置１０－２から初期値ｆ０－２と応答波形ｆｗ２とを、演算装置１０－３から初期値ｆ０－３と応答波形ｆｗ３とを、演算装置１０－４から初期値ｆ０－４と応答波形ｆｗ４とを、それぞれ取得する。
【0196】
すなわち、演算結果取得部２４０は、供給された初期時刻ｔ０に基づいて、複数の演算装置１０がそれぞれ算出する初期値ｆ０を取得する。
【0197】
（ステップＳ７０）統合部２５０は、ステップＳ６０において各演算装置１０から取得された、分割された応答波形ｆｗ（この具体例では、応答波形ｆｗ１～ｆｗ４）を１つの応答波形ｆｗに統合する。すなわち、統合部２５０は、演算結果取得部２４０が取得する複数の初期値ｆ０を統合する。
【0198】
［実施形態のまとめ］
以上説明したように、初期時刻t０が与えられた場合において、演算装置１０は、初期時刻ｔ０における応答波形ｆｗの初期値ｆ０を逆ラプラス変換法による離散的な時間応答解析によって算出する。
以上説明したように、初期時刻t０が与えられた場合において、演算装置１０は、初期時刻ｔ０における応答波形ｆｗの初期値ｆ０を逆ラプラス変換法による離散的な時間応答解析によって算出する。
【0199】
一般に、応答波形ｆｗの初期値ｆ０が求められれば、初期時刻ｔ０以降の応答波形ｆｗは、時間軸に沿った逐次演算（例えば、上述したＦＤＴＤ法）によって求めることができる。しかしながら、上述したＦＤＴＤ法によって、ある初期時刻ｔ０における初期値ｆ０を求める場合には、初期時刻ｔ０よりも前（過去）の時間軸に沿った逐次演算の結果が必要になる。つまり、従来の手法によると、任意の初期時刻ｔ０が与えられた場合において、この初期時刻ｔ０における応答波形ｆｗの値（つまり、初期値ｆ０）を求めることは困難である。
【0200】
本実施形態の演算装置１０は、ｓ領域において演算を行った結果を逆ラプラス変換によってｔ領域に逆変換することにより、離散的な時間応答解析を行う。このため、演算装置１０は、任意の初期時刻ｔ０が与えられた場合に、この初期時刻ｔ０における応答波形ｆｗの値（つまり、初期値ｆ０）を数値演算によって求めることができる。つまり、演算装置１０は、応答波形ｆｗの一部のみを数値演算によって求めることができる。
本実施形態の演算装置１０は、逆ラプラス変換法による離散的な時間応答解析によって算出するため、他の演算装置１０の演算結果を参照することなく初期値ｆ０を算出することができる。したがって、本実施形態の演算装置１０を複数台並列化した場合には、各演算装置１０が初期値ｆ０を他の演算装置１０と同時並列に算出することができる。つまり、本実施形態の演算装置１０によれば、応答波形ｆｗの初期値ｆ０の演算を同時並列化することができる。
本実施形態の演算装置１０は、逆ラプラス変換法による離散的な時間応答解析によって算出するため、他の演算装置１０の演算結果を参照することなく初期値ｆ０を算出することができる。したがって、本実施形態の演算装置１０を複数台並列化した場合には、各演算装置１０が初期値ｆ０を他の演算装置１０と同時並列に算出することができる。つまり、本実施形態の演算装置１０によれば、応答波形ｆｗの初期値ｆ０の演算を同時並列化することができる。
【0201】
また、本実施形態の演算装置１０は、逆ラプラス変換における数値演算についてＦＩＬＴ法を採用している。逆ラプラス変換は一般的には困難な演算であるが、本実施形態の演算装置１０によれば、逆ラプラス変換における数値演算を容易に行うことができる。
【0202】
また、本実施形態の演算装置１０によれば、演算装置１０の並列化を容易にすることができる。すなわち、本実施形態の演算装置１０は、上述したように、初期時刻ｔ０が与えられれば、この初期時刻ｔ０における初期値ｆ０（つまり、応答波形ｆｗの一部）を数値演算によって求めることができる。このため、複数の演算装置１０に対して、互いに異なる初期時刻ｔ０を与えることにより、各演算装置１０がこの演算装置１０に対する初期値ｆ０を算出することができる。例えば、演算装置１０－１～－４の４台の演算装置１０に対して、初期時刻ｔ０－１～－４を与えることにより、演算装置１０－１は、初期値ｆ０－１を、演算装置１０－２は初期値ｆ０－２を、演算装置１０－３は初期値ｆ０－３を、演算装置１０－４は初期値ｆ０－４をそれぞれ算出する。ここで、各演算装置１０は、逆ラプラス変換法による離散的な時間応答解析によって他の演算装置１０とは独立して初期値ｆ０を算出する。つまり、演算装置１０は、他の演算装置１０の演算結果を参照せずに初期値ｆ０を算出することができる。したがって、演算装置１０は、互いに他の演算装置１０が初期値ｆ０の演算中であったとしても、他の演算装置１０の演算状況によらず、初期値ｆ０を算出することができる。つまり、演算装置１０は、他の演算装置１０の演算状況によらず、個別に初期値ｆ０（すなわち、応答波形ｆｗの一部）を算出することができるため、演算装置１０の並列化を容易にすることができる。
【0203】
また、本実施形態の演算装置１０は、算出された初期値ｆ０に基づく時間軸に沿った逐次演算によって、応答波形ｆｗ（波形）についての所定の演算時間幅ｔｗぶんの演算結果ＲＳを算出する。このように構成することにより、演算装置１０は、初期値ｆ０のみならず、初期時刻ｔ０以降の時間幅の応答波形ｆｗをも算出することができる。また、演算装置１０は、初期時刻ｔ０以降の時間幅の応答波形ｆｗを、従来手法（例えば、ＦＤＴＤ法）によって算出するため、既存の演算ソフトウエアを流用することができ、演算ソフトウエアを流用できない場合に比べて簡便に演算ソフトウエアを構成することができる。
【0204】
また、本実施形態の演算装置１０において、演算時間幅ｔｗが、演算部１２０の演算性能ＣＰに基づいて定められている。したがって、応答波形ｆｗの算出を並列化した場合に、比較的高性能な演算装置１０にはより多くの演算時間幅ｔｗを割り当て、比較的低性能な演算装置１０にはより少ない演算時間幅ｔｗを割り当てることができる。ここで、高性能な演算装置１０は演算速度が相対的に速く、低性能な演算装置１０は演算速度が相対的に遅い。各演算装置１０は、演算性能ＣＰに基づく演算時間幅ｔｗを受け持つことにより、例えば、低性能な演算装置１０に対して多くの演算時間幅ｔｗが割り当てられることにより演算時間が長くなってしまうなど問題の発生を低減することができる。
【0205】
また、本実施形態の演算管理装置２０は、演算条件供給部２３０と、演算結果取得部２４０と、統合部２５０とを備えている。演算条件供給部２３０は、複数の演算装置１０に対して、演算装置１０ごとに互いに異ならせた初期時刻ｔ０を演算条件ＯＣとして供給する。演算結果取得部２４０は、供給された初期時刻ｔ０に基づいて、複数の演算装置１０がそれぞれ算出する初期値ｆ０を取得する。統合部２５０は、各演算装置１０が演算した複数の初期値ｆ０（すなわち、応答波形ｆｗの一部）を１つの応答波形ｆｗに統合する。この本実施形態の演算管理装置２０によれば、演算装置１０の並列化を容易にすることができる。
【0206】
また、本実施形態の演算管理装置２０は、演算時間幅算出部２２０を備えている。この演算時間幅算出部２２０は、それぞれの演算装置１０が備える演算部１２０の演算性能ＣＰに基づいて、初期時刻ｔ０からの演算時間幅ｔｗを演算装置１０ごとに算出する。この演算時間幅算出部２２０を備えることにより、応答波形ｆｗの算出を複数の演算装置１０によって並列化した場合に、比較的高性能な演算装置１０にはより多くの演算時間幅ｔｗを割り当て、比較的低性能な演算装置１０にはより少ない演算時間幅ｔｗを割り当てることができる。このように演算管理装置２０は、演算時間幅算出部２２０が各演算装置１０に割り当てる演算時間幅ｔｗを演算装置１０ごとに算出することにより、それぞれの演算装置１０による演算時間を演算装置１０ごとに調整することができる。例えば、演算管理装置２０は、低性能な演算装置１０に対して少ない演算時間幅ｔｗを割り当て、高性能な演算装置１０に対して多くの演算時間幅ｔｗを割り当てることにより、各演算装置１０による演算開始から演算終了までの時間（つまり、演算処理時間）を揃えることができる。
【0207】
以上説明したように、本実施形態の演算装置１０や演算管理装置２０によれば、応答波形ｆｗを求める演算を複数の演算装置１０によって同時並行に行うことが可能になり、１台の演算装置１０が応答波形ｆｗを求める場合に比べて、演算装置１０の１台あたりの演算負荷を低減することができる。このため本実施形態の演算装置１０や演算管理装置２０によれば、スーパーコンピュータのような極めて高性能な演算装置１０によらなくても妥当な演算処理時間によって応答波形ｆｗを算出することができる。
【0208】
［演算結果の例］
並列化された複数の演算装置１０によって算出される応答波形ｆｗの具体例について、図９を参照して説明する。
図９は、本実施形態の複数の演算装置１０によって算出される応答波形ｆｗの具体例を示す図である。この一例では、演算装置１０－１１～－１６（いずれも不図示）の６台の演算装置１０の並列演算によって応答波形ｆｗが算出される。演算時間幅ｔｗ１１（時刻ｔ１～時刻ｔ２）について、演算装置１０－１１が応答波形ｆｗ１１を算出する。演算時間幅ｔｗ１２（時刻ｔ２～時刻ｔ３）について、演算装置１０－１２が応答波形ｆｗ１２を算出する。以下同様にして、演算装置１０－１３～－１６が応答波形ｆｗ１３～１６をそれぞれ算出する。
演算管理装置２０の統合部２５０は、算出された応答波形ｆｗ１１～１６を応答波形ｆｗとして統合する。
並列化された複数の演算装置１０によって算出される応答波形ｆｗの具体例について、図９を参照して説明する。
図９は、本実施形態の複数の演算装置１０によって算出される応答波形ｆｗの具体例を示す図である。この一例では、演算装置１０－１１～－１６（いずれも不図示）の６台の演算装置１０の並列演算によって応答波形ｆｗが算出される。演算時間幅ｔｗ１１（時刻ｔ１～時刻ｔ２）について、演算装置１０－１１が応答波形ｆｗ１１を算出する。演算時間幅ｔｗ１２（時刻ｔ２～時刻ｔ３）について、演算装置１０－１２が応答波形ｆｗ１２を算出する。以下同様にして、演算装置１０－１３～－１６が応答波形ｆｗ１３～１６をそれぞれ算出する。
演算管理装置２０の統合部２５０は、算出された応答波形ｆｗ１１～１６を応答波形ｆｗとして統合する。
【0209】
［従来の演算結果との比較例］
図１０は、従来手法により１台の演算装置によって算出される応答波形ｗｏの具体例を示す図である。図１０の応答波形ｗｏは、１台の演算装置が従来手法（例えば、ＦＤＴＤ法）による時間軸に沿った逐次演算によって応答波形ｗｏを算出する。図１０に示す一例の場合、図９の応答波形ｆｗを求める際に与えられた演算対象情報ＭＤＬと同一の演算対象情報ＭＤＬが１台の演算装置に対して与えられている。つまり、図９の応答波形ｆｗと、図１０の応答波形ｗｏとは同一の演算対象情報ＭＤＬについての演算の結果である。図９に示した複数の演算装置１０による並列演算の結果である応答波形ｆｗと、１台の演算装置による応答波形ｗｏとがよく一致していることが示される。つまり、本実施形態の演算装置１０及び演算管理装置２０によれば、複数台の演算装置１０が同時並列演算を行ったとしても、１台の演算装置によって行う従来手法と同一の演算の結果が得られる。
図１０は、従来手法により１台の演算装置によって算出される応答波形ｗｏの具体例を示す図である。図１０の応答波形ｗｏは、１台の演算装置が従来手法（例えば、ＦＤＴＤ法）による時間軸に沿った逐次演算によって応答波形ｗｏを算出する。図１０に示す一例の場合、図９の応答波形ｆｗを求める際に与えられた演算対象情報ＭＤＬと同一の演算対象情報ＭＤＬが１台の演算装置に対して与えられている。つまり、図９の応答波形ｆｗと、図１０の応答波形ｗｏとは同一の演算対象情報ＭＤＬについての演算の結果である。図９に示した複数の演算装置１０による並列演算の結果である応答波形ｆｗと、１台の演算装置による応答波形ｗｏとがよく一致していることが示される。つまり、本実施形態の演算装置１０及び演算管理装置２０によれば、複数台の演算装置１０が同時並列演算を行ったとしても、１台の演算装置によって行う従来手法と同一の演算の結果が得られる。
【0210】
また、１台の演算装置が従来手法（例えば、ＦＤＴＤ法）による時間軸に沿った逐次演算によって応答波形ｗｏを算出するための演算処理時間は、９０分間であった。一方、これと同一の演算性能ＣＰの演算装置１０を利用した場合の、本実施形態の演算手法による実験結果は次の通りである。例えば、４台の演算装置１０による並列演算によって応答波形ｆｗを算出する場合、１台の演算装置１０が初期値ｆ０を算出するための演算処理時間が５分間であり、演算時間幅ｔｗぶんの応答波形ｆｗを算出するための演算処理時間が２８分間であった。つまり、本実施形態の演算手法によると、１台の演算装置が１０初期値ｆ０の算出と、演算時間幅ｔｗぶんの応答波形ｆｗの算出とに要する演算処理時間は、３２分間であった。本実施形態の演算装置１０は並列化されているため、応答波形ｆｗ全体を３２分間で算出することができる。したがって、この一例の場合において、本実施形態の演算手法によれば、従来手法の演算処理時間の約３５％の演算処理時間によって応答波形ｆｗを算出することができることが示された。
【0211】
［変形例］
なお、これまで演算装置１０は、初期値ｆ０と応答波形ｆｗとを算出するとして説明したが、これに限られない。演算装置１０は、初期値ｆ０のみを算出してもよい。演算装置１０は、演算対象情報ＭＤＬと初期時刻ｔ０とが与えられた場合、逆ラプラス変換法による離散的な時間応答解析によって初期値ｆ０を算出することができる。ここで、並列演算する演算装置１０の数は、初期値ｆ０の数に対応する。ｎ台の演算装置１０が初期値ｆ０－１～－ｎを算出し、算出した初期値ｆ０－１～－ｎをプロットしても応答波形ｆｗを得ることができる。初期値ｆ０－１～－ｎをプロットすることにより応答波形ｆｗを得れば、演算装置１０が応答波形ｆｗを算出しなくても、応答波形ｆｗを得ることができる。つまり、演算装置１０は、上述した逆ラプラス変換法による離散的な時間応答解析によって初期値ｆ０を算出すれば、算出された初期値ｆ０に基づく時間軸に沿った逐次演算を行わなくても応答波形ｆｗを得ることができる。
例えば、数万台の演算装置１０を並列化すれば、算出された数万点の初期値ｆ０を応答波形ｆｗとしてプロットすることができる。
なお、これまで演算装置１０は、初期値ｆ０と応答波形ｆｗとを算出するとして説明したが、これに限られない。演算装置１０は、初期値ｆ０のみを算出してもよい。演算装置１０は、演算対象情報ＭＤＬと初期時刻ｔ０とが与えられた場合、逆ラプラス変換法による離散的な時間応答解析によって初期値ｆ０を算出することができる。ここで、並列演算する演算装置１０の数は、初期値ｆ０の数に対応する。ｎ台の演算装置１０が初期値ｆ０－１～－ｎを算出し、算出した初期値ｆ０－１～－ｎをプロットしても応答波形ｆｗを得ることができる。初期値ｆ０－１～－ｎをプロットすることにより応答波形ｆｗを得れば、演算装置１０が応答波形ｆｗを算出しなくても、応答波形ｆｗを得ることができる。つまり、演算装置１０は、上述した逆ラプラス変換法による離散的な時間応答解析によって初期値ｆ０を算出すれば、算出された初期値ｆ０に基づく時間軸に沿った逐次演算を行わなくても応答波形ｆｗを得ることができる。
例えば、数万台の演算装置１０を並列化すれば、算出された数万点の初期値ｆ０を応答波形ｆｗとしてプロットすることができる。
【0212】
以上、本発明の実施形態を図面を参照して詳述してきたが、具体的な構成はこの実施形態に限られるものではなく、本発明の趣旨を逸脱しない範囲で適宜変更を加えることができる。上述した実施形態に記載の構成を組み合わせてもよい。
【0213】
なお、上記の実施形態における演算装置１０及び演算管理装置２０が備える各部は、専用のハードウェアにより実現されるものであってもよく、また、メモリおよびマイクロプロセッサにより実現させるものであってもよい。
【0214】
なお、演算装置１０及び演算管理装置２０が備える各部は、メモリおよびＣＰＵ（中央演算装置）により構成され、演算装置１０及び演算管理装置２０が備える各部の機能を実現するためのプログラムをメモリにロードして実行することによりその機能を実現させるものであってもよい。
【0215】
また、演算装置１０及び演算管理装置２０が備える各部の機能を実現するためのプログラムをコンピュータ読み取り可能な記録媒体に記録して、この記録媒体に記録されたプログラムをコンピュータシステムに読み込ませ、実行することにより、制御部が備える各部による処理を行ってもよい。なお、ここでいう「コンピュータシステム」とは、ＯＳや周辺機器等のハードウェアを含むものとする。
【0216】
また、「コンピュータシステム」は、ＷＷＷシステムを利用している場合であれば、ホームページ提供環境（あるいは表示環境）も含むものとする。
また、「コンピュータ読み取り可能な記録媒体」とは、フレキシブルディスク、光磁気ディスク、ＲＯＭ、ＣＤ－ＲＯＭ等の可搬媒体、コンピュータシステムに内蔵されるハードディスク等の記憶装置のことをいう。さらに「コンピュータ読み取り可能な記録媒体」とは、インターネット等のネットワークや電話回線等の通信回線を介してプログラムを送信する場合の通信線のように、短時間の間、動的にプログラムを保持するもの、その場合のサーバやクライアントとなるコンピュータシステム内部の揮発性メモリのように、一定時間プログラムを保持しているものも含むものとする。また上記プログラムは、前述した機能の一部を実現するためのものであってもよく、さらに前述した機能をコンピュータシステムにすでに記録されているプログラムとの組み合わせで実現できるものであってもよい。
また、「コンピュータ読み取り可能な記録媒体」とは、フレキシブルディスク、光磁気ディスク、ＲＯＭ、ＣＤ－ＲＯＭ等の可搬媒体、コンピュータシステムに内蔵されるハードディスク等の記憶装置のことをいう。さらに「コンピュータ読み取り可能な記録媒体」とは、インターネット等のネットワークや電話回線等の通信回線を介してプログラムを送信する場合の通信線のように、短時間の間、動的にプログラムを保持するもの、その場合のサーバやクライアントとなるコンピュータシステム内部の揮発性メモリのように、一定時間プログラムを保持しているものも含むものとする。また上記プログラムは、前述した機能の一部を実現するためのものであってもよく、さらに前述した機能をコンピュータシステムにすでに記録されているプログラムとの組み合わせで実現できるものであってもよい。
【符号の説明】
【0217】
１…演算システム、１０…演算装置、１１０…演算条件取得部、１２０…演算部、１３０…出力部、２０…演算管理装置、２１０…取得部、２２０…演算時間幅算出部、２３０…演算条件供給部、２４０…演算結果取得部、２５０…統合部、３０…記憶装置、ＭＤＬ…演算対象情報、ｔｗ…演算時間幅、ｆ０…初期値、ｔ０…初期時刻、ｆｗ…応答波形
【図1】
【図2】
【図3】
【図4】
【図5】
【図6】
【図7】
【図8】
【図9】
【図10】