演算子生成方法、演算子生成装置、及びシミュレーション装置
【課題】時間領域有限差分法により行う弾性体からなる有限な媒質中の地震波の伝播についての数値解析において、計算効率、精度に加え安定性を好適に実現することができる演算子を生成する方法、演算子生成装置、及びシミュレーション装置を提供する。
【解決手段】状態ベクトルに作用し各離散点での変位ベクトルの空間１階差分演算子をＰｉ、弾性定数行列をＣとして、剛性行列Ｈ′をＰｉＴＣＰｉ、又はその線形結合の形に設計する。また、質量行列Ｔ′に関しては、離散点での差分近似により導かれる正定値の対角質量行列をＴ、δＴｎｕｍ＝Ｔ′−Ｔとして、回転及び並進に対応する前記状態ベクトルに演算子δＴｎｕｍを作用させた結果が０となるように設計する。
【解決手段】状態ベクトルに作用し各離散点での変位ベクトルの空間１階差分演算子をＰｉ、弾性定数行列をＣとして、剛性行列Ｈ′をＰｉＴＣＰｉ、又はその線形結合の形に設計する。また、質量行列Ｔ′に関しては、離散点での差分近似により導かれる正定値の対角質量行列をＴ、δＴｎｕｍ＝Ｔ′−Ｔとして、回転及び並進に対応する前記状態ベクトルに演算子δＴｎｕｍを作用させた結果が０となるように設計する。
【発明の詳細な説明】
【技術分野】
【０００１】
本発明は、時間領域有限差分法により行う、弾性体からなる有限な媒質中の波動についての数値解析に関し、その演算子生成方法、演算子生成装置、及びシミュレーション装置に関する。
【背景技術】
【０００２】
弾性体の運動方程式を解く手法として時間領域差分法が用いられている。時間領域差分法を用いた数値解析では、安定性、計算効率、精度の３つが問題となる。古くから数多くの時間領域差分法による弾性波動計算手法は開発されてきたが、これら安定性、計算効率、精度を同時に好適に実現する手法はこれまでなかった。
【０００３】
例えば、非特許文献１に示される時間領域差分計算手法は、弾性波動を効率よく、高精度に計算する手法を提供したが、安定性は保証されていなかった。
【０００４】
ここでいう、「効率がよい」とは、特定の計算精度を得るための数値解析を行う計算機の使用リソースが少ないことを意味する。具体的には、計算時間が少なく、また当該計算に用いる差分演算子の幅が抑制される手法は、高い効率性であると評価される。また、計算手法における「精度」の評価に関しては、許容される効率性の範囲で実現できる「最適な精度」というものを考えることができる。以下、「最適精度を持つ演算子」とは非特許文献２に示される条件（後述の（９）式）を満たす演算子を意味する。
【非特許文献１】Takeuchi,N., and R.J. Geller, 2000, Physics of the Earth and Planetary Interiors, 119,99-131.
【非特許文献２】Geller,R.J., and N. Takeuchi, 1995, Geophysical Journal International, 123, 449-470.
【非特許文献３】Geller,R.J., and N. Takeuchi, 1998, Geophysical Journal International, 135, 48-62.
【発明の開示】
【発明が解決しようとする課題】
【０００５】
既に述べたように、時間領域有限差分法により行う弾性体からなる有限な媒質中の波動についての数値解析において、従来の計算手法は、安定性、計算効率、精度を同時に好適に実現することができないという問題がある。
【０００６】
本発明は上記問題点を解決するためになされたものであり、時間領域有限差分法により行う弾性体からなる有限な媒質中の波動についての数値解析において、安定性、計算効率、精度を同時に好適に実現することができる演算子を生成する方法、演算子生成装置、及びシミュレーション装置を提供することを目的とする。
【課題を解決するための手段】
【０００７】
本発明に係る演算子生成方法及び演算子生成装置は、予測子修正子法を用いた時間領域有限差分法により行う、弾性体からなる有限な媒質中の波動についての数値解析に関し、運動方程式の陰形式の差分近似である次式、
【数２２】
（ここで、Δｔは時間方向の差分ステップを表し、ｃ＋，ｃ０，ｃ−は解析対象空間に設定した複数の離散点での変位ベクトルの各成分からなる状態ベクトルであって、それぞれ時刻ｔ＋Δｔ，ｔ，ｔ−Δｔにおけるものを表し、ｆは前記各離散点での外力ベクトルを表し、Ｔ′は質量行列を表し、Ｈ′は剛性行列を表す。）に基づいて設定される予測子及び修正子にて用いる演算子を生成するものであって、当該演算子として前記Ｈ′、及び次式、
δＴｎｕｍ＝Ｔ′−Ｔ
で表されるδＴｎｕｍ（ここで、Ｔは前記離散点での差分近似により導かれる正定値の対角質量行列である。）を生成する方法又は手段を含み、前記Ｈ′を、前記状態ベクトルに作用し前記各離散点での前記変位ベクトルの空間１階差分演算子Ｐｉ（ここで、ｉは、Ｉを前記解析対象空間の次元に応じて決まる所定の自然数として、１≦ｉ≦Ｉなる整数である。）と、弾性定数行列Ｃとから次式、
【数２３】
（ここで、ａｉは正のスカラー係数であり、演算子ＰｉＴは演算子Ｐｉを表す行列の転置行列で表される演算子である。）に基づいて生成し、かつ、回転及び並進に対応する前記状態ベクトルに前記演算子δＴｎｕｍを作用させた結果が０となるように当該δＴｎｕｍを設定する。
【０００８】
本発明に係るシミュレーション装置は、予測子修正子法を用いた時間領域有限差分法により行う、弾性体からなる有限な媒質中の波動についての数値解析を行うものであって、運動方程式の陰形式の差分近似である次式、
【数２４】
（ここで、Δｔは時間方向の差分ステップを表し、ｃ＋，ｃ０，ｃ−は解析対象空間に設定した複数の離散点での変位ベクトルの各成分からなる状態ベクトルであって、それぞれ時刻ｔ＋Δｔ，ｔ，ｔ−Δｔにおけるものを表し、ｆは前記各離散点での外力ベクトルを表し、Ｔ′は質量行列を表し、Ｈ′は剛性行列を表す。）に基づいて設定される予測子及び修正子であって、前記離散点での差分近似により導かれる正定値の対角質量行列Ｔ及び、δＴｎｕｍ＝Ｔ′−Ｔで表される演算子δＴｎｕｍを用いて、次式、
【数２５】
で与えられる前記予測子と、次式、
【数２６】
で与えられる前記修正子とをそれぞれ計算して、次式、
【数２７】
で与えられる前記状態ベクトルｃ＋を時間発展的に順次生成する演算手段を有し、前記Ｈ′は、前記状態ベクトルに作用し前記各離散点での前記変位ベクトルの空間１階差分演算子Ｐｉ（ここで、ｉは、Ｉを前記解析対象空間の次元に応じて決まる所定の自然数として、１≦ｉ≦Ｉなる整数である。）と、弾性定数行列Ｃとから次式、
【数２８】
（ここで、ａｉは正のスカラー係数であり、演算子ＰｉＴは演算子Ｐｉを表す行列の転置行列で表される演算子である。）で表される演算子であり、前記演算子δＴｎｕｍは、回転及び並進に対応する前記状態ベクトルに当該演算子δＴｎｕｍを作用させた結果が０となる演算子である。
【発明の効果】
【０００９】
本発明によれば、時間領域有限差分法により行う弾性体からなる有限な媒質中の波動についての数値解析に関し、安定性、計算効率、精度を同時に好適に実現することができる演算子及びシミュレーション装置が得られる。
【発明を実施するための最良の形態】
【００１０】
以下、本発明の実施の形態（以下実施形態という）について説明する。
【００１１】
本発明の実施形態である計算機システムは、演算処理部、記憶部、入力部及び出力部を含んで構成され、例えば、地震波解析に用いられる。演算処理部は、ＣＰＵ（Central Processing Unit：中央演算処理部）として、各種のプログラムを実行する。具体的には、演算処理部は、後述する演算子生成方法に基づく演算子生成処理を行ったり、当該処理で生成される演算子を用いたシミュレーションを行う。これにより、計算機システムは、本発明に係る演算子生成装置、シミュレーション装置を実現する。また、演算処理部は、計算機システムの各部を制御する。
【００１２】
記憶部は、演算処理部で実行される各種プログラムを格納する。また記憶部は、プログラム実行のための作業領域として使用され、計算過程でのデータを一時的に保持するためにも用いられる。また、記憶部は、演算子生成処理の結果として得られる演算子を保持したり、時間発展的なシミュレーションで得られる各時間ステップの解析結果を保持する。
【００１３】
入力部は、キーボードやマウス等を含み、ユーザは入力部を操作して、演算処理部において演算子生成処理やシミュレーションに用いられる各種パラメータ等の設定を行う。
【００１４】
出力部は、プリンタやディスプレイ等を含み、例えば、生成された演算子をプリントアウトしたり、シミュレーション結果をグラフィック表示したりすることができる。
【００１５】
本実施形態に係る演算子生成方法、演算子生成装置及びシミュレーション装置は、時間領域有限差分法により行う、弾性体からなる有限な媒質中の波動についての数値解析に関する。
【００１６】
ここで、演算子生成方法、演算子生成装置及びシミュレーション装置は相互に関連している。すなわち、演算子生成方法は、シミュレーション装置にて実行される時間領域有限差分法による数値解析処理に対応した演算子を生成する方法であり、演算子生成装置は当該演算子生成方法を実行して演算子を生成する。そして、シミュレーション装置は生成された演算子を用いて数値解析を行う。そのため、本発明に係る演算子生成及びシミュレーションの内容を理解する上では、それら全体を１つの計算スキームとして説明することが好適である。そこで、以下、本発明に係る演算子生成及びシミュレーションの内容を含んだ弾性波動解析スキームについて説明する。
【００１７】
［運動方程式］
本発明に関するシミュレーションの対象となる弾性体中の弾性波動伝播に関する運動方程式は以下のような式で表される。
【数２９】
【００１８】
ここで、ρは密度、Ｃｉｊｋｌは弾性定数、ｕｉは変位ベクトルであり、各添え字ｉ，ｊ，ｋ，ｌはそれぞれｘ，ｙ又はｚであり、三次元空間におけるｘｙｚ直交座標系の各軸方向の成分を表す。また、ｕｉ，ｊ はｕｉ のｊ方向の微分（ｊ＝ｘ，ｙ，ｚ）を表す。ｆｉは外力ベクトルのｉ方向成分である。なお、この式は、アインシュタインの縮約規則を用いて表現している。ちなみに当該規則は以下の数式表現においても適宜採用する。
【００１９】
［１ステップの陰形式の差分近似］
時間領域差分法による解析では運動方程式を最適計算スキームに離散化する。離散化により、上記運動方程式の陰形式の差分近似である次式が得られる。
【数３０】
【００２０】
ここで、Δｔは時間ステップを表す。ｃ＋，ｃ０，ｃ−は解析対象空間に設定した複数の離散点に対する変位ベクトルの各成分からなる状態ベクトルであって、それぞれ、時刻ｔ＋Δｔ，ｔ，ｔ−Δｔにおけるものを表す。また、ｆは前記各離散点での外力ベクトル、Ｔ′は質量行列、Ｈ′は剛性行列を表す。
【００２１】
この離散化手法が時間発展に対して完全に安定であるためには、Ｔ′が正定値、Ｈ′が半正定値であることが必要条件となる。ここで、Ｈ′は以下の形式で表現できれば、半正定値性の条件を満たす。
【数３１】
【００２２】
ここで、Ｐは運動方程式中に現れる空間１階微分の差分近似であり、上付き文字Ｔ は行列の転置を表す。また、Ｃは弾性定数行列（空間の不均質性及び異方性も含む）である。また、（３）式右辺に現れる半正定値演算子の線形結合である次式で表されるＨ′も半正定値性の条件を満たす。
【数３２】
【００２３】
ここで、ｉ = １，…，Ｉ（Ｉは整数）で、ａｉは正の係数で、Ｐｉはそれぞれ１階差分演算子である。
【００２４】
［予測子修正子法］
上記陰形式の方程式はそのままの形では計算機システムを用いた処理に適さず、実際には予測子修正子法を用いて計算機システムによるシミュレーションが行われる。予測子修正子法では、（２）式に基づいて、それぞれ陽形式で表現される予測子及び修正子が設定される。予測子及び修正子としてそれぞれ下記（５），（６）式を定義する。
【数３３】
【数３４】
【００２５】
ここで、Ｔは対角質量行列である。また、δＴｎｕｍは、
【数３５】
となる質量行列の残差である。
【００２６】
（５）式左辺括弧内の第１項は、時刻ｔ＋Δｔにおける状態ベクトルｃ＋の予測値である。シミュレーションでは、まず、予測子を用いて既知の状態ベクトルｃ０，ｃ−からｃ＋の予測値を計算する。次いで、修正子を用いて、予測値に対する修正値δｃ＋を計算する。そして、次式に示すように予測子と修正子とから得られた結果を加算して、次のタイムステップの状態ベクトルｃ＋が得られる。
【数３６】
【００２７】
［安定性の実現］
この予測子修正子法の安定性を保証するために、δＴｎｕｍ＝Ｔ′−Ｔを、以下の条件を満たすように設計する。すなわち、修正子で用いられる演算子δＴｎｕｍは、それをもとの運動方程式の固有周波数が０となる固有ベクトル（並進および回転に対応する）に作用させたとき、その結果が０となるという条件と、Ｔ′が正定値となるという条件とを満たすように設計する。
【００２８】
本発明に係る演算子生成方法、演算子生成装置は、このδＴｎｕｍに関する条件と、上述の（３）式又は（４）式の条件とを満たすように、（５），（６）式で用いる演算子Ｈ′，δＴｎｕｍを生成する。そして、予測子修正子法によるシミュレーションは、このように生成された演算子Ｈ′，δＴｎｕｍを用いて行われ、これにより安定性を確保することができる。
【００２９】
［最適精度の実現］
以上、シミュレーションにおいて安定性を実現する演算子の生成について述べてきたが、次に安定性に加え、さらに上述の「最適な精度」が得られる演算子の生成について説明する。
【００３０】
厳密な演算子の表現となる質量行列及び剛性行列をそれぞれＴ０，Ｈ０と定義し、演算子の残差を以下のように定義する。
δＴ＝Ｔ′−Ｔ０
δＨ＝Ｈ′−Ｈ０
【００３１】
また、正確な運動方程式の固有ベクトルをｕｍ（ここでｍは複数存在し得る固有モードの識別子である。）、当該ｕｍに対応する固有周波数をωｍとして、固有ベクトルに対するδＴ及びδＨの行列要素の対角成分をδＴｍｍ及びδＨｍｍと定義する。
【００３２】
Ｔ′及びＨ′が非特許文献２の最適精度演算子となる条件は、次式が空間についての差分ステップの２次以上の精度で成立することである。
【数３７】
【００３３】
この条件を満たす演算子を得るために、離散化した運動方程式に現れる１階差分演算子を以下のように設計する。Ｄｘ，Ｄｙ，Ｄｚを所定の微分定義位置でのそれぞれｘ，ｙ，ｚ軸方向に対する２点１階差分演算子とする。また、Ｅｘｒ，Ｅｙｓ，Ｅｚｔを前記微分定義位置でのｘ軸方向、ｙ軸方向及びｚ軸方向それぞれに関する３点差分近似による、最適精度を持つように設計した２次精度０階差分演算子とする。ここでｒ，ｓ，ｔはそれぞれ１又は２で、３点差分近似に用いる３つの離散点の微分定義位置に対する配置の相違を示すインデックスである。例えば、ｒ＝１は３つの離散点がｘ軸の負方向にシフトした配置、ｒ＝２は正方向にシフトした配置を表す。ここで、各離散点のｘ，ｙ，ｚ各方向における中点の座標をそれぞれｘｃ，ｙｃ，ｚｃとし、離散点の間隔をｈとする。例えば、微分定義位置の座標を（ｘｃ，ｙｃ，ｚｃ）とした場合、ｒ＝１は具体的には、３つの離散点のうち中央の点のｘ座標がｘｃ−ｈ/２となる配置を示し、ｒ＝２は中央の点のｘ座標がｘｃ＋ｈ/２となる配置を示す。他のインデックスｓ，ｔについても同様である。
【００３４】
これらＤｘ，Ｄｙ，Ｄｚ、及びＥｘｒ，Ｅｙｓ，Ｅｚｔを用いて、１階差分演算子Ｄｘｒｓｔ，Ｄｙｒｓｔ，Ｄｚｒｓｔを次式で定義する。
【数３８】
【００３５】
ちなみに、例えば、差分演算子Ｄｘｒｓｔは、Ｄｘを微分方向であるｘ軸とは別の他の２つの軸、すなわちｙ軸及びｚ軸両方向、に拡散した１階微分の差分近似である。Ｄｙｒｓｔ，Ｄｚｒｓｔについても同様である。
【００３６】
例えば、上式中に現れるＥｘ１，Ｅｘ２の具体形を、格子点上の任意連続関数ｆに作用する形で示したものが次に示す（１３），（１４）式それぞれの第１行、第２行の式である。Ｅｙｓ，Ｅｚｔも同様の具体形を有する。
【数３９】
【００３７】
なお、（１３），（１４）式それぞれの第３行の式は第１行、第２行で表される式をテイラー展開した結果を示している。このようにＥｘｒ，Ｅｙｓ，Ｅｚｔは、０階微分（当該テイラー展開の第１項）に特定の係数を有する２階微分成分（当該テイラー展開の第２項）が付加されたものとなり、これにより（９）式で示した最適精度の条件が満たされるようになる。
【００３８】
最適精度条件が満たされる場合には、（４）式の差分演算子Ｐｉは、Ｄｘｒｓｔ，Ｄｙｒｓｔ，Ｄｚｒｓｔを用いてｒ，ｓ，ｔの８通り（２３通り）の組み合わせのそれぞれに対応して定められる。すなわち、（４）式の説明で述べたＩは８であり、ｉは１≦ｉ≦Ｉなる整数である。また、各ａｉは１／８となる。
【００３９】
具体的には、Ｐｒｓｔを、任意なベクトルｖ，ｗを用いた以下の双二次形式で定義する。
【数４０】
【００４０】
ここでＰｒｓｔはＤｉｒｓｔ（ｉ＝ｘ，ｙ，ｚ）の組み合わせを行列表記した１階差分演算子であり、ｒ，ｓ，ｔの組み合わせに応じて番号付けを行い、各Ｐｉと対応させる。
【００４１】
２次元(Ｐ−ＳＶ問題およびＳＨ問題) の場合（解析対象空間をｘｚ平面とする）は、３次元の場合のＤｘｒｓｔ，Ｄｙｒｓｔ，Ｄｚｒｓｔに相当する差分演算子Ｄｘｒｓ，Ｄｚｒｓが次式で定義される。ここで、ｒとｓは上述の３次元の場合のｒとｔに対応するインデックスであり、それぞれの値は１又は２である。
【数４１】
【００４２】
ちなみに、例えば、差分演算子Ｄｘｒｓは、Ｄｘを微分方向であるｘ軸とは別の他の１つの軸、すなわちｚ軸方向に拡散した１階微分の差分近似である。Ｄｚｒｓについても同様である。
【００４３】
最適精度条件が満たされる場合には、（４）式の演算子Ｐｉは、Ｄｘｒｓ，Ｄｚｒｓを用いてｒ，ｓの４通り（２２通り）の組み合わせのそれぞれに対応して定められる。すなわち、２次元の場合、１≦ｉ≦Ｉなる整数ｉの上限値Ｉは４であり、また、各ａｉは１／４となる。
【００４４】
具体的には、２次元の場合のＰｒｓを以下の表式を用いて定義する。
【数４２】
【００４５】
ここでＰｒｓはＤｉｒｓ（ｉ＝ｘ，ｚ）の組み合わせを行列表記した１階差分演算子であり、ｒ，ｓの組で番号付けを行うことで各Ｐｉに対応させる。
【００４６】
この演算子が作用する状態ベクトルとして、ｙ成分の変位のみを持つようにとることでＳＨ問題は定義され、ｘ，ｚ成分のみを持つようにとることでＰ−ＳＶ問題は定義される。
【００４７】
上記のようにＨ′を複数のＰｉの和で表現することにより、数値演算子の一般固有値問題を解いたとき、並進・回転モード以外の固有周波数が０となる固有モードの存在を防ぐことができる。
【００４８】
［非物理現象の抑制］
上記のように定義されたＨ′を用いてシミュレーションを行うと、３次元、２次元のいずれの解析においても、空間一格子ごとに振動するような成分を有する高波数の波が低い振動数を持つようになる現象が生じ得る。物理的にいえば、高波数（短波長）の波はより高い振動数を持つのが自然であり、上述の現象は非物理的な現象である。この現象の発生を抑制する演算子について以下説明する。
【００４９】
非物理現象の発生を抑制するためには、Ｆｉｒで定義されるｉ軸（ｉ＝ｘ，ｙ，ｚ）方向の３点２階差分演算子を用いて演算子Ｑを設計し、上述のＨ′を以下のように修正する。ただし、ｒは、上述のＥｘｒの定義と同様に、差分演算子を定義する離散点の配置のシフト方向を示すインデックスであり、その値は１又は２である。
【数４３】
【００５０】
ここで、ｊ＝１，…，Ｊ（Ｊは問題の次元により定義される所定の整数）であり、εｊは正の係数である。また、Ｑｊは各ｊに対応する２階差分演算子である。（１９）式は２次元、３次元問題共通である。具体的なＦｉｒ，Ｑｊ，εｊの表式は後に記載する。
【００５１】
ここまでで、安定でかつ最適精度の条件を満たし、非物理現象の励起を抑制する剛性行列Ｈ′の定義が完了した。
【００５２】
［質量行列の定義］
次に、安定でかつ最適精度条件に適合し得る正定値質量行列Ｔ′を示す。最適精度条件に適合し得る質量行列Ｔ′を、予測子計算に用いる対角質量行列Ｔと、残差δＴｎｕｍとを用いて（７）式と同じ形式の次式で表す。
【数４４】
【００５３】
前述のように、予測子修正子法の安定性のために、δＴｎｕｍは、離散化された剛性行列Ｈ′と同様に並進及び回転モードに作用させたとき、その値が０となる必要がある。この条件を満たしつつ、最適精度の条件も満たすために、Ｈ′の表式を利用して、以下のようにδＴｎｕｍを定義する。
【数４５】
【００５４】
この式の右辺は剛性行列の表現における弾性行列ＣをＸで置き換えたものである。Ｘは、最適精度条件を満たすように、密度および格子間隔Δｘａ（ａはｘ，ｙ，ｚのいずれか）を用いて定義される値で、
【数４６】
となるものである。ここで、ρ(α) はα番目の微分定義点における密度の値であり、δａｂ，δｃｄ，δａｃはクロネッカーのデルタであり、ａ，ｂ，ｃ，ｄはそれぞれ座標軸ｘ，ｙ，ｚに対応する添え字である。
【００５５】
ここまでで、予測子修正子法を用いた手法を安定にかつ、最適な精度で計算するための、剛性行列Ｈ′と、質量行列Ｔ′（及びδＴｎｕｍ）を導出することができた。
【００５６】
これらの演算子を用い、予測子修正子法を用いて時間発展計算を行うことで、時刻格子、空間格子上における任意の点における状態ベクトルで表せる変位場を計算することが可能となる。
【００５７】
［シミュレーションにおける計算効率の向上］
予測子修正子法におけるほとんどの計算負荷は、Ｔ′，Ｈ′による空間差分演算部分になる。この空間差分演算を以下の手順で行うことにより、シミュレーション実行時の計算効率を向上させることができる。
【００５８】
上述したように、Ｈ′，Ｔ′は（３）式のような形式の行列の積を含んでいる。例えば、（３）式で表される演算子を状態ベクトルに作用させる演算では、
ｓｔｅｐ１：差分演算子行列Ｐと格子点における変位の値に対応する状態ベクトルｕとの積、
ｓｔｅｐ２：step１の差分演算結果の各要素を弾性定数倍する演算、
ｓｔｅｐ３：step２の演算結果を、差分演算子の転置であるＰＴ倍する演算、
を連続的に、効率よく行う必要がある。なお、弾性定数は微分定義位置において定義されているものとする。
【００５９】
３次元上の格子点上の変位を、通常の有限要素法で用いられるように番号付けを行い、１次元ベクトルに置き換えたものをｕｊとし、Ｐの要素をＰｉｊとすると、上記step１の演算は、
【数４７】
と書き表される。ｕｉは格子点での値で定義されるベクトルであるのに対し、ｖｉは微分定義位置における演算結果を１次元的に番号付けして整列したベクトルである。（２３）式は、或る微分定義位置(i) におけるｖｉを計算するために微分定義位置ｉ周辺の各ｕｊの値にＰｉｊで定義される重みを乗算し、和をとる形になっている。図１は、２次元の場合の（２３）式の演算を説明する図であり、ステンシルの模式図である。図１において、白丸（○）が格子点であり、黒四角（■）が微分定義位置を表す。各格子点の近くに示す数値はＰｉｊを示す。なお、一般的には微分定義位置に対して図１（ａ）のように６点の格子点が存在し得るが、微分定義位置が解析対象空間の境界に隣接する場合、図１（ｂ）のように４点からなるステンシルが設定され得る。
【００６０】
ｓｔｅｐ２の演算は、単純なスカラー倍で書き表されるので特に工夫は必要でない。
【００６１】
ｓｔｅｐ３の演算は、
【数４８】
と書き表される。（２４）式において、ｗｋは格子点における演算結果を表すベクトルである。また（２４）式において、ｗｋ（ｊ）＝Ｐｊｋｖｊとおいた。ここで、ｗｋ（ｊ）はｖｊで表されるベクトルのｊ番目の微分定義点の値を、ｋ番目の格子点に分配することで計算される。すなわち、右辺のｖｉが上述のように微分定義位置での値で定義されるベクトルであるのに対し、ｗｋはｕｉと同じく格子点における演算結果を表すベクトルである。図２は、２次元の場合の（２４）式の演算を説明する図であり、ステンシルの模式図である。図２は図１と同様、白丸（○）が格子点であり、黒四角（■）が微分定義位置を表す。また、各格子点の近くに示す数値はＰｉｊを示す。なお、図２（ａ），（ｂ）でのステンシルの違いも図１と同様である。
【００６２】
すべての微分定義位置について、上記step１〜３の演算を繰り返し行うことで、全ての行列演算が完了する。この計算方法によれば巨大な行列とベクトルの積の計算を部分的な演算の和で効率よく実行することができる。
【００６３】
［Ｑの具体的な表現］
以下に上述のＦｉｒ，Ｑｊ及びεｊの具体形について表記する。まず３次元の場合について記述し、その後で２次元の場合について記述する。
【００６４】
はじめに、２階差分演算子Ｆｉｒを０階差分演算子を定義したときと同様に、連続関数ｆに対する離散化表現を用いて定義する。以下にＦｘ１の場合の具体形を示す。
【数４９】
【００６５】
他のＦｉｒについても同様に定義される。これらのＦｉｒの定義を用い、（１５）式により表されるＰｒｓｔの定義式におけるＤｉｒｓｔを以下の６パターンで置き換え、Ｐｒｓｔに相当する部分を抽出することによってＱｊｒｓｔ を定義する。
【数５０】
【００６６】
これらの６パターンの置き換えにより生成される２階差分演算子Ｑｊｒｓｔを用いて、以下のように剛性行列の要素Ｈｒｓｔを定義する。
【数５１】
【００６７】
そして、最終的に剛性行列Ｈ′を
【数５２】
で計算する。上式ではｊ＝１，２，３に対してεｊ＝５/１４４，ｊ＝４，５，６に対してεｊ＝ (５/１４４)２を用いた。５/１４４以上の値を使用した手法も可能であるが、その場合時間ステップΔt をより小さくとる必要があるため、５/１４４を用いるのが好適である。
【００６８】
２次元の場合の２階差分演算子Ｑｊｒｓの表現は３次元の場合と同様に（１８）式により表されるＰｒｓの定義式におけるＤｉｒｓを以下のように置き換えたものに対して、Ｐｒｓに相当する演算子を抽出することにより計算することができる。
【数５３】
【００６９】
これらを用いて構成される２階差分演算子Ｑｊｒｓを用いて、２次元の場合の剛性行列の要素Ｈｒｓを定義する。
【数５４】
【００７０】
上式ではεｊ＝５/１４４，（ｊ＝１，２）を用いるが、２次元の場合は４/１４４ 以上であれば、非物理的な現象を抑制することができる。これらの剛性行列の要素Ｈｒｓを用いて最終的に剛性行列Ｈ′は
【数５５】
と書き表される。
【００７１】
［上述の計算スキームの応用例等］
本計算スキームでは、各微分定義点に弾性定数(異方性を含む)および密度の値を与えることで、全ての行列を定義する。任意不均質構造は、弾性定数、密度の値を各微分定義点に与えることで実現でき、また、任意の表面地形形状は、空中の弾性定数および密度に０を与えることで導入できる。
【００７２】
また、本計算スキームは陽解法を用いて時間を更新していく手法であり、時間ステップΔｔのとりうる最大値は、固有値解析により求められる最大固有値により決定することが可能である。これは一般的な手法であるので、ここでは特に記述しない。
【００７３】
上述の演算子生成方法により定義される安定かつ最適精度を持つ空間差分演算子は、既存の吸収境界条件を用いる手法と組み合わせて使用することもできる。よって、このような手法も本発明の演算子生成方法・装置及びシミュレーション装置の範囲に含まれる。
【００７４】
また、上記により定義される安定かつ最適精度を持つ空間差分演算子は、既存の非弾性減衰を表現する手法と組み合わせて使用することができる。このような手法も本発明の範囲に含まれる。
【００７５】
本発明においては３次元および２次元問題についての手法についてのみ扱う。これは、１次元の場合は、非特許文献３によって定義される手法と等価であるためである。
【００７６】
「時間領域差分法」には、(a)「規則的な格子上で離散化し、変位もしくは速度、加速度を変数とし、時間発展を計算する手法」と、(b)「半格子幅での食い違いをもつ２種類の格子上で離散化し、変位もしくは速度、加速度、および応力を変数とし、時間発展を計算する手法」とが含まれる。ここまでは(a) について説明してきたが、本手法で用いた空間差分演算子を(b) の手法に適用することが可能である。すなわち、変位および速度の微分にＰ（もしくはＰＴ ）を用い、応力の微分にはＰＴ（もしくはＰ）を用いて、時間発展計算を行うことで、安定性が保障されることになる。このような手法の場合、非特許文献３で示されているように、予測子修正子法は適用できないため、高精度性は保障できないが、予測子のみを用いて時間発展計算を行うことで、安定化の実現が可能である。
【００７７】
なお、一般に解析対象空間に設定される格子点の数は非常に大きく、例えば、行列形式で表される離散化した演算子の規模も巨大となる。そのため、当該巨大な演算子行列の全体を生成して記憶部に保持したり、記憶部に保持した巨大な行列と、巨大な状態ベクトルとを直接、乗算することは、計算効率の観点から得策ではない。この点については、既に１つの工夫を、図１及び図２を用いて説明した。このように、本発明に係るシミュレーション装置では、状態ベクトルに演算子行列を作用させる計算を、解析対象空間の一部分ずつ行い得る。そして、その場合には、演算子行列のうち、これから計算を行う一部分に対応した行列要素が記憶部上に存在すれば処理を実行することができる。本発明に係る演算子生成方法、演算子生成装置は、このように必要とする一部分ずつを逐次的にシミュレーション装置に提供する演算子生成の仕方を含んでいる。
【００７８】
また、演算子Ｈ′，δＴｎｕｍは上述のように演算子（行列）の積を含む形式で表され得る。この場合に、演算子を生成するとは、その行列の積を計算した結果得られる行列Ｈ′や行列δＴｎｕｍの行列要素を算出することには限定されず、積を構成する個々の演算子行列の要素を求めることも含む。例えば、図１及び図２を用いて説明した処理方法では、本発明の演算子生成では、Ｈ′を構成する行列の積ＰｉＴＣＰｉの１階差分演算子ＰｉＴ，Ｐｉの行列要素は直接的・具体的に定めるが、Ｈ′の行列要素の具体的な値までは計算しない。しかし、この場合においてもＰｉＴ，Ｐｉに基づいてＨ′は（４）式により定義され、この意味においてＨ′は生成されているということができる。
【産業上の利用可能性】
【００７９】
本発明に係る演算子生成方法、演算子生成装置、及びシミュレーション装置は、地震波の解析に用いることができる。地球内部を伝播してきた地震波形記録が持つ情報を十分に活用することができれば、より詳細で正確な内部構造推定が可能となる。高精度かつ高解像度の地球内部構造モデルの推定は、例えば、防災や、油田開発における物理探査などに利用することができる。
【図面の簡単な説明】
【００８０】
【図１】ＰＴＣＰの形式の行列の積と状態ベクトルｕとの積の演算方法を２次元の場合について説明するステンシルの模式図である。
【図２】ＰＴＣＰの形式の行列の積と状態ベクトルｕとの積の演算方法を２次元の場合について説明するステンシルの模式図である。
【技術分野】
【０００１】
本発明は、時間領域有限差分法により行う、弾性体からなる有限な媒質中の波動についての数値解析に関し、その演算子生成方法、演算子生成装置、及びシミュレーション装置に関する。
【背景技術】
【０００２】
弾性体の運動方程式を解く手法として時間領域差分法が用いられている。時間領域差分法を用いた数値解析では、安定性、計算効率、精度の３つが問題となる。古くから数多くの時間領域差分法による弾性波動計算手法は開発されてきたが、これら安定性、計算効率、精度を同時に好適に実現する手法はこれまでなかった。
【０００３】
例えば、非特許文献１に示される時間領域差分計算手法は、弾性波動を効率よく、高精度に計算する手法を提供したが、安定性は保証されていなかった。
【０００４】
ここでいう、「効率がよい」とは、特定の計算精度を得るための数値解析を行う計算機の使用リソースが少ないことを意味する。具体的には、計算時間が少なく、また当該計算に用いる差分演算子の幅が抑制される手法は、高い効率性であると評価される。また、計算手法における「精度」の評価に関しては、許容される効率性の範囲で実現できる「最適な精度」というものを考えることができる。以下、「最適精度を持つ演算子」とは非特許文献２に示される条件（後述の（９）式）を満たす演算子を意味する。
【非特許文献１】Takeuchi,N., and R.J. Geller, 2000, Physics of the Earth and Planetary Interiors, 119,99-131.
【非特許文献２】Geller,R.J., and N. Takeuchi, 1995, Geophysical Journal International, 123, 449-470.
【非特許文献３】Geller,R.J., and N. Takeuchi, 1998, Geophysical Journal International, 135, 48-62.
【発明の開示】
【発明が解決しようとする課題】
【０００５】
既に述べたように、時間領域有限差分法により行う弾性体からなる有限な媒質中の波動についての数値解析において、従来の計算手法は、安定性、計算効率、精度を同時に好適に実現することができないという問題がある。
【０００６】
本発明は上記問題点を解決するためになされたものであり、時間領域有限差分法により行う弾性体からなる有限な媒質中の波動についての数値解析において、安定性、計算効率、精度を同時に好適に実現することができる演算子を生成する方法、演算子生成装置、及びシミュレーション装置を提供することを目的とする。
【課題を解決するための手段】
【０００７】
本発明に係る演算子生成方法及び演算子生成装置は、予測子修正子法を用いた時間領域有限差分法により行う、弾性体からなる有限な媒質中の波動についての数値解析に関し、運動方程式の陰形式の差分近似である次式、
【数２２】
（ここで、Δｔは時間方向の差分ステップを表し、ｃ＋，ｃ０，ｃ−は解析対象空間に設定した複数の離散点での変位ベクトルの各成分からなる状態ベクトルであって、それぞれ時刻ｔ＋Δｔ，ｔ，ｔ−Δｔにおけるものを表し、ｆは前記各離散点での外力ベクトルを表し、Ｔ′は質量行列を表し、Ｈ′は剛性行列を表す。）に基づいて設定される予測子及び修正子にて用いる演算子を生成するものであって、当該演算子として前記Ｈ′、及び次式、
δＴｎｕｍ＝Ｔ′−Ｔ
で表されるδＴｎｕｍ（ここで、Ｔは前記離散点での差分近似により導かれる正定値の対角質量行列である。）を生成する方法又は手段を含み、前記Ｈ′を、前記状態ベクトルに作用し前記各離散点での前記変位ベクトルの空間１階差分演算子Ｐｉ（ここで、ｉは、Ｉを前記解析対象空間の次元に応じて決まる所定の自然数として、１≦ｉ≦Ｉなる整数である。）と、弾性定数行列Ｃとから次式、
【数２３】
（ここで、ａｉは正のスカラー係数であり、演算子ＰｉＴは演算子Ｐｉを表す行列の転置行列で表される演算子である。）に基づいて生成し、かつ、回転及び並進に対応する前記状態ベクトルに前記演算子δＴｎｕｍを作用させた結果が０となるように当該δＴｎｕｍを設定する。
【０００８】
本発明に係るシミュレーション装置は、予測子修正子法を用いた時間領域有限差分法により行う、弾性体からなる有限な媒質中の波動についての数値解析を行うものであって、運動方程式の陰形式の差分近似である次式、
【数２４】
（ここで、Δｔは時間方向の差分ステップを表し、ｃ＋，ｃ０，ｃ−は解析対象空間に設定した複数の離散点での変位ベクトルの各成分からなる状態ベクトルであって、それぞれ時刻ｔ＋Δｔ，ｔ，ｔ−Δｔにおけるものを表し、ｆは前記各離散点での外力ベクトルを表し、Ｔ′は質量行列を表し、Ｈ′は剛性行列を表す。）に基づいて設定される予測子及び修正子であって、前記離散点での差分近似により導かれる正定値の対角質量行列Ｔ及び、δＴｎｕｍ＝Ｔ′−Ｔで表される演算子δＴｎｕｍを用いて、次式、
【数２５】
で与えられる前記予測子と、次式、
【数２６】
で与えられる前記修正子とをそれぞれ計算して、次式、
【数２７】
で与えられる前記状態ベクトルｃ＋を時間発展的に順次生成する演算手段を有し、前記Ｈ′は、前記状態ベクトルに作用し前記各離散点での前記変位ベクトルの空間１階差分演算子Ｐｉ（ここで、ｉは、Ｉを前記解析対象空間の次元に応じて決まる所定の自然数として、１≦ｉ≦Ｉなる整数である。）と、弾性定数行列Ｃとから次式、
【数２８】
（ここで、ａｉは正のスカラー係数であり、演算子ＰｉＴは演算子Ｐｉを表す行列の転置行列で表される演算子である。）で表される演算子であり、前記演算子δＴｎｕｍは、回転及び並進に対応する前記状態ベクトルに当該演算子δＴｎｕｍを作用させた結果が０となる演算子である。
【発明の効果】
【０００９】
本発明によれば、時間領域有限差分法により行う弾性体からなる有限な媒質中の波動についての数値解析に関し、安定性、計算効率、精度を同時に好適に実現することができる演算子及びシミュレーション装置が得られる。
【発明を実施するための最良の形態】
【００１０】
以下、本発明の実施の形態（以下実施形態という）について説明する。
【００１１】
本発明の実施形態である計算機システムは、演算処理部、記憶部、入力部及び出力部を含んで構成され、例えば、地震波解析に用いられる。演算処理部は、ＣＰＵ（Central Processing Unit：中央演算処理部）として、各種のプログラムを実行する。具体的には、演算処理部は、後述する演算子生成方法に基づく演算子生成処理を行ったり、当該処理で生成される演算子を用いたシミュレーションを行う。これにより、計算機システムは、本発明に係る演算子生成装置、シミュレーション装置を実現する。また、演算処理部は、計算機システムの各部を制御する。
【００１２】
記憶部は、演算処理部で実行される各種プログラムを格納する。また記憶部は、プログラム実行のための作業領域として使用され、計算過程でのデータを一時的に保持するためにも用いられる。また、記憶部は、演算子生成処理の結果として得られる演算子を保持したり、時間発展的なシミュレーションで得られる各時間ステップの解析結果を保持する。
【００１３】
入力部は、キーボードやマウス等を含み、ユーザは入力部を操作して、演算処理部において演算子生成処理やシミュレーションに用いられる各種パラメータ等の設定を行う。
【００１４】
出力部は、プリンタやディスプレイ等を含み、例えば、生成された演算子をプリントアウトしたり、シミュレーション結果をグラフィック表示したりすることができる。
【００１５】
本実施形態に係る演算子生成方法、演算子生成装置及びシミュレーション装置は、時間領域有限差分法により行う、弾性体からなる有限な媒質中の波動についての数値解析に関する。
【００１６】
ここで、演算子生成方法、演算子生成装置及びシミュレーション装置は相互に関連している。すなわち、演算子生成方法は、シミュレーション装置にて実行される時間領域有限差分法による数値解析処理に対応した演算子を生成する方法であり、演算子生成装置は当該演算子生成方法を実行して演算子を生成する。そして、シミュレーション装置は生成された演算子を用いて数値解析を行う。そのため、本発明に係る演算子生成及びシミュレーションの内容を理解する上では、それら全体を１つの計算スキームとして説明することが好適である。そこで、以下、本発明に係る演算子生成及びシミュレーションの内容を含んだ弾性波動解析スキームについて説明する。
【００１７】
［運動方程式］
本発明に関するシミュレーションの対象となる弾性体中の弾性波動伝播に関する運動方程式は以下のような式で表される。
【数２９】
【００１８】
ここで、ρは密度、Ｃｉｊｋｌは弾性定数、ｕｉは変位ベクトルであり、各添え字ｉ，ｊ，ｋ，ｌはそれぞれｘ，ｙ又はｚであり、三次元空間におけるｘｙｚ直交座標系の各軸方向の成分を表す。また、ｕｉ，ｊ はｕｉ のｊ方向の微分（ｊ＝ｘ，ｙ，ｚ）を表す。ｆｉは外力ベクトルのｉ方向成分である。なお、この式は、アインシュタインの縮約規則を用いて表現している。ちなみに当該規則は以下の数式表現においても適宜採用する。
【００１９】
［１ステップの陰形式の差分近似］
時間領域差分法による解析では運動方程式を最適計算スキームに離散化する。離散化により、上記運動方程式の陰形式の差分近似である次式が得られる。
【数３０】
【００２０】
ここで、Δｔは時間ステップを表す。ｃ＋，ｃ０，ｃ−は解析対象空間に設定した複数の離散点に対する変位ベクトルの各成分からなる状態ベクトルであって、それぞれ、時刻ｔ＋Δｔ，ｔ，ｔ−Δｔにおけるものを表す。また、ｆは前記各離散点での外力ベクトル、Ｔ′は質量行列、Ｈ′は剛性行列を表す。
【００２１】
この離散化手法が時間発展に対して完全に安定であるためには、Ｔ′が正定値、Ｈ′が半正定値であることが必要条件となる。ここで、Ｈ′は以下の形式で表現できれば、半正定値性の条件を満たす。
【数３１】
【００２２】
ここで、Ｐは運動方程式中に現れる空間１階微分の差分近似であり、上付き文字Ｔ は行列の転置を表す。また、Ｃは弾性定数行列（空間の不均質性及び異方性も含む）である。また、（３）式右辺に現れる半正定値演算子の線形結合である次式で表されるＨ′も半正定値性の条件を満たす。
【数３２】
【００２３】
ここで、ｉ = １，…，Ｉ（Ｉは整数）で、ａｉは正の係数で、Ｐｉはそれぞれ１階差分演算子である。
【００２４】
［予測子修正子法］
上記陰形式の方程式はそのままの形では計算機システムを用いた処理に適さず、実際には予測子修正子法を用いて計算機システムによるシミュレーションが行われる。予測子修正子法では、（２）式に基づいて、それぞれ陽形式で表現される予測子及び修正子が設定される。予測子及び修正子としてそれぞれ下記（５），（６）式を定義する。
【数３３】
【数３４】
【００２５】
ここで、Ｔは対角質量行列である。また、δＴｎｕｍは、
【数３５】
となる質量行列の残差である。
【００２６】
（５）式左辺括弧内の第１項は、時刻ｔ＋Δｔにおける状態ベクトルｃ＋の予測値である。シミュレーションでは、まず、予測子を用いて既知の状態ベクトルｃ０，ｃ−からｃ＋の予測値を計算する。次いで、修正子を用いて、予測値に対する修正値δｃ＋を計算する。そして、次式に示すように予測子と修正子とから得られた結果を加算して、次のタイムステップの状態ベクトルｃ＋が得られる。
【数３６】
【００２７】
［安定性の実現］
この予測子修正子法の安定性を保証するために、δＴｎｕｍ＝Ｔ′−Ｔを、以下の条件を満たすように設計する。すなわち、修正子で用いられる演算子δＴｎｕｍは、それをもとの運動方程式の固有周波数が０となる固有ベクトル（並進および回転に対応する）に作用させたとき、その結果が０となるという条件と、Ｔ′が正定値となるという条件とを満たすように設計する。
【００２８】
本発明に係る演算子生成方法、演算子生成装置は、このδＴｎｕｍに関する条件と、上述の（３）式又は（４）式の条件とを満たすように、（５），（６）式で用いる演算子Ｈ′，δＴｎｕｍを生成する。そして、予測子修正子法によるシミュレーションは、このように生成された演算子Ｈ′，δＴｎｕｍを用いて行われ、これにより安定性を確保することができる。
【００２９】
［最適精度の実現］
以上、シミュレーションにおいて安定性を実現する演算子の生成について述べてきたが、次に安定性に加え、さらに上述の「最適な精度」が得られる演算子の生成について説明する。
【００３０】
厳密な演算子の表現となる質量行列及び剛性行列をそれぞれＴ０，Ｈ０と定義し、演算子の残差を以下のように定義する。
δＴ＝Ｔ′−Ｔ０
δＨ＝Ｈ′−Ｈ０
【００３１】
また、正確な運動方程式の固有ベクトルをｕｍ（ここでｍは複数存在し得る固有モードの識別子である。）、当該ｕｍに対応する固有周波数をωｍとして、固有ベクトルに対するδＴ及びδＨの行列要素の対角成分をδＴｍｍ及びδＨｍｍと定義する。
【００３２】
Ｔ′及びＨ′が非特許文献２の最適精度演算子となる条件は、次式が空間についての差分ステップの２次以上の精度で成立することである。
【数３７】
【００３３】
この条件を満たす演算子を得るために、離散化した運動方程式に現れる１階差分演算子を以下のように設計する。Ｄｘ，Ｄｙ，Ｄｚを所定の微分定義位置でのそれぞれｘ，ｙ，ｚ軸方向に対する２点１階差分演算子とする。また、Ｅｘｒ，Ｅｙｓ，Ｅｚｔを前記微分定義位置でのｘ軸方向、ｙ軸方向及びｚ軸方向それぞれに関する３点差分近似による、最適精度を持つように設計した２次精度０階差分演算子とする。ここでｒ，ｓ，ｔはそれぞれ１又は２で、３点差分近似に用いる３つの離散点の微分定義位置に対する配置の相違を示すインデックスである。例えば、ｒ＝１は３つの離散点がｘ軸の負方向にシフトした配置、ｒ＝２は正方向にシフトした配置を表す。ここで、各離散点のｘ，ｙ，ｚ各方向における中点の座標をそれぞれｘｃ，ｙｃ，ｚｃとし、離散点の間隔をｈとする。例えば、微分定義位置の座標を（ｘｃ，ｙｃ，ｚｃ）とした場合、ｒ＝１は具体的には、３つの離散点のうち中央の点のｘ座標がｘｃ−ｈ/２となる配置を示し、ｒ＝２は中央の点のｘ座標がｘｃ＋ｈ/２となる配置を示す。他のインデックスｓ，ｔについても同様である。
【００３４】
これらＤｘ，Ｄｙ，Ｄｚ、及びＥｘｒ，Ｅｙｓ，Ｅｚｔを用いて、１階差分演算子Ｄｘｒｓｔ，Ｄｙｒｓｔ，Ｄｚｒｓｔを次式で定義する。
【数３８】
【００３５】
ちなみに、例えば、差分演算子Ｄｘｒｓｔは、Ｄｘを微分方向であるｘ軸とは別の他の２つの軸、すなわちｙ軸及びｚ軸両方向、に拡散した１階微分の差分近似である。Ｄｙｒｓｔ，Ｄｚｒｓｔについても同様である。
【００３６】
例えば、上式中に現れるＥｘ１，Ｅｘ２の具体形を、格子点上の任意連続関数ｆに作用する形で示したものが次に示す（１３），（１４）式それぞれの第１行、第２行の式である。Ｅｙｓ，Ｅｚｔも同様の具体形を有する。
【数３９】
【００３７】
なお、（１３），（１４）式それぞれの第３行の式は第１行、第２行で表される式をテイラー展開した結果を示している。このようにＥｘｒ，Ｅｙｓ，Ｅｚｔは、０階微分（当該テイラー展開の第１項）に特定の係数を有する２階微分成分（当該テイラー展開の第２項）が付加されたものとなり、これにより（９）式で示した最適精度の条件が満たされるようになる。
【００３８】
最適精度条件が満たされる場合には、（４）式の差分演算子Ｐｉは、Ｄｘｒｓｔ，Ｄｙｒｓｔ，Ｄｚｒｓｔを用いてｒ，ｓ，ｔの８通り（２３通り）の組み合わせのそれぞれに対応して定められる。すなわち、（４）式の説明で述べたＩは８であり、ｉは１≦ｉ≦Ｉなる整数である。また、各ａｉは１／８となる。
【００３９】
具体的には、Ｐｒｓｔを、任意なベクトルｖ，ｗを用いた以下の双二次形式で定義する。
【数４０】
【００４０】
ここでＰｒｓｔはＤｉｒｓｔ（ｉ＝ｘ，ｙ，ｚ）の組み合わせを行列表記した１階差分演算子であり、ｒ，ｓ，ｔの組み合わせに応じて番号付けを行い、各Ｐｉと対応させる。
【００４１】
２次元(Ｐ−ＳＶ問題およびＳＨ問題) の場合（解析対象空間をｘｚ平面とする）は、３次元の場合のＤｘｒｓｔ，Ｄｙｒｓｔ，Ｄｚｒｓｔに相当する差分演算子Ｄｘｒｓ，Ｄｚｒｓが次式で定義される。ここで、ｒとｓは上述の３次元の場合のｒとｔに対応するインデックスであり、それぞれの値は１又は２である。
【数４１】
【００４２】
ちなみに、例えば、差分演算子Ｄｘｒｓは、Ｄｘを微分方向であるｘ軸とは別の他の１つの軸、すなわちｚ軸方向に拡散した１階微分の差分近似である。Ｄｚｒｓについても同様である。
【００４３】
最適精度条件が満たされる場合には、（４）式の演算子Ｐｉは、Ｄｘｒｓ，Ｄｚｒｓを用いてｒ，ｓの４通り（２２通り）の組み合わせのそれぞれに対応して定められる。すなわち、２次元の場合、１≦ｉ≦Ｉなる整数ｉの上限値Ｉは４であり、また、各ａｉは１／４となる。
【００４４】
具体的には、２次元の場合のＰｒｓを以下の表式を用いて定義する。
【数４２】
【００４５】
ここでＰｒｓはＤｉｒｓ（ｉ＝ｘ，ｚ）の組み合わせを行列表記した１階差分演算子であり、ｒ，ｓの組で番号付けを行うことで各Ｐｉに対応させる。
【００４６】
この演算子が作用する状態ベクトルとして、ｙ成分の変位のみを持つようにとることでＳＨ問題は定義され、ｘ，ｚ成分のみを持つようにとることでＰ−ＳＶ問題は定義される。
【００４７】
上記のようにＨ′を複数のＰｉの和で表現することにより、数値演算子の一般固有値問題を解いたとき、並進・回転モード以外の固有周波数が０となる固有モードの存在を防ぐことができる。
【００４８】
［非物理現象の抑制］
上記のように定義されたＨ′を用いてシミュレーションを行うと、３次元、２次元のいずれの解析においても、空間一格子ごとに振動するような成分を有する高波数の波が低い振動数を持つようになる現象が生じ得る。物理的にいえば、高波数（短波長）の波はより高い振動数を持つのが自然であり、上述の現象は非物理的な現象である。この現象の発生を抑制する演算子について以下説明する。
【００４９】
非物理現象の発生を抑制するためには、Ｆｉｒで定義されるｉ軸（ｉ＝ｘ，ｙ，ｚ）方向の３点２階差分演算子を用いて演算子Ｑを設計し、上述のＨ′を以下のように修正する。ただし、ｒは、上述のＥｘｒの定義と同様に、差分演算子を定義する離散点の配置のシフト方向を示すインデックスであり、その値は１又は２である。
【数４３】
【００５０】
ここで、ｊ＝１，…，Ｊ（Ｊは問題の次元により定義される所定の整数）であり、εｊは正の係数である。また、Ｑｊは各ｊに対応する２階差分演算子である。（１９）式は２次元、３次元問題共通である。具体的なＦｉｒ，Ｑｊ，εｊの表式は後に記載する。
【００５１】
ここまでで、安定でかつ最適精度の条件を満たし、非物理現象の励起を抑制する剛性行列Ｈ′の定義が完了した。
【００５２】
［質量行列の定義］
次に、安定でかつ最適精度条件に適合し得る正定値質量行列Ｔ′を示す。最適精度条件に適合し得る質量行列Ｔ′を、予測子計算に用いる対角質量行列Ｔと、残差δＴｎｕｍとを用いて（７）式と同じ形式の次式で表す。
【数４４】
【００５３】
前述のように、予測子修正子法の安定性のために、δＴｎｕｍは、離散化された剛性行列Ｈ′と同様に並進及び回転モードに作用させたとき、その値が０となる必要がある。この条件を満たしつつ、最適精度の条件も満たすために、Ｈ′の表式を利用して、以下のようにδＴｎｕｍを定義する。
【数４５】
【００５４】
この式の右辺は剛性行列の表現における弾性行列ＣをＸで置き換えたものである。Ｘは、最適精度条件を満たすように、密度および格子間隔Δｘａ（ａはｘ，ｙ，ｚのいずれか）を用いて定義される値で、
【数４６】
となるものである。ここで、ρ(α) はα番目の微分定義点における密度の値であり、δａｂ，δｃｄ，δａｃはクロネッカーのデルタであり、ａ，ｂ，ｃ，ｄはそれぞれ座標軸ｘ，ｙ，ｚに対応する添え字である。
【００５５】
ここまでで、予測子修正子法を用いた手法を安定にかつ、最適な精度で計算するための、剛性行列Ｈ′と、質量行列Ｔ′（及びδＴｎｕｍ）を導出することができた。
【００５６】
これらの演算子を用い、予測子修正子法を用いて時間発展計算を行うことで、時刻格子、空間格子上における任意の点における状態ベクトルで表せる変位場を計算することが可能となる。
【００５７】
［シミュレーションにおける計算効率の向上］
予測子修正子法におけるほとんどの計算負荷は、Ｔ′，Ｈ′による空間差分演算部分になる。この空間差分演算を以下の手順で行うことにより、シミュレーション実行時の計算効率を向上させることができる。
【００５８】
上述したように、Ｈ′，Ｔ′は（３）式のような形式の行列の積を含んでいる。例えば、（３）式で表される演算子を状態ベクトルに作用させる演算では、
ｓｔｅｐ１：差分演算子行列Ｐと格子点における変位の値に対応する状態ベクトルｕとの積、
ｓｔｅｐ２：step１の差分演算結果の各要素を弾性定数倍する演算、
ｓｔｅｐ３：step２の演算結果を、差分演算子の転置であるＰＴ倍する演算、
を連続的に、効率よく行う必要がある。なお、弾性定数は微分定義位置において定義されているものとする。
【００５９】
３次元上の格子点上の変位を、通常の有限要素法で用いられるように番号付けを行い、１次元ベクトルに置き換えたものをｕｊとし、Ｐの要素をＰｉｊとすると、上記step１の演算は、
【数４７】
と書き表される。ｕｉは格子点での値で定義されるベクトルであるのに対し、ｖｉは微分定義位置における演算結果を１次元的に番号付けして整列したベクトルである。（２３）式は、或る微分定義位置(i) におけるｖｉを計算するために微分定義位置ｉ周辺の各ｕｊの値にＰｉｊで定義される重みを乗算し、和をとる形になっている。図１は、２次元の場合の（２３）式の演算を説明する図であり、ステンシルの模式図である。図１において、白丸（○）が格子点であり、黒四角（■）が微分定義位置を表す。各格子点の近くに示す数値はＰｉｊを示す。なお、一般的には微分定義位置に対して図１（ａ）のように６点の格子点が存在し得るが、微分定義位置が解析対象空間の境界に隣接する場合、図１（ｂ）のように４点からなるステンシルが設定され得る。
【００６０】
ｓｔｅｐ２の演算は、単純なスカラー倍で書き表されるので特に工夫は必要でない。
【００６１】
ｓｔｅｐ３の演算は、
【数４８】
と書き表される。（２４）式において、ｗｋは格子点における演算結果を表すベクトルである。また（２４）式において、ｗｋ（ｊ）＝Ｐｊｋｖｊとおいた。ここで、ｗｋ（ｊ）はｖｊで表されるベクトルのｊ番目の微分定義点の値を、ｋ番目の格子点に分配することで計算される。すなわち、右辺のｖｉが上述のように微分定義位置での値で定義されるベクトルであるのに対し、ｗｋはｕｉと同じく格子点における演算結果を表すベクトルである。図２は、２次元の場合の（２４）式の演算を説明する図であり、ステンシルの模式図である。図２は図１と同様、白丸（○）が格子点であり、黒四角（■）が微分定義位置を表す。また、各格子点の近くに示す数値はＰｉｊを示す。なお、図２（ａ），（ｂ）でのステンシルの違いも図１と同様である。
【００６２】
すべての微分定義位置について、上記step１〜３の演算を繰り返し行うことで、全ての行列演算が完了する。この計算方法によれば巨大な行列とベクトルの積の計算を部分的な演算の和で効率よく実行することができる。
【００６３】
［Ｑの具体的な表現］
以下に上述のＦｉｒ，Ｑｊ及びεｊの具体形について表記する。まず３次元の場合について記述し、その後で２次元の場合について記述する。
【００６４】
はじめに、２階差分演算子Ｆｉｒを０階差分演算子を定義したときと同様に、連続関数ｆに対する離散化表現を用いて定義する。以下にＦｘ１の場合の具体形を示す。
【数４９】
【００６５】
他のＦｉｒについても同様に定義される。これらのＦｉｒの定義を用い、（１５）式により表されるＰｒｓｔの定義式におけるＤｉｒｓｔを以下の６パターンで置き換え、Ｐｒｓｔに相当する部分を抽出することによってＱｊｒｓｔ を定義する。
【数５０】
【００６６】
これらの６パターンの置き換えにより生成される２階差分演算子Ｑｊｒｓｔを用いて、以下のように剛性行列の要素Ｈｒｓｔを定義する。
【数５１】
【００６７】
そして、最終的に剛性行列Ｈ′を
【数５２】
で計算する。上式ではｊ＝１，２，３に対してεｊ＝５/１４４，ｊ＝４，５，６に対してεｊ＝ (５/１４４)２を用いた。５/１４４以上の値を使用した手法も可能であるが、その場合時間ステップΔt をより小さくとる必要があるため、５/１４４を用いるのが好適である。
【００６８】
２次元の場合の２階差分演算子Ｑｊｒｓの表現は３次元の場合と同様に（１８）式により表されるＰｒｓの定義式におけるＤｉｒｓを以下のように置き換えたものに対して、Ｐｒｓに相当する演算子を抽出することにより計算することができる。
【数５３】
【００６９】
これらを用いて構成される２階差分演算子Ｑｊｒｓを用いて、２次元の場合の剛性行列の要素Ｈｒｓを定義する。
【数５４】
【００７０】
上式ではεｊ＝５/１４４，（ｊ＝１，２）を用いるが、２次元の場合は４/１４４ 以上であれば、非物理的な現象を抑制することができる。これらの剛性行列の要素Ｈｒｓを用いて最終的に剛性行列Ｈ′は
【数５５】
と書き表される。
【００７１】
［上述の計算スキームの応用例等］
本計算スキームでは、各微分定義点に弾性定数(異方性を含む)および密度の値を与えることで、全ての行列を定義する。任意不均質構造は、弾性定数、密度の値を各微分定義点に与えることで実現でき、また、任意の表面地形形状は、空中の弾性定数および密度に０を与えることで導入できる。
【００７２】
また、本計算スキームは陽解法を用いて時間を更新していく手法であり、時間ステップΔｔのとりうる最大値は、固有値解析により求められる最大固有値により決定することが可能である。これは一般的な手法であるので、ここでは特に記述しない。
【００７３】
上述の演算子生成方法により定義される安定かつ最適精度を持つ空間差分演算子は、既存の吸収境界条件を用いる手法と組み合わせて使用することもできる。よって、このような手法も本発明の演算子生成方法・装置及びシミュレーション装置の範囲に含まれる。
【００７４】
また、上記により定義される安定かつ最適精度を持つ空間差分演算子は、既存の非弾性減衰を表現する手法と組み合わせて使用することができる。このような手法も本発明の範囲に含まれる。
【００７５】
本発明においては３次元および２次元問題についての手法についてのみ扱う。これは、１次元の場合は、非特許文献３によって定義される手法と等価であるためである。
【００７６】
「時間領域差分法」には、(a)「規則的な格子上で離散化し、変位もしくは速度、加速度を変数とし、時間発展を計算する手法」と、(b)「半格子幅での食い違いをもつ２種類の格子上で離散化し、変位もしくは速度、加速度、および応力を変数とし、時間発展を計算する手法」とが含まれる。ここまでは(a) について説明してきたが、本手法で用いた空間差分演算子を(b) の手法に適用することが可能である。すなわち、変位および速度の微分にＰ（もしくはＰＴ ）を用い、応力の微分にはＰＴ（もしくはＰ）を用いて、時間発展計算を行うことで、安定性が保障されることになる。このような手法の場合、非特許文献３で示されているように、予測子修正子法は適用できないため、高精度性は保障できないが、予測子のみを用いて時間発展計算を行うことで、安定化の実現が可能である。
【００７７】
なお、一般に解析対象空間に設定される格子点の数は非常に大きく、例えば、行列形式で表される離散化した演算子の規模も巨大となる。そのため、当該巨大な演算子行列の全体を生成して記憶部に保持したり、記憶部に保持した巨大な行列と、巨大な状態ベクトルとを直接、乗算することは、計算効率の観点から得策ではない。この点については、既に１つの工夫を、図１及び図２を用いて説明した。このように、本発明に係るシミュレーション装置では、状態ベクトルに演算子行列を作用させる計算を、解析対象空間の一部分ずつ行い得る。そして、その場合には、演算子行列のうち、これから計算を行う一部分に対応した行列要素が記憶部上に存在すれば処理を実行することができる。本発明に係る演算子生成方法、演算子生成装置は、このように必要とする一部分ずつを逐次的にシミュレーション装置に提供する演算子生成の仕方を含んでいる。
【００７８】
また、演算子Ｈ′，δＴｎｕｍは上述のように演算子（行列）の積を含む形式で表され得る。この場合に、演算子を生成するとは、その行列の積を計算した結果得られる行列Ｈ′や行列δＴｎｕｍの行列要素を算出することには限定されず、積を構成する個々の演算子行列の要素を求めることも含む。例えば、図１及び図２を用いて説明した処理方法では、本発明の演算子生成では、Ｈ′を構成する行列の積ＰｉＴＣＰｉの１階差分演算子ＰｉＴ，Ｐｉの行列要素は直接的・具体的に定めるが、Ｈ′の行列要素の具体的な値までは計算しない。しかし、この場合においてもＰｉＴ，Ｐｉに基づいてＨ′は（４）式により定義され、この意味においてＨ′は生成されているということができる。
【産業上の利用可能性】
【００７９】
本発明に係る演算子生成方法、演算子生成装置、及びシミュレーション装置は、地震波の解析に用いることができる。地球内部を伝播してきた地震波形記録が持つ情報を十分に活用することができれば、より詳細で正確な内部構造推定が可能となる。高精度かつ高解像度の地球内部構造モデルの推定は、例えば、防災や、油田開発における物理探査などに利用することができる。
【図面の簡単な説明】
【００８０】
【図１】ＰＴＣＰの形式の行列の積と状態ベクトルｕとの積の演算方法を２次元の場合について説明するステンシルの模式図である。
【図２】ＰＴＣＰの形式の行列の積と状態ベクトルｕとの積の演算方法を２次元の場合について説明するステンシルの模式図である。
【特許請求の範囲】
【請求項１】
予測子修正子法を用いた時間領域有限差分法により行う、弾性体からなる有限な媒質中の波動についての数値解析に関し、運動方程式の陰形式の差分近似である次式、
【数１】
（ここで、Δｔは時間方向の差分ステップを表し、ｃ＋，ｃ０，ｃ−は解析対象空間に設定した複数の離散点での変位ベクトルの各成分からなる状態ベクトルであって、それぞれ時刻ｔ＋Δｔ，ｔ，ｔ−Δｔにおけるものを表し、ｆは前記各離散点での外力ベクトルを表し、Ｔ′は質量行列を表し、Ｈ′は剛性行列を表す。）に基づいて設定される予測子及び修正子にて用いる演算子を生成する演算子生成方法であって、
当該演算子として前記Ｈ′、及び次式、
δＴｎｕｍ＝Ｔ′−Ｔ
で表されるδＴｎｕｍ（ここで、Ｔは前記離散点での差分近似により導かれる正定値の対角質量行列である。）を生成する方法を含み、
前記Ｈ′を、前記状態ベクトルに作用し前記各離散点での前記変位ベクトルの空間１階差分演算子Ｐｉ（ここで、ｉは、Ｉを前記解析対象空間の次元に応じて決まる所定の自然数として、１≦ｉ≦Ｉなる整数である。）と、弾性定数行列Ｃとから次式、
【数２】
（ここで、ａｉは正のスカラー係数であり、演算子ＰｉＴは演算子Ｐｉを表す行列の転置行列で表される演算子である。）に基づいて生成し、
回転及び並進に対応する前記状態ベクトルに前記演算子δＴｎｕｍを作用させた結果が０となるように当該δＴｎｕｍを設定すること、
を特徴とする演算子生成方法。
【請求項２】
請求項１に記載の演算子生成方法において、
前記解析対象空間はｘ軸、ｙ軸及びｚ軸を座標軸とする３次元空間であって、
前記各離散点に対応して設定される微分定義位置でのｘ軸方向、ｙ軸方向及びｚ軸方向それぞれに関する２点１階差分演算子をＤｘ，Ｄｙ，Ｄｚとし、
前記微分定義位置でのｘ軸方向、ｙ軸方向及びｚ軸方向それぞれに関する３点差分近似による、最適精度を持つように設計した２次精度０階差分演算子をＥｘｒ，Ｅｙｓ，Ｅｚｔ（ここでｒ，ｓ，ｔは１又は２で、差分に用いる３つの前記離散点の前記微分定義位置に対する配置の相違に応じた当該演算子の種類を示す。）とし、
次式、
【数３】
により、前記Ｄｘ，Ｄｙ，Ｄｚを微分方向の軸とは別の他の２つの軸方向に拡散した１階差分演算子Ｄｘｒｓｔ，Ｄｙｒｓｔ，Ｄｚｒｓｔを前記各微分定義位置にて定義した場合に、
前記各演算子Ｐｉ（ここで、前記Ｉは８であり、ｉは１≦ｉ≦Ｉなる整数である。）は、前記Ｄｘｒｓｔ，Ｄｙｒｓｔ，Ｄｚｒｓｔを用いて前記ｒ，ｓ，ｔの８通りの組み合わせのそれぞれに対応して定められ、前記各ａｉは１／８であること、
を特徴とする演算子生成方法。
【請求項３】
請求項１に記載の演算子生成方法において、
前記解析対象空間はｘ軸及びｚ軸を座標軸とする２次元空間であって、
前記各離散点に対応して設定される微分定義位置でのｘ軸方向及びｚ軸方向それぞれに関する２点１階差分演算子をＤｘ，Ｄｚとし、
前記微分定義位置でのｘ軸方向及びｚ軸方向それぞれに関する３点差分近似による、最適精度を持つように設計した２次精度０階差分演算子をＥｘｒ，Ｅｚｓ（ここでｒ，ｓは１又は２で、差分に用いる３つの前記離散点の前記微分定義位置に対する配置の相違に応じた当該演算子の種類を示す。）とし、
次式、
【数４】
により、前記Ｄｘ，Ｄｚを微分方向の軸とは別の軸方向に拡散した１階差分演算子Ｄｘｒｓ，Ｄｚｒｓを前記各微分定義位置にて定義した場合に、
前記各演算子Ｐｉ（ここで、前記Ｉは４であり、ｉは１≦ｉ≦Ｉなる整数である。）は、前記Ｄｘｒｓ，Ｄｚｒｓを用いて前記ｒ，ｓの４通りの組み合わせのそれぞれに対応して定められ、前記各ａｉは１／４であること、
を特徴とする演算子生成方法。
【請求項４】
請求項２又は請求項３に記載の演算子生成方法において、
前記Ｈ′を、前記状態ベクトルに作用し前記各離散点での前記変位ベクトルの空間２階差分演算子Ｑｊ（ここで、ｊは、Ｊを前記解析対象空間の次元に応じて決まる所定の自然数として、１≦ｊ≦Ｊなる整数である。）を用いて修正した次式、
【数５】
（ここで、εｊは正のスカラー係数であり、演算子ＱｊＴは演算子Ｑｊを表す行列の転置行列で表される演算子である。）に基づいて生成し、
前記εｊは、当該Ｈ′の前記微分定義位置におけるテイラー展開の４次の項の係数が０以上となるように定められること、
を特徴とする演算子生成方法。
【請求項５】
請求項４に記載の演算子生成方法において、
前記質量行列を表す正確な演算子Ｔ０に対する前記Ｔ′の誤差をδＴとし、前記剛性行列を表す正確な演算子Ｈ０に対する前記Ｈ′の誤差をδＨとし、正確な運動方程式の固有ベクトルｕｍ（ここでｍは複数存在し得る固有モードの識別子である。）、当該ｕｍに対応する固有周波数をωｍとして、δＴｍｍ及びδＨｍｍを固有ベクトルに対するδＴ及びδＨの行列要素の対角成分とした場合に、
前記Ｔ′及び前記Ｈ′が全ての前記固有モードについて関係式、
ωｍ２δＴｍｍ−δＨｍｍ＝０
を空間についての差分間隔の２次以上の精度で近似的に満たすという条件の下で、対角行列Ｘを選択して、次式、
【数６】
で表される前記δＴｎｕｍを生成すること、
を特徴とする演算子生成方法。
【請求項６】
予測子修正子法を用いた時間領域有限差分法により行う、弾性体からなる有限な媒質中の波動についての数値解析に関し、運動方程式の陰形式の差分近似である次式、
【数７】
（ここで、Δｔは時間方向の差分ステップを表し、ｃ＋，ｃ０，ｃ−は解析対象空間に設定した複数の離散点での変位ベクトルの各成分からなる状態ベクトルであって、それぞれ時刻ｔ＋Δｔ，ｔ，ｔ−Δｔにおけるものを表し、ｆは前記各離散点での外力ベクトルを表し、Ｔ′は質量行列を表し、Ｈ′は剛性行列を表す。）に基づいて設定される予測子及び修正子にて用いる演算子を生成する演算子生成装置であって、
当該演算子として前記Ｈ′、及び次式、
δＴｎｕｍ＝Ｔ′−Ｔ
で表されるδＴｎｕｍ（ここで、Ｔは前記離散点での差分近似により導かれる正定値の対角質量行列である。）を生成する手段を有し、
前記Ｈ′を、前記状態ベクトルに作用し前記各離散点での前記変位ベクトルの空間１階差分演算子Ｐｉ（ここで、ｉは、Ｉを前記解析対象空間の次元に応じて決まる所定の自然数として、１≦ｉ≦Ｉなる整数である。）と、弾性定数行列Ｃとから次式、
【数８】
（ここで、ａｉは正のスカラー係数であり、演算子ＰｉＴは演算子Ｐｉを表す行列の転置行列で表される演算子である。）に基づいて生成し、
回転及び並進に対応する前記状態ベクトルに前記演算子δＴｎｕｍを作用させた結果が０となるように当該δＴｎｕｍを設定すること、
を特徴とする演算子生成装置。
【請求項７】
請求項６に記載の演算子生成装置において、
前記解析対象空間はｘ軸、ｙ軸及びｚ軸を座標軸とする３次元空間であって、
前記各離散点に対応して設定される微分定義位置でのｘ軸方向、ｙ軸方向及びｚ軸方向それぞれに関する２点１階差分演算子をＤｘ，Ｄｙ，Ｄｚとし、
前記微分定義位置でのｘ軸方向、ｙ軸方向及びｚ軸方向それぞれに関する３点差分近似による、最適精度を持つように設計した２次精度０階差分演算子をＥｘｒ，Ｅｙｓ，Ｅｚｔ（ここでｒ，ｓ，ｔは１又は２で、差分に用いる３つの前記離散点の前記微分定義位置に対する配置の相違に応じた当該演算子の種類を示す。）とし、
次式、
【数９】
により、前記Ｄｘ，Ｄｙ，Ｄｚを微分方向の軸とは別の他の２つの軸方向に拡散した１階差分演算子Ｄｘｒｓｔ，Ｄｙｒｓｔ，Ｄｚｒｓｔを前記各微分定義位置にて定義した場合に、
前記各演算子Ｐｉ（ここで、前記Ｉは８であり、ｉは１≦ｉ≦Ｉなる整数である。）は、前記Ｄｘｒｓｔ，Ｄｙｒｓｔ，Ｄｚｒｓｔを用いて前記ｒ，ｓ，ｔの８通りの組み合わせのそれぞれに対応して定められ、前記各ａｉは１／８であること、
を特徴とする演算子生成装置。
【請求項８】
請求項６に記載の演算子生成装置において、
前記解析対象空間はｘ軸及びｚ軸を座標軸とする２次元空間であって、
前記各離散点に対応して設定される微分定義位置でのｘ軸方向及びｚ軸方向それぞれに関する２点１階差分演算子をＤｘ，Ｄｚとし、
前記微分定義位置でのｘ軸方向及びｚ軸方向それぞれに関する３点差分近似による、最適精度を持つように設計した２次精度０階差分演算子をＥｘｒ，Ｅｚｓ（ここでｒ，ｓは１又は２で、差分に用いる３つの前記離散点の前記微分定義位置に対する配置の相違に応じた当該演算子の種類を示す。）とし、
次式、
【数１０】
により、前記Ｄｘ，Ｄｚを微分方向の軸とは別の軸方向に拡散した１階差分演算子Ｄｘｒｓ，Ｄｚｒｓを前記各微分定義位置にて定義した場合に、
前記各演算子Ｐｉ（ここで、前記Ｉは４であり、ｉは１≦ｉ≦Ｉなる整数である。）は、前記Ｄｘｒｓ，Ｄｚｒｓを用いて前記ｒ，ｓの４通りの組み合わせのそれぞれに対応して定められ、前記各ａｉは１／４であること、
を特徴とする演算子生成装置。
【請求項９】
請求項７又は請求項８に記載の演算子生成装置において、
前記Ｈ′を、前記状態ベクトルに作用し前記各離散点での前記変位ベクトルの空間２階差分演算子Ｑｊ（ここで、ｊは、Ｊを前記解析対象空間の次元に応じて決まる所定の自然数として、１≦ｊ≦Ｊなる整数である。）を用いて修正した次式、
【数１１】
（ここで、εｊは正のスカラー係数であり、演算子ＱｊＴは演算子Ｑｊを表す行列の転置行列で表される演算子である。）に基づいて生成し、
前記εｊは、当該Ｈ′の前記微分定義位置におけるテイラー展開の４次の項の係数が０以上となるように定められること、
を特徴とする演算子生成装置。
【請求項１０】
請求項９に記載の演算子生成装置において、
前記質量行列を表す正確な演算子Ｔ０に対する前記Ｔ′の誤差をδＴとし、前記剛性行列を表す正確な演算子Ｈ０に対する前記Ｈ′の誤差をδＨとし、正確な運動方程式の固有ベクトルｕｍ（ここでｍは複数存在し得る固有モードの識別子である。）、当該ｕｍに対応する固有周波数をωｍとして、δＴｍｍ及びδＨｍｍを固有ベクトルに対するδＴ及びδＨの行列要素の対角成分とした場合に、
前記Ｔ′及び前記Ｈ′が全ての前記固有モードについて関係式、
ωｍ２δＴｍｍ−δＨｍｍ＝０
を空間についての差分間隔の２次以上の精度で近似的に満たすという条件の下で、対角行列Ｘを選択して、次式、
【数１２】
で表される前記δＴｎｕｍを生成すること、
を特徴とする演算子生成装置。
【請求項１１】
予測子修正子法を用いた時間領域有限差分法により行う、弾性体からなる有限な媒質中の波動についての数値解析を行うシミュレーション装置であって、
運動方程式の陰形式の差分近似である次式、
【数１３】
（ここで、Δｔは時間方向の差分ステップを表し、ｃ＋，ｃ０，ｃ−は解析対象空間に設定した複数の離散点での変位ベクトルの各成分からなる状態ベクトルであって、それぞれ時刻ｔ＋Δｔ，ｔ，ｔ−Δｔにおけるものを表し、ｆは前記各離散点での外力ベクトルを表し、Ｔ′は質量行列を表し、Ｈ′は剛性行列を表す。）に基づいて設定される予測子及び修正子であって、前記離散点での差分近似により導かれる正定値の対角質量行列Ｔ及び、δＴｎｕｍ＝Ｔ′−Ｔで表される演算子δＴｎｕｍを用いて、次式、
【数１４】
で与えられる前記予測子と、次式、
【数１５】
で与えられる前記修正子とをそれぞれ計算して、次式、
【数１６】
で与えられる前記状態ベクトルｃ＋を時間発展的に順次生成する演算手段を有し、
前記Ｈ′は、前記状態ベクトルに作用し前記各離散点での前記変位ベクトルの空間１階差分演算子Ｐｉ（ここで、ｉは、Ｉを前記解析対象空間の次元に応じて決まる所定の自然数として、１≦ｉ≦Ｉなる整数である。）と、弾性定数行列Ｃとから次式、
【数１７】
（ここで、ａｉは正のスカラー係数であり、演算子ＰｉＴは演算子Ｐｉを表す行列の転置行列で表される演算子である。）で表される演算子であり、
前記演算子δＴｎｕｍは、回転及び並進に対応する前記状態ベクトルに当該演算子δＴｎｕｍを作用させた結果が０となる演算子であること、
を特徴とするシミュレーション装置。
【請求項１２】
請求項１１に記載のシミュレーション装置において、
前記解析対象空間はｘ軸、ｙ軸及びｚ軸を座標軸とする３次元空間であって、
前記各離散点に対応して設定される微分定義位置でのｘ軸方向、ｙ軸方向及びｚ軸方向それぞれに関する２点１階差分演算子をＤｘ，Ｄｙ，Ｄｚとし、
前記微分定義位置でのｘ軸方向、ｙ軸方向及びｚ軸方向それぞれに関する３点差分近似による、最適精度を持つように設計した２次精度０階差分演算子をＥｘｒ，Ｅｙｓ，Ｅｚｔ（ここでｒ，ｓ，ｔは１又は２で、差分に用いる３つの前記離散点の前記微分定義位置に対する配置の相違に応じた当該演算子の種類を示す。）とし、
次式、
【数１８】
により、前記Ｄｘ，Ｄｙ，Ｄｚを微分方向の軸とは別の他の２つの軸方向に拡散した１階差分演算子Ｄｘｒｓｔ，Ｄｙｒｓｔ，Ｄｚｒｓｔを前記各微分定義位置にて定義した場合に、
前記Ｈ′は、前記Ｄｘｒｓｔ，Ｄｙｒｓｔ，Ｄｚｒｓｔを用いて前記ｒ，ｓ，ｔの８通りの組み合わせのそれぞれに対応して定められた前記演算子Ｐｉ（ここで、前記Ｉは８であり、ｉは１≦ｉ≦Ｉなる整数である。）を用い、前記各ａｉを１／８として生成されるものであること、
を特徴とするシミュレーション装置。
【請求項１３】
請求項１１に記載のシミュレーション装置において、
前記解析対象空間はｘ軸及びｚ軸を座標軸とする２次元空間であって、
前記各離散点に対応して設定される微分定義位置でのｘ軸方向及びｚ軸方向それぞれに関する２点１階差分演算子をＤｘ，Ｄｚとし、
前記微分定義位置でのｘ軸方向及びｚ軸方向それぞれに関する３点差分近似による、最適精度を持つように設計した２次精度０階差分演算子をＥｘｒ，Ｅｚｓ（ここでｒ，ｓは１又は２で、差分に用いる３つの前記離散点の前記微分定義位置に対する配置の相違に応じた当該演算子の種類を示す。）とし、
次式、
【数１９】
により、前記Ｄｘ，Ｄｚを微分方向の軸とは別の軸方向に拡散した１階差分演算子Ｄｘｒｓ，Ｄｚｒｓを前記各微分定義位置にて定義した場合に、
前記Ｈ′は、前記Ｄｘｒｓ，Ｄｚｒｓを用いて前記ｒ，ｓの４通りの組み合わせのそれぞれに対応して定められた前記演算子Ｐｉ（ここで、前記Ｉは４であり、ｉは１≦ｉ≦Ｉなる整数である。）を用い、前記各ａｉを１／４として生成されるものであること、
を特徴とするシミュレーション装置。
【請求項１４】
請求項１２又は請求項１３に記載のシミュレーション装置において、
前記Ｈ′は、前記状態ベクトルに作用し前記各離散点での前記変位ベクトルの空間２階差分演算子Ｑｊ（ここで、ｊは、Ｊを前記解析対象空間の次元に応じて決まる所定の自然数として、１≦ｊ≦Ｊなる整数である。）を用いて修正した次式、
【数２０】
（ここで、εｊは正のスカラー係数であり、演算子ＱｊＴは演算子Ｑｊを表す行列の転置行列で表される演算子である。）で表される演算子であり、
前記εｊは、当該Ｈ′の前記微分定義位置におけるテイラー展開の４次の項の係数が０以上となるように定められていること、
を特徴とするシミュレーション装置。
【請求項１５】
請求項１４に記載のシミュレーション装置において、
前記質量行列を表す正確な演算子Ｔ０に対する前記Ｔ′の誤差をδＴとし、前記剛性行列を表す正確な演算子Ｈ０に対する前記Ｈ′の誤差をδＨとし、正確な運動方程式の固有ベクトルｕｍ（ここでｍは複数存在し得る固有モードの識別子である。）、当該ｕｍに対応する固有周波数をωｍとして、δＴｍｍ及びδＨｍｍを固有ベクトルに対するδＴ及びδＨの行列要素の対角成分とした場合に、
前記Ｔ′及び前記Ｈ′は、全ての前記固有モードについて関係式、
ωｍ２δＴｍｍ−δＨｍｍ＝０
を空間についての差分間隔の２次以上の精度で近似的に満たし、
前記δＴｎｕｍは、対角行列Ｘを用いて、次式、
【数２１】
で表されるものであること、
を特徴とするシミュレーション装置。
【請求項１６】
速度もしくは変位、及び応力を互いに食い違った格子上で定義するスタッガードグリッド法を用いた時間領域有限差分法により弾性体からなる有限な媒質中の波動についての数値解析を行うシミュレーション装置であって、
前記速度もしくは変位の空間１階差分演算子としてＰを用い、
前記応力の空間１階差分演算子として前記演算子Ｐを表す行列の転置行列で表される演算子ＰＴを用い、
解析対象空間における自由表面を、垂直速度成分を持つ格子上に配置すること、
を特徴とするシミュレーション装置。
【請求項１】
予測子修正子法を用いた時間領域有限差分法により行う、弾性体からなる有限な媒質中の波動についての数値解析に関し、運動方程式の陰形式の差分近似である次式、
【数１】
（ここで、Δｔは時間方向の差分ステップを表し、ｃ＋，ｃ０，ｃ−は解析対象空間に設定した複数の離散点での変位ベクトルの各成分からなる状態ベクトルであって、それぞれ時刻ｔ＋Δｔ，ｔ，ｔ−Δｔにおけるものを表し、ｆは前記各離散点での外力ベクトルを表し、Ｔ′は質量行列を表し、Ｈ′は剛性行列を表す。）に基づいて設定される予測子及び修正子にて用いる演算子を生成する演算子生成方法であって、
当該演算子として前記Ｈ′、及び次式、
δＴｎｕｍ＝Ｔ′−Ｔ
で表されるδＴｎｕｍ（ここで、Ｔは前記離散点での差分近似により導かれる正定値の対角質量行列である。）を生成する方法を含み、
前記Ｈ′を、前記状態ベクトルに作用し前記各離散点での前記変位ベクトルの空間１階差分演算子Ｐｉ（ここで、ｉは、Ｉを前記解析対象空間の次元に応じて決まる所定の自然数として、１≦ｉ≦Ｉなる整数である。）と、弾性定数行列Ｃとから次式、
【数２】
（ここで、ａｉは正のスカラー係数であり、演算子ＰｉＴは演算子Ｐｉを表す行列の転置行列で表される演算子である。）に基づいて生成し、
回転及び並進に対応する前記状態ベクトルに前記演算子δＴｎｕｍを作用させた結果が０となるように当該δＴｎｕｍを設定すること、
を特徴とする演算子生成方法。
【請求項２】
請求項１に記載の演算子生成方法において、
前記解析対象空間はｘ軸、ｙ軸及びｚ軸を座標軸とする３次元空間であって、
前記各離散点に対応して設定される微分定義位置でのｘ軸方向、ｙ軸方向及びｚ軸方向それぞれに関する２点１階差分演算子をＤｘ，Ｄｙ，Ｄｚとし、
前記微分定義位置でのｘ軸方向、ｙ軸方向及びｚ軸方向それぞれに関する３点差分近似による、最適精度を持つように設計した２次精度０階差分演算子をＥｘｒ，Ｅｙｓ，Ｅｚｔ（ここでｒ，ｓ，ｔは１又は２で、差分に用いる３つの前記離散点の前記微分定義位置に対する配置の相違に応じた当該演算子の種類を示す。）とし、
次式、
【数３】
により、前記Ｄｘ，Ｄｙ，Ｄｚを微分方向の軸とは別の他の２つの軸方向に拡散した１階差分演算子Ｄｘｒｓｔ，Ｄｙｒｓｔ，Ｄｚｒｓｔを前記各微分定義位置にて定義した場合に、
前記各演算子Ｐｉ（ここで、前記Ｉは８であり、ｉは１≦ｉ≦Ｉなる整数である。）は、前記Ｄｘｒｓｔ，Ｄｙｒｓｔ，Ｄｚｒｓｔを用いて前記ｒ，ｓ，ｔの８通りの組み合わせのそれぞれに対応して定められ、前記各ａｉは１／８であること、
を特徴とする演算子生成方法。
【請求項３】
請求項１に記載の演算子生成方法において、
前記解析対象空間はｘ軸及びｚ軸を座標軸とする２次元空間であって、
前記各離散点に対応して設定される微分定義位置でのｘ軸方向及びｚ軸方向それぞれに関する２点１階差分演算子をＤｘ，Ｄｚとし、
前記微分定義位置でのｘ軸方向及びｚ軸方向それぞれに関する３点差分近似による、最適精度を持つように設計した２次精度０階差分演算子をＥｘｒ，Ｅｚｓ（ここでｒ，ｓは１又は２で、差分に用いる３つの前記離散点の前記微分定義位置に対する配置の相違に応じた当該演算子の種類を示す。）とし、
次式、
【数４】
により、前記Ｄｘ，Ｄｚを微分方向の軸とは別の軸方向に拡散した１階差分演算子Ｄｘｒｓ，Ｄｚｒｓを前記各微分定義位置にて定義した場合に、
前記各演算子Ｐｉ（ここで、前記Ｉは４であり、ｉは１≦ｉ≦Ｉなる整数である。）は、前記Ｄｘｒｓ，Ｄｚｒｓを用いて前記ｒ，ｓの４通りの組み合わせのそれぞれに対応して定められ、前記各ａｉは１／４であること、
を特徴とする演算子生成方法。
【請求項４】
請求項２又は請求項３に記載の演算子生成方法において、
前記Ｈ′を、前記状態ベクトルに作用し前記各離散点での前記変位ベクトルの空間２階差分演算子Ｑｊ（ここで、ｊは、Ｊを前記解析対象空間の次元に応じて決まる所定の自然数として、１≦ｊ≦Ｊなる整数である。）を用いて修正した次式、
【数５】
（ここで、εｊは正のスカラー係数であり、演算子ＱｊＴは演算子Ｑｊを表す行列の転置行列で表される演算子である。）に基づいて生成し、
前記εｊは、当該Ｈ′の前記微分定義位置におけるテイラー展開の４次の項の係数が０以上となるように定められること、
を特徴とする演算子生成方法。
【請求項５】
請求項４に記載の演算子生成方法において、
前記質量行列を表す正確な演算子Ｔ０に対する前記Ｔ′の誤差をδＴとし、前記剛性行列を表す正確な演算子Ｈ０に対する前記Ｈ′の誤差をδＨとし、正確な運動方程式の固有ベクトルｕｍ（ここでｍは複数存在し得る固有モードの識別子である。）、当該ｕｍに対応する固有周波数をωｍとして、δＴｍｍ及びδＨｍｍを固有ベクトルに対するδＴ及びδＨの行列要素の対角成分とした場合に、
前記Ｔ′及び前記Ｈ′が全ての前記固有モードについて関係式、
ωｍ２δＴｍｍ−δＨｍｍ＝０
を空間についての差分間隔の２次以上の精度で近似的に満たすという条件の下で、対角行列Ｘを選択して、次式、
【数６】
で表される前記δＴｎｕｍを生成すること、
を特徴とする演算子生成方法。
【請求項６】
予測子修正子法を用いた時間領域有限差分法により行う、弾性体からなる有限な媒質中の波動についての数値解析に関し、運動方程式の陰形式の差分近似である次式、
【数７】
（ここで、Δｔは時間方向の差分ステップを表し、ｃ＋，ｃ０，ｃ−は解析対象空間に設定した複数の離散点での変位ベクトルの各成分からなる状態ベクトルであって、それぞれ時刻ｔ＋Δｔ，ｔ，ｔ−Δｔにおけるものを表し、ｆは前記各離散点での外力ベクトルを表し、Ｔ′は質量行列を表し、Ｈ′は剛性行列を表す。）に基づいて設定される予測子及び修正子にて用いる演算子を生成する演算子生成装置であって、
当該演算子として前記Ｈ′、及び次式、
δＴｎｕｍ＝Ｔ′−Ｔ
で表されるδＴｎｕｍ（ここで、Ｔは前記離散点での差分近似により導かれる正定値の対角質量行列である。）を生成する手段を有し、
前記Ｈ′を、前記状態ベクトルに作用し前記各離散点での前記変位ベクトルの空間１階差分演算子Ｐｉ（ここで、ｉは、Ｉを前記解析対象空間の次元に応じて決まる所定の自然数として、１≦ｉ≦Ｉなる整数である。）と、弾性定数行列Ｃとから次式、
【数８】
（ここで、ａｉは正のスカラー係数であり、演算子ＰｉＴは演算子Ｐｉを表す行列の転置行列で表される演算子である。）に基づいて生成し、
回転及び並進に対応する前記状態ベクトルに前記演算子δＴｎｕｍを作用させた結果が０となるように当該δＴｎｕｍを設定すること、
を特徴とする演算子生成装置。
【請求項７】
請求項６に記載の演算子生成装置において、
前記解析対象空間はｘ軸、ｙ軸及びｚ軸を座標軸とする３次元空間であって、
前記各離散点に対応して設定される微分定義位置でのｘ軸方向、ｙ軸方向及びｚ軸方向それぞれに関する２点１階差分演算子をＤｘ，Ｄｙ，Ｄｚとし、
前記微分定義位置でのｘ軸方向、ｙ軸方向及びｚ軸方向それぞれに関する３点差分近似による、最適精度を持つように設計した２次精度０階差分演算子をＥｘｒ，Ｅｙｓ，Ｅｚｔ（ここでｒ，ｓ，ｔは１又は２で、差分に用いる３つの前記離散点の前記微分定義位置に対する配置の相違に応じた当該演算子の種類を示す。）とし、
次式、
【数９】
により、前記Ｄｘ，Ｄｙ，Ｄｚを微分方向の軸とは別の他の２つの軸方向に拡散した１階差分演算子Ｄｘｒｓｔ，Ｄｙｒｓｔ，Ｄｚｒｓｔを前記各微分定義位置にて定義した場合に、
前記各演算子Ｐｉ（ここで、前記Ｉは８であり、ｉは１≦ｉ≦Ｉなる整数である。）は、前記Ｄｘｒｓｔ，Ｄｙｒｓｔ，Ｄｚｒｓｔを用いて前記ｒ，ｓ，ｔの８通りの組み合わせのそれぞれに対応して定められ、前記各ａｉは１／８であること、
を特徴とする演算子生成装置。
【請求項８】
請求項６に記載の演算子生成装置において、
前記解析対象空間はｘ軸及びｚ軸を座標軸とする２次元空間であって、
前記各離散点に対応して設定される微分定義位置でのｘ軸方向及びｚ軸方向それぞれに関する２点１階差分演算子をＤｘ，Ｄｚとし、
前記微分定義位置でのｘ軸方向及びｚ軸方向それぞれに関する３点差分近似による、最適精度を持つように設計した２次精度０階差分演算子をＥｘｒ，Ｅｚｓ（ここでｒ，ｓは１又は２で、差分に用いる３つの前記離散点の前記微分定義位置に対する配置の相違に応じた当該演算子の種類を示す。）とし、
次式、
【数１０】
により、前記Ｄｘ，Ｄｚを微分方向の軸とは別の軸方向に拡散した１階差分演算子Ｄｘｒｓ，Ｄｚｒｓを前記各微分定義位置にて定義した場合に、
前記各演算子Ｐｉ（ここで、前記Ｉは４であり、ｉは１≦ｉ≦Ｉなる整数である。）は、前記Ｄｘｒｓ，Ｄｚｒｓを用いて前記ｒ，ｓの４通りの組み合わせのそれぞれに対応して定められ、前記各ａｉは１／４であること、
を特徴とする演算子生成装置。
【請求項９】
請求項７又は請求項８に記載の演算子生成装置において、
前記Ｈ′を、前記状態ベクトルに作用し前記各離散点での前記変位ベクトルの空間２階差分演算子Ｑｊ（ここで、ｊは、Ｊを前記解析対象空間の次元に応じて決まる所定の自然数として、１≦ｊ≦Ｊなる整数である。）を用いて修正した次式、
【数１１】
（ここで、εｊは正のスカラー係数であり、演算子ＱｊＴは演算子Ｑｊを表す行列の転置行列で表される演算子である。）に基づいて生成し、
前記εｊは、当該Ｈ′の前記微分定義位置におけるテイラー展開の４次の項の係数が０以上となるように定められること、
を特徴とする演算子生成装置。
【請求項１０】
請求項９に記載の演算子生成装置において、
前記質量行列を表す正確な演算子Ｔ０に対する前記Ｔ′の誤差をδＴとし、前記剛性行列を表す正確な演算子Ｈ０に対する前記Ｈ′の誤差をδＨとし、正確な運動方程式の固有ベクトルｕｍ（ここでｍは複数存在し得る固有モードの識別子である。）、当該ｕｍに対応する固有周波数をωｍとして、δＴｍｍ及びδＨｍｍを固有ベクトルに対するδＴ及びδＨの行列要素の対角成分とした場合に、
前記Ｔ′及び前記Ｈ′が全ての前記固有モードについて関係式、
ωｍ２δＴｍｍ−δＨｍｍ＝０
を空間についての差分間隔の２次以上の精度で近似的に満たすという条件の下で、対角行列Ｘを選択して、次式、
【数１２】
で表される前記δＴｎｕｍを生成すること、
を特徴とする演算子生成装置。
【請求項１１】
予測子修正子法を用いた時間領域有限差分法により行う、弾性体からなる有限な媒質中の波動についての数値解析を行うシミュレーション装置であって、
運動方程式の陰形式の差分近似である次式、
【数１３】
（ここで、Δｔは時間方向の差分ステップを表し、ｃ＋，ｃ０，ｃ−は解析対象空間に設定した複数の離散点での変位ベクトルの各成分からなる状態ベクトルであって、それぞれ時刻ｔ＋Δｔ，ｔ，ｔ−Δｔにおけるものを表し、ｆは前記各離散点での外力ベクトルを表し、Ｔ′は質量行列を表し、Ｈ′は剛性行列を表す。）に基づいて設定される予測子及び修正子であって、前記離散点での差分近似により導かれる正定値の対角質量行列Ｔ及び、δＴｎｕｍ＝Ｔ′−Ｔで表される演算子δＴｎｕｍを用いて、次式、
【数１４】
で与えられる前記予測子と、次式、
【数１５】
で与えられる前記修正子とをそれぞれ計算して、次式、
【数１６】
で与えられる前記状態ベクトルｃ＋を時間発展的に順次生成する演算手段を有し、
前記Ｈ′は、前記状態ベクトルに作用し前記各離散点での前記変位ベクトルの空間１階差分演算子Ｐｉ（ここで、ｉは、Ｉを前記解析対象空間の次元に応じて決まる所定の自然数として、１≦ｉ≦Ｉなる整数である。）と、弾性定数行列Ｃとから次式、
【数１７】
（ここで、ａｉは正のスカラー係数であり、演算子ＰｉＴは演算子Ｐｉを表す行列の転置行列で表される演算子である。）で表される演算子であり、
前記演算子δＴｎｕｍは、回転及び並進に対応する前記状態ベクトルに当該演算子δＴｎｕｍを作用させた結果が０となる演算子であること、
を特徴とするシミュレーション装置。
【請求項１２】
請求項１１に記載のシミュレーション装置において、
前記解析対象空間はｘ軸、ｙ軸及びｚ軸を座標軸とする３次元空間であって、
前記各離散点に対応して設定される微分定義位置でのｘ軸方向、ｙ軸方向及びｚ軸方向それぞれに関する２点１階差分演算子をＤｘ，Ｄｙ，Ｄｚとし、
前記微分定義位置でのｘ軸方向、ｙ軸方向及びｚ軸方向それぞれに関する３点差分近似による、最適精度を持つように設計した２次精度０階差分演算子をＥｘｒ，Ｅｙｓ，Ｅｚｔ（ここでｒ，ｓ，ｔは１又は２で、差分に用いる３つの前記離散点の前記微分定義位置に対する配置の相違に応じた当該演算子の種類を示す。）とし、
次式、
【数１８】
により、前記Ｄｘ，Ｄｙ，Ｄｚを微分方向の軸とは別の他の２つの軸方向に拡散した１階差分演算子Ｄｘｒｓｔ，Ｄｙｒｓｔ，Ｄｚｒｓｔを前記各微分定義位置にて定義した場合に、
前記Ｈ′は、前記Ｄｘｒｓｔ，Ｄｙｒｓｔ，Ｄｚｒｓｔを用いて前記ｒ，ｓ，ｔの８通りの組み合わせのそれぞれに対応して定められた前記演算子Ｐｉ（ここで、前記Ｉは８であり、ｉは１≦ｉ≦Ｉなる整数である。）を用い、前記各ａｉを１／８として生成されるものであること、
を特徴とするシミュレーション装置。
【請求項１３】
請求項１１に記載のシミュレーション装置において、
前記解析対象空間はｘ軸及びｚ軸を座標軸とする２次元空間であって、
前記各離散点に対応して設定される微分定義位置でのｘ軸方向及びｚ軸方向それぞれに関する２点１階差分演算子をＤｘ，Ｄｚとし、
前記微分定義位置でのｘ軸方向及びｚ軸方向それぞれに関する３点差分近似による、最適精度を持つように設計した２次精度０階差分演算子をＥｘｒ，Ｅｚｓ（ここでｒ，ｓは１又は２で、差分に用いる３つの前記離散点の前記微分定義位置に対する配置の相違に応じた当該演算子の種類を示す。）とし、
次式、
【数１９】
により、前記Ｄｘ，Ｄｚを微分方向の軸とは別の軸方向に拡散した１階差分演算子Ｄｘｒｓ，Ｄｚｒｓを前記各微分定義位置にて定義した場合に、
前記Ｈ′は、前記Ｄｘｒｓ，Ｄｚｒｓを用いて前記ｒ，ｓの４通りの組み合わせのそれぞれに対応して定められた前記演算子Ｐｉ（ここで、前記Ｉは４であり、ｉは１≦ｉ≦Ｉなる整数である。）を用い、前記各ａｉを１／４として生成されるものであること、
を特徴とするシミュレーション装置。
【請求項１４】
請求項１２又は請求項１３に記載のシミュレーション装置において、
前記Ｈ′は、前記状態ベクトルに作用し前記各離散点での前記変位ベクトルの空間２階差分演算子Ｑｊ（ここで、ｊは、Ｊを前記解析対象空間の次元に応じて決まる所定の自然数として、１≦ｊ≦Ｊなる整数である。）を用いて修正した次式、
【数２０】
（ここで、εｊは正のスカラー係数であり、演算子ＱｊＴは演算子Ｑｊを表す行列の転置行列で表される演算子である。）で表される演算子であり、
前記εｊは、当該Ｈ′の前記微分定義位置におけるテイラー展開の４次の項の係数が０以上となるように定められていること、
を特徴とするシミュレーション装置。
【請求項１５】
請求項１４に記載のシミュレーション装置において、
前記質量行列を表す正確な演算子Ｔ０に対する前記Ｔ′の誤差をδＴとし、前記剛性行列を表す正確な演算子Ｈ０に対する前記Ｈ′の誤差をδＨとし、正確な運動方程式の固有ベクトルｕｍ（ここでｍは複数存在し得る固有モードの識別子である。）、当該ｕｍに対応する固有周波数をωｍとして、δＴｍｍ及びδＨｍｍを固有ベクトルに対するδＴ及びδＨの行列要素の対角成分とした場合に、
前記Ｔ′及び前記Ｈ′は、全ての前記固有モードについて関係式、
ωｍ２δＴｍｍ−δＨｍｍ＝０
を空間についての差分間隔の２次以上の精度で近似的に満たし、
前記δＴｎｕｍは、対角行列Ｘを用いて、次式、
【数２１】
で表されるものであること、
を特徴とするシミュレーション装置。
【請求項１６】
速度もしくは変位、及び応力を互いに食い違った格子上で定義するスタッガードグリッド法を用いた時間領域有限差分法により弾性体からなる有限な媒質中の波動についての数値解析を行うシミュレーション装置であって、
前記速度もしくは変位の空間１階差分演算子としてＰを用い、
前記応力の空間１階差分演算子として前記演算子Ｐを表す行列の転置行列で表される演算子ＰＴを用い、
解析対象空間における自由表面を、垂直速度成分を持つ格子上に配置すること、
を特徴とするシミュレーション装置。
【図１】
【図２】
【図２】
【公開番号】特開２０１０−１２３０５６（Ｐ２０１０−１２３０５６Ａ）
【公開日】平成２２年６月３日（２０１０．６．３）
【国際特許分類】
【出願番号】特願２００８−２９８２０７（Ｐ２００８−２９８２０７）
【出願日】平成２０年１１月２１日（２００８．１１．２１）
【出願人】（５０４１３７９１２）国立大学法人 東京大学 (1,942)
【Ｆターム（参考）】
【公開日】平成２２年６月３日（２０１０．６．３）
【国際特許分類】
【出願日】平成２０年１１月２１日（２００８．１１．２１）
【出願人】（５０４１３７９１２）国立大学法人 東京大学 (1,942)
【Ｆターム（参考）】
[ Back to top ]