JP5279016B2 - 地球内部を伝播する地震波についての数値解析に用いる演算子生成方法、及び演算子生成装置、並びに地球内部を伝播する地震波についての数値解析を行うシミュレーション装置 - Google Patents

Description

本発明は、時間領域有限差分法により行う、弾性体からなる有限な媒質中の波動についての数値解析に関し、その演算子生成方法、演算子生成装置、及びシミュレーション装置に関する。

弾性体の運動方程式を解く手法として時間領域差分法が用いられている。時間領域差分法を用いた数値解析では、安定性、計算効率、精度の３つが問題となる。古くから数多くの時間領域差分法による弾性波動計算手法は開発されてきたが、これら安定性、計算効率、精度を同時に好適に実現する手法はこれまでなかった。

例えば、非特許文献１に示される時間領域差分計算手法は、弾性波動を効率よく、高精度に計算する手法を提供したが、安定性は保証されていなかった。

ここでいう、「効率がよい」とは、特定の計算精度を得るための数値解析を行う計算機の使用リソースが少ないことを意味する。具体的には、計算時間が少なく、また当該計算に用いる差分演算子の幅が抑制される手法は、高い効率性であると評価される。また、計算手法における「精度」の評価に関しては、許容される効率性の範囲で実現できる「最適な精度」というものを考えることができる。以下、「最適精度を持つ演算子」とは非特許文献２に示される条件（後述の（９）式）を満たす演算子を意味する。
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) の手法に適用することが可能である。すなわち、変位および速度の微分にＰ（もしくはＰ^Ｔ ）を用い、応力の微分にはＰ^Ｔ（もしくはＰ）を用いて、時間発展計算を行うことで、安定性が保障されることになる。このような手法の場合、非特許文献３で示されているように、予測子修正子法は適用できないため、高精度性は保障できないが、予測子のみを用いて時間発展計算を行うことで、安定化の実現が可能である。

なお、一般に解析対象空間に設定される格子点の数は非常に大きく、例えば、行列形式で表される離散化した演算子の規模も巨大となる。そのため、当該巨大な演算子行列の全体を生成して記憶部に保持したり、記憶部に保持した巨大な行列と、巨大な状態ベクトルとを直接、乗算することは、計算効率の観点から得策ではない。この点については、既に１つの工夫を、図１及び図２を用いて説明した。このように、本発明に係るシミュレーション装置では、状態ベクトルに演算子行列を作用させる計算を、解析対象空間の一部分ずつ行い得る。そして、その場合には、演算子行列のうち、これから計算を行う一部分に対応した行列要素が記憶部上に存在すれば処理を実行することができる。本発明に係る演算子生成方法、演算子生成装置は、このように必要とする一部分ずつを逐次的にシミュレーション装置に提供する演算子生成の仕方を含んでいる。

また、演算子Ｈ′，δＴ^ｎｕｍは上述のように演算子（行列）の積を含む形式で表され得る。この場合に、演算子を生成するとは、その行列の積を計算した結果得られる行列Ｈ′や行列δＴ^ｎｕｍの行列要素を算出することには限定されず、積を構成する個々の演算子行列の要素を求めることも含む。例えば、図１及び図２を用いて説明した処理方法では、本発明の演算子生成では、Ｈ′を構成する行列の積Ｐ_ｉ ^ＴＣＰ_ｉの１階差分演算子Ｐ_ｉ ^Ｔ，Ｐ_ｉの行列要素は直接的・具体的に定めるが、Ｈ′の行列要素の具体的な値までは計算しない。しかし、この場合においてもＰ_ｉ ^Ｔ，Ｐ_ｉに基づいてＨ′は（４）式により定義され、この意味においてＨ′は生成されているということができる。

本発明に係る演算子生成方法、演算子生成装置、及びシミュレーション装置は、地震波の解析に用いることができる。地球内部を伝播してきた地震波形記録が持つ情報を十分に活用することができれば、より詳細で正確な内部構造推定が可能となる。高精度かつ高解像度の地球内部構造モデルの推定は、例えば、防災や、油田開発における物理探査などに利用することができる。

Ｐ^ＴＣＰの形式の行列の積と状態ベクトルｕとの積の演算方法を２次元の場合について説明するステンシルの模式図である。 Ｐ^ＴＣＰの形式の行列の積と状態ベクトルｕとの積の演算方法を２次元の場合について説明するステンシルの模式図である。

Claims (15)

1. 予測子修正子法を用いた時間領域有限差分法により行う、地球内部を伝播する地震波についての数値解析に関し、運動方程式の陰形式の差分近似である次式、
   

   

   （ここで、Δｔは時間方向の差分ステップを表し、ｃ_＋，ｃ_０，ｃ₋は前記地球内部の解析対象空間に設定した複数の離散点での変位ベクトルの各成分からなる状態ベクトルであって、それぞれ時刻ｔ＋Δｔ，ｔ，ｔ−Δｔにおけるものを表し、ｆは前記各離散点での外力ベクトルを表し、Ｔ′は前記各離散点での密度に応じた質量行列を表し、Ｈ′は剛性行列を表す。）に基づいて設定される予測子及び修正子にて用いる演算子を生成する演算子生成方法であって、
   当該演算子として前記Ｈ′、及び次式、
   δＴ^ｎｕｍ＝Ｔ′−Ｔ
   で表されるδＴ^ｎｕｍ（ここで、Ｔは前記離散点での差分近似により導かれる正定値の対角質量行列である。）を生成する方法を含み、
   前記Ｈ′を、前記状態ベクトルに作用し前記各離散点での前記変位ベクトルの空間１階差分演算子Ｐ_ｉ（ここで、ｉは、Ｉを前記解析対象空間の次元に応じて決まる所定の自然数として、１≦ｉ≦Ｉなる整数である。）と、前記各離散点での弾性定数に応じた弾性定数行列Ｃとから次式、
   

   

   （ここで、ａ_ｉは正のスカラー係数であり、演算子Ｐ_ｉ ^Ｔは演算子Ｐ_ｉを表す行列の転置行列で表される演算子である。）に基づいて生成し、
   回転及び並進に対応する前記状態ベクトルに前記演算子δＴ^ｎｕｍを作用させた結果が０となるように当該δＴ^ｎｕｍを設定すること、
   を特徴とする地球内部を伝播する地震波についての数値解析に用いる演算子生成方法。
2. 請求項１に記載の演算子生成方法において、
   前記解析対象空間はｘ軸、ｙ軸及びｚ軸を座標軸とする３次元空間であって、
   前記各離散点に対応して設定される微分定義位置でのｘ軸方向、ｙ軸方向及びｚ軸方向それぞれに関する２点１階差分演算子をＤ_ｘ，Ｄ_ｙ，Ｄ_ｚとし、
   前記微分定義位置でのｘ軸方向、ｙ軸方向及びｚ軸方向それぞれに関する３点差分近似による、最適精度を持つように設計した２次精度０階差分演算子をＥ_ｘ ^ｒ，Ｅ_ｙ ^ｓ，Ｅ_ｚ ^ｔ（ここでｒ，ｓ，ｔは１又は２で、差分に用いる３つの前記離散点の前記微分定義位置に対する配置の相違に応じた当該演算子の種類を示す。）とし、
   次式、
   

   

   により、前記Ｄ_ｘ，Ｄ_ｙ，Ｄ_ｚを微分方向の軸とは別の他の２つの軸方向に拡散した１階差分演算子Ｄ_ｘ ^ｒｓｔ，Ｄ_ｙ ^ｒｓｔ，Ｄ_ｚ ^ｒｓｔを前記各微分定義位置にて定義した場合に、
   前記各演算子Ｐ_ｉ（ここで、前記Ｉは８であり、ｉは１≦ｉ≦Ｉなる整数である。）は、前記Ｄ_ｘ ^ｒｓｔ，Ｄ_ｙ ^ｒｓｔ，Ｄ_ｚ ^ｒｓｔを用いて前記ｒ，ｓ，ｔの８通りの組み合わせのそれぞれに対応して定められ、前記各ａ_ｉは１／８であること、
   を特徴とする演算子生成方法。
3. 請求項１に記載の演算子生成方法において、
   前記解析対象空間はｘ軸及びｚ軸を座標軸とする２次元空間であって、
   前記各離散点に対応して設定される微分定義位置でのｘ軸方向及びｚ軸方向それぞれに関する２点１階差分演算子をＤ_ｘ，Ｄ_ｚとし、
   前記微分定義位置でのｘ軸方向及びｚ軸方向それぞれに関する３点差分近似による、最適精度を持つように設計した２次精度０階差分演算子をＥ_ｘ ^ｒ，Ｅ_ｚ ^ｓ（ここでｒ，ｓは１又は２で、差分に用いる３つの前記離散点の前記微分定義位置に対する配置の相違に応じた当該演算子の種類を示す。）とし、
   次式、
   

   

   により、前記Ｄ_ｘ，Ｄ_ｚを微分方向の軸とは別の軸方向に拡散した１階差分演算子Ｄ_ｘ ^ｒｓ，Ｄ_ｚ ^ｒｓを前記各微分定義位置にて定義した場合に、
   前記各演算子Ｐ_ｉ（ここで、前記Ｉは４であり、ｉは１≦ｉ≦Ｉなる整数である。）は、前記Ｄ_ｘ ^ｒｓ，Ｄ_ｚ ^ｒｓを用いて前記ｒ，ｓの４通りの組み合わせのそれぞれに対応して定められ、前記各ａ_ｉは１／４であること、
   を特徴とする演算子生成方法。
4. 請求項２又は請求項３に記載の演算子生成方法において、
   前記Ｈ′を、前記状態ベクトルに作用し前記各離散点での前記変位ベクトルの空間２階差分演算子Ｑ_ｊ（ここで、ｊは、Ｊを前記解析対象空間の次元に応じて決まる所定の自然数として、１≦ｊ≦Ｊなる整数である。）を用いて修正した次式、
   

   

   （ここで、ε_ｊは正のスカラー係数であり、演算子Ｑ_ｊ ^Ｔは演算子Ｑ_ｊを表す行列の転置行列で表される演算子である。）に基づいて生成し、
   前記ε_ｊは、当該Ｈ′の前記微分定義位置におけるテイラー展開の４次の項の係数が０以上となるように定められること、
   を特徴とする演算子生成方法。
5. 請求項４に記載の演算子生成方法において、
   前記質量行列を表す正確な演算子Ｔ^０に対する前記Ｔ′の誤差をδＴとし、前記剛性行列を表す正確な演算子Ｈ^０に対する前記Ｈ′の誤差をδＨとし、正確な運動方程式の固有ベクトルｕ_ｍ（ここでｍは複数存在し得る固有モードの識別子である。）、当該ｕ_ｍに対応する固有周波数をω_ｍとして、δＴ_ｍｍ及びδＨ_ｍｍを固有ベクトルに対するδＴ及びδＨの行列要素の対角成分とした場合に、
   前記Ｔ′及び前記Ｈ′が全ての前記固有モードについて関係式、
   ω_ｍ ^２δＴ_ｍｍ−δＨ_ｍｍ＝０
   を空間についての差分間隔の２次以上の精度で近似的に満たすという条件の下で、対角行列Ｘを選択して、次式、
   

   

   で表される前記δＴ^ｎｕｍを生成すること、
   を特徴とする演算子生成方法。
6. 予測子修正子法を用いた時間領域有限差分法により行う、地球内部を伝播する地震波についての数値解析に関し、運動方程式の陰形式の差分近似である次式、
   

   

   （ここで、Δｔは時間方向の差分ステップを表し、ｃ_＋，ｃ_０，ｃ₋は前記地球内部の解析対象空間に設定した複数の離散点での変位ベクトルの各成分からなる状態ベクトルであって、それぞれ時刻ｔ＋Δｔ，ｔ，ｔ−Δｔにおけるものを表し、ｆは前記各離散点での外力ベクトルを表し、Ｔ′は前記各離散点での密度に応じた質量行列を表し、Ｈ′は剛性行列を表す。）に基づいて設定される予測子及び修正子にて用いる演算子を生成する演算子生成装置であって、
   当該演算子として前記Ｈ′、及び次式、
   δＴ^ｎｕｍ＝Ｔ′−Ｔ
   で表されるδＴ^ｎｕｍ（ここで、Ｔは前記離散点での差分近似により導かれる正定値の対角質量行列である。）を生成する手段を有し、
   前記Ｈ′を、前記状態ベクトルに作用し前記各離散点での前記変位ベクトルの空間１階差分演算子Ｐ_ｉ（ここで、ｉは、Ｉを前記解析対象空間の次元に応じて決まる所定の自然数として、１≦ｉ≦Ｉなる整数である。）と、前記各離散点での弾性定数に応じた弾性定数行列Ｃとから次式、
   

   

   （ここで、ａ_ｉは正のスカラー係数であり、演算子Ｐ_ｉ ^Ｔは演算子Ｐ_ｉを表す行列の転置行列で表される演算子である。）に基づいて生成し、
   回転及び並進に対応する前記状態ベクトルに前記演算子δＴ^ｎｕｍを作用させた結果が０となるように当該δＴ^ｎｕｍを設定すること、
   を特徴とする地球内部を伝播する地震波についての数値解析に用いる演算子生成装置。
7. 請求項６に記載の演算子生成装置において、
   前記解析対象空間はｘ軸、ｙ軸及びｚ軸を座標軸とする３次元空間であって、
   前記各離散点に対応して設定される微分定義位置でのｘ軸方向、ｙ軸方向及びｚ軸方向それぞれに関する２点１階差分演算子をＤ_ｘ，Ｄ_ｙ，Ｄ_ｚとし、
   前記微分定義位置でのｘ軸方向、ｙ軸方向及びｚ軸方向それぞれに関する３点差分近似による、最適精度を持つように設計した２次精度０階差分演算子をＥ_ｘ ^ｒ，Ｅ_ｙ ^ｓ，Ｅ_ｚ ^ｔ（ここでｒ，ｓ，ｔは１又は２で、差分に用いる３つの前記離散点の前記微分定義位置に対する配置の相違に応じた当該演算子の種類を示す。）とし、
   次式、
   

   

   により、前記Ｄ_ｘ，Ｄ_ｙ，Ｄ_ｚを微分方向の軸とは別の他の２つの軸方向に拡散した１階差分演算子Ｄ_ｘ ^ｒｓｔ，Ｄ_ｙ ^ｒｓｔ，Ｄ_ｚ ^ｒｓｔを前記各微分定義位置にて定義した場合に、
   前記各演算子Ｐ_ｉ（ここで、前記Ｉは８であり、ｉは１≦ｉ≦Ｉなる整数である。）は、前記Ｄ_ｘ ^ｒｓｔ，Ｄ_ｙ ^ｒｓｔ，Ｄ_ｚ ^ｒｓｔを用いて前記ｒ，ｓ，ｔの８通りの組み合わせのそれぞれに対応して定められ、前記各ａ_ｉは１／８であること、
   を特徴とする演算子生成装置。
8. 請求項６に記載の演算子生成装置において、
   前記解析対象空間はｘ軸及びｚ軸を座標軸とする２次元空間であって、
   前記各離散点に対応して設定される微分定義位置でのｘ軸方向及びｚ軸方向それぞれに関する２点１階差分演算子をＤ_ｘ，Ｄ_ｚとし、
   前記微分定義位置でのｘ軸方向及びｚ軸方向それぞれに関する３点差分近似による、最適精度を持つように設計した２次精度０階差分演算子をＥ_ｘ ^ｒ，Ｅ_ｚ ^ｓ（ここでｒ，ｓは１又は２で、差分に用いる３つの前記離散点の前記微分定義位置に対する配置の相違に応じた当該演算子の種類を示す。）とし、
   次式、
   

   

   により、前記Ｄ_ｘ，Ｄ_ｚを微分方向の軸とは別の軸方向に拡散した１階差分演算子Ｄ_ｘ ^ｒｓ，Ｄ_ｚ ^ｒｓを前記各微分定義位置にて定義した場合に、
   前記各演算子Ｐ_ｉ（ここで、前記Ｉは４であり、ｉは１≦ｉ≦Ｉなる整数である。）は、前記Ｄ_ｘ ^ｒｓ，Ｄ_ｚ ^ｒｓを用いて前記ｒ，ｓの４通りの組み合わせのそれぞれに対応して定められ、前記各ａ_ｉは１／４であること、
   を特徴とする演算子生成装置。
9. 請求項７又は請求項８に記載の演算子生成装置において、
   前記Ｈ′を、前記状態ベクトルに作用し前記各離散点での前記変位ベクトルの空間２階差分演算子Ｑ_ｊ（ここで、ｊは、Ｊを前記解析対象空間の次元に応じて決まる所定の自然数として、１≦ｊ≦Ｊなる整数である。）を用いて修正した次式、
   

   

   （ここで、ε_ｊは正のスカラー係数であり、演算子Ｑ_ｊ ^Ｔは演算子Ｑ_ｊを表す行列の転置行列で表される演算子である。）に基づいて生成し、
   前記ε_ｊは、当該Ｈ′の前記微分定義位置におけるテイラー展開の４次の項の係数が０以上となるように定められること、
   を特徴とする演算子生成装置。
10. 請求項９に記載の演算子生成装置において、
    前記質量行列を表す正確な演算子Ｔ^０に対する前記Ｔ′の誤差をδＴとし、前記剛性行列を表す正確な演算子Ｈ^０に対する前記Ｈ′の誤差をδＨとし、正確な運動方程式の固有ベクトルｕ_ｍ（ここでｍは複数存在し得る固有モードの識別子である。）、当該ｕ_ｍに対応する固有周波数をω_ｍとして、δＴ_ｍｍ及びδＨ_ｍｍを固有ベクトルに対するδＴ及びδＨの行列要素の対角成分とした場合に、
    前記Ｔ′及び前記Ｈ′が全ての前記固有モードについて関係式、
    ω_ｍ ^２δＴ_ｍｍ−δＨ_ｍｍ＝０
    を空間についての差分間隔の２次以上の精度で近似的に満たすという条件の下で、対角行列Ｘを選択して、次式、
    

    

    で表される前記δＴ^ｎｕｍを生成すること、
    を特徴とする演算子生成装置。
11. 予測子修正子法を用いた時間領域有限差分法により行う、地球内部を伝播する地震波についての数値解析を行うシミュレーション装置であって、
    運動方程式の陰形式の差分近似である次式、
    

    

    （ここで、Δｔは時間方向の差分ステップを表し、ｃ_＋，ｃ_０，ｃ₋は前記地球内部の解析対象空間に設定した複数の離散点での変位ベクトルの各成分からなる状態ベクトルであって、それぞれ時刻ｔ＋Δｔ，ｔ，ｔ−Δｔにおけるものを表し、ｆは前記各離散点での外力ベクトルを表し、Ｔ′は前記各離散点での密度に応じた質量行列を表し、Ｈ′は剛性行列を表す。）に基づいて設定される予測子及び修正子であって、前記離散点での差分近似により導かれる正定値の対角質量行列Ｔ及び、δＴ^ｎｕｍ＝Ｔ′−Ｔで表される演算子δＴ^ｎｕｍを用いて、次式、
    

    

    で与えられる前記予測子と、次式、
    

    

    で与えられる前記修正子とをそれぞれ計算して、次式、
    

    

    で与えられる前記状態ベクトルｃ_＋を時間発展的に順次生成する演算手段を有し、
    前記Ｈ′は、前記状態ベクトルに作用し前記各離散点での前記変位ベクトルの空間１階差分演算子Ｐ_ｉ（ここで、ｉは、Ｉを前記解析対象空間の次元に応じて決まる所定の自然数として、１≦ｉ≦Ｉなる整数である。）と、前記各離散点での弾性定数に応じた弾性定数行列Ｃとから次式、
    

    

    （ここで、ａ_ｉは正のスカラー係数であり、演算子Ｐ_ｉ ^Ｔは演算子Ｐ_ｉを表す行列の転置行列で表される演算子である。）で表される演算子であり、
    前記演算子δＴ^ｎｕｍは、回転及び並進に対応する前記状態ベクトルに当該演算子δＴ^ｎｕｍを作用させた結果が０となる演算子であること、
    を特徴とするシミュレーション装置。
12. 請求項１１に記載のシミュレーション装置において、
    前記解析対象空間はｘ軸、ｙ軸及びｚ軸を座標軸とする３次元空間であって、
    前記各離散点に対応して設定される微分定義位置でのｘ軸方向、ｙ軸方向及びｚ軸方向それぞれに関する２点１階差分演算子をＤ_ｘ，Ｄ_ｙ，Ｄ_ｚとし、
    前記微分定義位置でのｘ軸方向、ｙ軸方向及びｚ軸方向それぞれに関する３点差分近似による、最適精度を持つように設計した２次精度０階差分演算子をＥ_ｘ ^ｒ，Ｅ_ｙ ^ｓ，Ｅ_ｚ ^ｔ（ここでｒ，ｓ，ｔは１又は２で、差分に用いる３つの前記離散点の前記微分定義位置に対する配置の相違に応じた当該演算子の種類を示す。）とし、
    次式、
    

    

    により、前記Ｄ_ｘ，Ｄ_ｙ，Ｄ_ｚを微分方向の軸とは別の他の２つの軸方向に拡散した１階差分演算子Ｄ_ｘ ^ｒｓｔ，Ｄ_ｙ ^ｒｓｔ，Ｄ_ｚ ^ｒｓｔを前記各微分定義位置にて定義した場合に、
    前記Ｈ′は、前記Ｄ_ｘ ^ｒｓｔ，Ｄ_ｙ ^ｒｓｔ，Ｄ_ｚ ^ｒｓｔを用いて前記ｒ，ｓ，ｔの８通りの組み合わせのそれぞれに対応して定められた前記演算子Ｐ_ｉ（ここで、前記Ｉは８であり、ｉは１≦ｉ≦Ｉなる整数である。）を用い、前記各ａ_ｉを１／８として生成されるものであること、
    を特徴とするシミュレーション装置。
13. 請求項１１に記載のシミュレーション装置において、
    前記解析対象空間はｘ軸及びｚ軸を座標軸とする２次元空間であって、
    前記各離散点に対応して設定される微分定義位置でのｘ軸方向及びｚ軸方向それぞれに関する２点１階差分演算子をＤ_ｘ，Ｄ_ｚとし、
    前記微分定義位置でのｘ軸方向及びｚ軸方向それぞれに関する３点差分近似による、最適精度を持つように設計した２次精度０階差分演算子をＥ_ｘ ^ｒ，Ｅ_ｚ ^ｓ（ここでｒ，ｓは１又は２で、差分に用いる３つの前記離散点の前記微分定義位置に対する配置の相違に応じた当該演算子の種類を示す。）とし、
    次式、
    

    

    により、前記Ｄ_ｘ，Ｄ_ｚを微分方向の軸とは別の軸方向に拡散した１階差分演算子Ｄ_ｘ ^ｒｓ，Ｄ_ｚ ^ｒｓを前記各微分定義位置にて定義した場合に、
    前記Ｈ′は、前記Ｄ_ｘ ^ｒｓ，Ｄ_ｚ ^ｒｓを用いて前記ｒ，ｓの４通りの組み合わせのそれぞれに対応して定められた前記演算子Ｐ_ｉ（ここで、前記Ｉは４であり、ｉは１≦ｉ≦Ｉなる整数である。）を用い、前記各ａ_ｉを１／４として生成されるものであること、
    を特徴とするシミュレーション装置。
14. 請求項１２又は請求項１３に記載のシミュレーション装置において、
    前記Ｈ′は、前記状態ベクトルに作用し前記各離散点での前記変位ベクトルの空間２階差分演算子Ｑ_ｊ（ここで、ｊは、Ｊを前記解析対象空間の次元に応じて決まる所定の自然数として、１≦ｊ≦Ｊなる整数である。）を用いて修正した次式、
    

    

    （ここで、ε_ｊは正のスカラー係数であり、演算子Ｑ_ｊ ^Ｔは演算子Ｑ_ｊを表す行列の転置行列で表される演算子である。）で表される演算子であり、
    前記ε_ｊは、当該Ｈ′の前記微分定義位置におけるテイラー展開の４次の項の係数が０以上となるように定められていること、
    を特徴とするシミュレーション装置。
15. 請求項１４に記載のシミュレーション装置において、
    前記質量行列を表す正確な演算子Ｔ^０に対する前記Ｔ′の誤差をδＴとし、前記剛性行列を表す正確な演算子Ｈ^０に対する前記Ｈ′の誤差をδＨとし、正確な運動方程式の固有ベクトルｕ_ｍ（ここでｍは複数存在し得る固有モードの識別子である。）、当該ｕ_ｍに対応する固有周波数をω_ｍとして、δＴ_ｍｍ及びδＨ_ｍｍを固有ベクトルに対するδＴ及びδＨの行列要素の対角成分とした場合に、
    前記Ｔ′及び前記Ｈ′は、全ての前記固有モードについて関係式、
    ω_ｍ ^２δＴ_ｍｍ−δＨ_ｍｍ＝０
    を空間についての差分間隔の２次以上の精度で近似的に満たし、
    前記δＴ^ｎｕｍは、対角行列Ｘを用いて、次式、
    

    

    で表されるものであること、
    を特徴とするシミュレーション装置。

Images