(57)【要約】      （修正有）【課題】対象系の地震観測記録，常時微動観測記録あるいは振動試験記録などを使用して系の動的特性を示す固有値および固有ベクトルを演算の安定性および収束性を確保しながら精度良く算出する。【解決手段】時間領域のモード解析手法は，対象系が運動方程式に従うことを仮定して系の観測記録の真値を同式の解で表現し，その解に対して非線形の関係にある固有値は微小修正量について線形化した式を，線形の関係にある固有ベクトルは運動方程式の解をそれぞれ回帰式として使用して，最小二乗法に従って正規方程式を作成し，これを解いて算出される精度の高い安定した収束演算の手段を有する。【選択図】図１
2014021956.tif 000002
【特許請求の範囲】
【請求項１】
自由度Ｎの解析対象系の運動が運動方程式３９に従うことを仮定して，系の観測記録の真値を同式の解である式４０によって表現し，式中の固有値λｒおよび固有ベクトルｕｒをパラメータとする最小二乗法に基づいて定式化した時間領域のモード解析手法．
【数３９】
2014021956.tif 000057
ここに，Ｍ，ＣおよびＫは，それぞれ質量，減衰，剛性の（Ｎ×Ｎ）行列であり，ｘ（ｔ），ｆ（ｔ）は相対変位，外力を意味するＮ元ベクトルである．また，・は時間ｔに関する微分を意味する．
【数４０】
2014021956.tif 000058
ここに，λｒ，ｕｒは，第ｒ次固有値および固有ベクトルであり，ηｒ（ｔ）＝ｅｘｐ（λｒｔ）である．また，αｒｏは相対変位ｘ（ｔ）を固有ベクトルの線形展開式に表したときの展開係数αｒ（ｔ）の初期値であり，αｒ（ｔ）は式４１の微分方程式を満たす．
【数４１】
2014021956.tif 000059
また，ベクトルβｒ（ｔ）は式４２の微分方程式を満たす．
【数４２】
2014021956.tif 000060
尚，式４０は，モード定数λｒ，ｕｒの共役複素数を考慮すれば式４３のようになる．
【数４３】
2014021956.tif 000061
ここに，Ｒｅａｌは実数部を意味する．
【請求項２】
2014021956.tif 000062
い，最小二乗法により式４４を回帰式として求められる請求項１に記載の時間領域のモード解析手法．
【数４４】
2014021956.tif 000063
2014021956.tif 000064
は時間ｔ＝ｔｉを意味する．αｒｏ，ｕｒｊを含む新たな未知変数ａｒｊ，Ｕｒｊは，上式の理論値と観測値間の差の二乗和を評価して求められる．このとき誤差評価関数は，全ての観測データが等精度の場合は式４５のようになる．
【数４５】
2014021956.tif 000065
2014021956.tif 000066
式４６のεｊ最小化必要条件式より正規方程式を立て，これを解いて求められる．
【数４６】
2014021956.tif 000067
【請求項３】
2014021956.tif 000068
2014021956.tif 000069
化した式４７を回帰式として最小二乗法により求められる請求項１に記載の時間領域のモード解析手法．
【数４７】
2014021956.tif 000070
2014021956.tif 000071
の理論値である．未知変数Δλｒは，式４７の値と観測値間の差の二乗和を評価して求められる．このとき，誤差評価関数εは式４８のようになる．
【数４８】
2014021956.tif 000072
ここに，記号Σは観測点および解析対象時間内の総和である．Δλｒは式４９のεの最小化必要条件式より正規方程式を立て，これを解いて求められる．そして，固有値λｒはΔλｒが加えられ更新される．
【数４９】
2014021956.tif 000073
【請求項４】
2014021956.tif 000074
び相対変位ｘがそれぞれ式５０，５１より求められ，また，対象系のｊ点の伝達関数Ｔｊ（ω）が，入力が加速度（慣性力）の場合は式５２より，単位力加振の場合は式５３より，各々求められる請求項１に記載の時間領域のモード解析手法．
【数５０】
2014021956.tif 000075
【数５１】
2014021956.tif 000076
【数５２】
2014021956.tif 000077
【数５３】
2014021956.tif 000078
ここに，Ｉは単位ベクトルである．ところで，刺激係数ｕｒＴＭＩは系の定数であり，Ｍ単独では任意の値を採り得る．つまり，本解析から算出されるＵｒｊは適当なＭに応じた値となるが，刺激関数ＵｒｊＴＭＩ（＝ｕｒＴＭＩｕｒｊ）はＭに依らず一定である．また，ｅｋは観測点ｋに単位力が働いた場合，第ｋ要素に１．０，他の要素に０．０を持つベクトルである．加振点が複数の場合は，対応する要素に１．０を与えればよい．
【請求項５】
自由度Ｎの解析対象系の運動が運動力程式５４に従うことを仮定して，系の観測記録の真値を同式の解である式５５によって表現し，式中の固有値λｒおよび固有ベクトルｕｒをパラメータとする最小二乗法に基づいて定式化した時間領域のモード解析プログラム．
【数５４】
2014021956.tif 000079
ここに，Ｍ，ＣおよびＫは，それぞれ質量，減衰，剛性の（Ｎ×Ｎ）行列であり，ｘ（ｔ），ｆ（ｔ）は相対変位，外力を意味するＮ元ベクトルである．また，・は時間ｔに関する微分を意味する．
【数５５】
2014021956.tif 000080
ここに，λｒ，ｕｒは，第ｒ次固有値および固有ベクトルであり，ηｒ（ｔ）＝ｅｘｐ（λｒｔ）である．また，αｒｏは相対変位ｘ（ｔ）を固有ベクトルの線形展開式に表したときの展開係数αｒ（ｔ）の初期値であり，αｒ（ｔ）は式５６の微分方程式を満たす．
【数５６】
2014021956.tif 000081
また，ベクトルβｒ（ｔ）は式５７の微分方程式を満たす．
【数５７】
2014021956.tif 000082
尚，式５５は，モード定数λｒ，ｕｒの共役複素数を考慮すれば式５８のようになる．
【数５８】
2014021956.tif 000083
ここに，Ｒｅａｌは実数部を意味する．
【請求項６】
2014021956.tif 000084
い，最小二乗法により式５９を回帰式として求められる請求項５に記載の時間領域のモード解析プログラム．
【数５９】
2014021956.tif 000085
2014021956.tif 000086
は時間ｔ＝ｔｉを意味する．αｒｏ，ｕｒｊを含む新たな未知変数ａｒｊ，Ｕｒｊは，上式の理論値と観測値間の差の二乗和を評価して求められる．このとき誤差評価関数は，全ての観測データが等精度の場合は式６０のようになる．
【数６０】
2014021956.tif 000087
2014021956.tif 000088
式６１のεｊの最小化必要条件式より正規方程式を立て，これを解いて求められる．
【数６１】
2014021956.tif 000089
【請求項７】
2014021956.tif 000090
化した式６２を回帰式として最小二乗法により求められる請求項５に記載の時間領域のモード解析プログラム．
【数６２】
2014021956.tif 000091
2014021956.tif 000092
2014021956.tif 000093
の理論値である．未知変数Δλｒは，式６２の値と観測値間の差の二乗和を評価して求められる．このとき，誤差評価関数εは式６３のようになる．
【数６３】
2014021956.tif 000094
ここに，記号Σは観測点および解析対象時間内の総和である．Δλｒは式６４のεの最小化必要条件式より正規方程式を立て，これを解いて求められる．そして，固有値λｒはΔλｒが加えられ更新される．
【数６４】
2014021956.tif 000095
【請求項８】
2014021956.tif 000096
び相対変位ｘがそれぞれ式６５，６６より求められ，また，対象系のｊ点の伝達関数Ｔｊ（ω）が，入力が加速度（慣性力）の場合は式６７より，単位力加振の場合は式６８より，各々求められる請求項５に記載の時間領域のモード解析プログラム．
【数６５】
2014021956.tif 000097
【数６６】
2014021956.tif 000098
【数６７】
2014021956.tif 000099
【数６８】
2014021956.tif 000100
ここに，Ｉは単位ベクトルである．ところで，刺激係数ｕｒＴＭＩは系の定数であり，Ｍ単独では任意の値を採り得る．つまり，本解析から算出されるＵｒｊは適当なＭに応じた値となるが，刺激関数ＵｒｊＴＭＩ（＝ｕｒＴＭＩｕｒｊ）はＭに依らず一定である．また，ｅｋは観測点ｋに単位力が働いた場合，第ｋ要素に１．０，他の要素に０．０を持つベクトルである．加振点が複数の場合は，対応する要素に１．０を与えればよい．
【請求項９】
自由度Ｎの解析対象系の運動が運動方程式６９に従うことを仮定して，系の観測記録の真値を同式の解である式７０によって表現し，式中の固有値λｒおよび固有ベクトルｕｒをパラメータとする最小二乗法に基づいて定式化した時間領域のモード解析プログラムを記録したコンピュータ読み取り可能な記録媒体．
【数６９】
2014021956.tif 000101
ここに，Ｍ，ＣおよびＫは，それぞれ質量，減衰，剛性の（Ｎ×Ｎ）行列であり，ｘ（ｔ），ｆ（ｔ）は相対変位，外力を意味するＮ元ベクトルである．また，・は時間ｔに関する微分を意味する．
【数７０】
2014021956.tif 000102
ここに，λｒ，ｕｒは，第ｒ次固有値および固有ベクトルであり，ηｒ（ｔ）＝ｅｘｐ（λｒｔ）である．また，αｒｏは相対変位ｘ（ｔ）を固有ベクトルの線形展開式に表したときの展開係数αｒ（ｔ）の初期値であり，αｒ（ｔ）は式７１の微分方程式を満たす．
【数７１】
2014021956.tif 000103
また，ベクトルβｒ（ｔ）は式７２の微分方程式を満たす．
【数７２】
2014021956.tif 000104
尚，式７０は，モード定数λｒ，ｕｒの共役複素数を考慮すれば式７３のようになる．
【数７３】
2014021956.tif 000105
ここに，Ｒｅａｌは実数部を意味する．
【請求項１０】
2014021956.tif 000106
い，最小二乗法により式７４を回帰式として求められる請求項９に記載の時間領域のモード解析プログラムを記録したコンピュータ読み取り可能な記録媒体．
【数７４】
2014021956.tif 000107
2014021956.tif 000108
は時間ｔ＝ｔｉを意味する．αｒｏ，ｕｒｊを含む新たな未知変数ａｒｊ，Ｕｒｊは，上式の理論値と観測値間の差の二乗和を評価して求められる．このとき誤差評価関数は，全ての観測データが等精度の場合は式７５のようになる．
【数７５】
2014021956.tif 000109
2014021956.tif 000110
式７６のεｊの最小化必要条件式より正規方程式を立て，これを解いて求められる．
【数７６】
2014021956.tif 000111
【請求項１１】
2014021956.tif 000112
化した式７７を回帰式として最小二乗法により求められる請求項９に記載の時間領域のモード解析プログラムを記録したコンピュータ読み取り可能な記録媒体．
【数７７】
2014021956.tif 000113
2014021956.tif 000114
の理論値である．未知変数Δλｒは，式７７の値と観測値間の差の二乗和を評価して求められる．このとき，誤差評価関数εは式７８のようになる．
【数７８】
2014021956.tif 000115
ここに，記号Σは観測点および解析対象時間内の総和である．Δλｒは式７９のεの最小化必要条件式より正規方程式を立て，これを解いて求められる．そして，固有値λｒはΔλｒが加えられ更新される．
【数７９】
2014021956.tif 000116
【請求項１２】
2014021956.tif 000117
び相対変位ｘがそれぞれ式８０，８１より求められ，また，対象系のｊ点の伝達関数Ｔｊ（ω）が，入力が加速度（慣性力）の場合は式８２より，単位力加振の場合は式８３より，各々求められる請求項９に記載の時間領域のモード解析プログラムを記録したコンピュータ読み取り可能な記録媒体．
【数８０】
2014021956.tif 000118
【数８１】
2014021956.tif 000119
【数８２】
2014021956.tif 000120
【数８３】
2014021956.tif 000121
ここに，Ｉは単位ベクトルである．ところで，刺激係数ｕｒＴＭＩは系の定数であり，Ｍ単独では任意の値を採り得る．つまり，本解析から算出されるＵｒｊは適当なＭに応じた値となるが，刺激関数ＵｒｊＴＭＩ（＝ｕｒＴＭＩｕｒｊ）はＭに依らず一定である．また，ｅｋは観測点ｋに単位力が働いた場合，第ｋ要素に１．０，他の要素に０．０を持つベクトルである．加振点が複数の場合は，対応する要素に１．０を与えればよい．
【発明の詳細な説明】
【技術分野】
【０００１】  本発明は，土木，建築および機械系構造物の動的特性を把握するための技術に係り，対象構造物の地震観測記録，常時微動観測記録あるいは振動試験記録を使用して構造物の動的特性を示す固有値および固有ベクトルを求めるための時間領域のモード解析に関する．
【背景技術】
【０００２】  時間領域のモード解析手法には，ベック手法（非特許文献１）がある．同手法は，対象系の観測加速度記録を使用して，系のモード定数，即ち，固有値および固有ベクトルを次数（モード）毎に推定する方法である．自由度Ｎの運動方程式は式１のように表されるが，対象系がこの式に従うものと仮定して，ベック手法について述べていく．
【数１】
2014021956.tif 000003
ここに，Ｍ，ＣおよびＫは，それぞれ質量，減衰，剛性の（Ｎ×Ｎ）行列であり，ｘ（ｔ），ｆ（ｔ丿は相対変位，外力を意味するＮ元ベクトルである．入力が加速度であれば外力ｆ（ｔ）は慣性力となる．また，・は時間ｔに関する微分を意味する．【０００３】
2014021956.tif 000004
ベクトルｕｒの線形展開式に表される．
【数２】
2014021956.tif 000005
ここに，λｒ，ｕｒは第ｒ次固有値，固有ベクトルであり，αｒ（ｔ）は相対変位ｘ（ｔ）をｕｒの線形展開式に表したときの展開係数である．また，Ｒｅａｌは実数部を，Σは全次数についての和を意味する．【０００４】  時間領域のモード解析は逆解析であり，観測加速度記録の誤差分布を正規分布と仮定して最小二乗法に基づいて検討する場合は，式２を回帰式に使用して全次数の固有値および固有ベクトルを同時にパラメータとして処理しなければならない．ところが，ベック手法は，上記のように，特定の次数，例えば，第ｓ次モードを対象とした場合，固有値λｓおよび固有ベクトルｕｓをパラメータとして扱い，式３を観測方程式とする．
【数３】
2014021956.tif 000006
2014021956.tif 000007
他のモードの和を意味する．【０００５】
2014021956.tif 000008
形の関係にあるが，ベック手法では最小二乗法に従って解いているのは固有ベクトルのみである．固有値λｓは，式４の関係式を利用して固有周波数ｆｓと固有減衰定数ｈｓを試行的に現在値の周辺で微小変化させ，同時に誤差の二乗和，即ち，式３の両辺の差の時間ｔと全
2014021956.tif 000009
これら二変数によって求められる固有値をその最確値とする．
【数４】
2014021956.tif 000010
ここに，ｉは虚数単位である．【０００６】  解析では，λｓ，ｕｓを繰り返し交互に計算して，式３の二乗誤差の収束を図る．そして，この演算を全次数について実施して，最終的に収束して得られた固有値λｓ，固有ベクトルｕｓを系のモード定数とする．
【先行技術文献】
【非特許文献】
【０００７】
【非特許文献１】Ｂｅｃｋ，Ｊ．Ｌ．：Ｄｅｔｅｒｍｉｎｉｎｇ  Ｍｏｄｅｌｓ  ｏｆ  Ｓｔｒｕｃｔｕｒｅｓ  ｆｒｏｍ  Ｅａｒｔｈｑｕａｋｅ  Ｒｅｃｏｒｄｓ，１９７８，Ｃ．Ｉ．Ｔ．，Ｐａｓａｄｅｎａ
【発明の概要】
【発明が解決しようとする課題】
【０００８】  従来の解析手法のベック手法は，対象系が運動方程式５に従うことを仮定して，最小二乗法に基づいて系の観測記録より次数毎に固有値および固有ベクトルを推定する方法である．例えば，第ｓ次の固有値λｓと固有ベクトルｕｓを対象とした場合，式６を観測方程式とする．
【数５】
2014021956.tif 000011
ここに，Ｍ，ＣおよびＫは，それぞれ質量，減衰，剛性の行列であり，ｘ（ｔ），ｆ（ｔ）は相対変位，外力を意味するベクトルである．入力が加速度であれば外力ｆ（ｔ）は慣性力となる．また，・は時間ｔに関する微分を意味する．
【数６】
2014021956.tif 000012
ここに，λｒ，ｕｒは第ｒ次固有値，固有ベクトルであり，αｒ（ｔ）は相対変位ｘ（ｔ）をｕｒの線形
2014021956.tif 000013
る．そして，Ｒｅａｌは実数部を，Σは第ｓ次モードを除いた他のモードの和を意味する．【０００９】  式６は，右辺の第ｓ次モード定数λｓ，ｕｓと同様に変化しつつあるパラメータλｒ，ｕｒを真値とみなして設定した観測方程式である．最小二乗法に従って誤差を評価する場合は，回帰式の値を母集合の中の標本に対応させて行うが，式６は，標本に当たる左辺が未知変数を含むため標本自身が種々に変動していき，解の収束性に疑問が持たれる．従って，一自由度系を除いて式６の右辺を回帰式として使用するのは適当ではない．そして，相対加
2014021956.tif 000014
るが，最小二乗法に従って解いているのは固有ベクトルのみである．ただし，固有ベクトルｕｓは式６に基づいて計算をしているため，上記と同様の問題が存在する．−方，固有値λｓは，試行的に現在値の周辺に微小変化させ，その中で最小の誤差を与える値を最確値
2014021956.tif 000015
固有ベクトルｕｒの線形展開式で表現した運動方程式５の解を回帰式とし，式７を観測方程式とする．
【数７】
2014021956.tif 000016
上式より，本発明に対して誤差分布を正規分布と仮定した最小二乗法の本来のアルゴリズムが適用される．即ち，全次数の固有値λｒおよび固有ベクトルｕｒをパラメータに設定して，非線形パラメータλｒと線形パラメータｕｒを個別に計算することにより，安定した収束演算が期待される，また，固有値λｒの初期値の推定は容易であり，通常，刺激係数の大きいモードから順に算出される．【００１０】  以上，本発明は，対象系が運動方程式に従うことを仮定して最小二乗法に基づいて観測方程式を定式化し，演算の安定性および収束性を確保して，系の動的特性を示す固有値と固有ベクトルを精度良く算出する時間領域のモード解析手法，時間領域のモード解析プログラムおよび時間領域のモード解析プログラムを記録したコンピュータ読み取り可能な記録媒体を提供することを課題とする．
【課題を解決するための手段】
【００１１】  自由度Ｎの振動系の運動方程式は式８のように表される．対象系の運動は，この方程式に従うことを仮定する．
【数８】
2014021956.tif 000017
ここに，Ｍ，ＣおよびＫは，それぞれ質量，減衰，剛性の（Ｎ×Ｎ）行列であり，ｘ（ｔ），ｆ（ｔ）は相対変位，外力を意味するＮ元ベクトルである．入力が加速度であれば，外力ｆ（ｔ）は慣性力となる．以下では観測記録として加速度記録を扱う．これら入出力関係は式８によって定義される．【００１２】  時間領域のモード解析は，対象系の動的特性を示す固有値と固有ベクトルの２つのモード定数を観測記録より逆解析的に算出することを目的とする．そこで，はじめに観測記録の真値となる式８の解について，モード定数の陽関数に表現するために対象系の固有値問題を取り上げ，系の基本ベクトルを構成する固有ベクトルによる線形展開式について考えていく．【００１３】  ここで，以下の展開の便宜のために，式８を次のように改める．
【数９】
2014021956.tif 000018
ここに，Ａ，Ｂ，Ｘ（ｔ），Ｆ（ｔ）は，次の通りである．
【数１０】
2014021956.tif 000019
ここに，Ａ，Ｂは（２Ｎ×２Ｎ）の行列であり，Ｘ（ｔ），Ｆ（ｔ）は（２Ｎ×１）のベクトルである．【００１４】  式９の固有値問題は，Ｆ（ｔ）＝０と置いた閉鎖系において論じられる．このとき解Ｘ（ｔ）をＸ（ｔ）＝Ｖｅｘｐ（λｔ）と仮定すれば，式９は次の固有方程式へ移行する．
【数１１】
2014021956.tif 000020
ここに，λｒ，Ｖｒはそれぞれ第ｒ次の固有値および固有ベクトルであり，以下の内容を持つ．
【数１２】
2014021956.tif 000021
ここに，ｆｒ，ｈｒはそれぞれ第ｒ次固有周波数および固有減衰定数であり，ｉは虚数単位である．
【数１３】
2014021956.tif 000022
ここに，ベクトルｕｒは第ｒ次変位ベクトルである．【００１５】  また，ベクトルＶｒは次の直交関係式を満たす．
【数１４】
2014021956.tif 000023
【数１５】
2014021956.tif 000024
ただし，上の両式は，ｒ＝ｓのとき２λｓｕｓＴＭｕｓ＋ｕｓＴＣｕｓ，−λｓ２ｕｓＴＭｕｓ＋ｕｓＴＫｕｓとなるが，これらの値が，それぞれ１．０，λｓとなるように規格化されたＶｓについて成立する．以下ではこのように規格化されたＶｓあるいはｕｓを使用する．【００１６】  式１１を解けば閉鎖系の運動が明らかになる．そして，式９のべクトルＸ（ｔ）は，基本ベクトルを構成する固有ベクトルＶｒによって式１６のように線形展開式に表される．
【数１６】
2014021956.tif 000025
ここに，αｒ（ｔ）は展開係数である．また，Σはモードの総和を意味し，以後明示しない場合はｒ＝１〜２Ｎであるとする．２Ｎ個のモード定数の内のＮ個は，他のＮ個の共役複素数である．【００１７】
2014021956.tif 000026
【数１７】
2014021956.tif 000027
ここに，展開係数αｒ（ｔ）は，式１８の微分方程式を満たす．
【数１８】
2014021956.tif 000028
上式は，式１６を式９に代入し，式１４，１５の直交性より得られる．よって，上式を解けば相対加速度が求まり，式１９のようになる．
【数１９】
2014021956.tif 000029
ここに，αｒｏはαｒ（ｔ）の初期値であり，ηｒ（ｔ）＝ｅｘｐ（λｒｔ）である．【００１８】  上式の｛｝内のベクトルβｒ（ｔ）は式２０の微分方程式を満たす．
【数２０】
2014021956.tif 000030
尚，式１９は，モード定数λｒ，ｕｒの共役複素数を考慮すれば次のようになる．
【数２１】
2014021956.tif 000031
ここに，Ｒｅａｌは実数部を意味する．【００１９】  時間領域のモード解析の２つのパラメータ固有値λｒ，固有ベクトルｕｒは，式２１と観測
2014021956.tif 000032
有値は非線形の関係に，固有ベクトルは線形の関係にあるため，これらを個別に扱う．【００２０】  先ず，固有ベクトルｕｒは，固有値λｒを既知として扱い，観測点毎に計算する．式２１よ
2014021956.tif 000033
【数２２】
2014021956.tif 000034
式２２を固有ベクトル計算時の回帰式とする．未知変数αｒｏ，ｕｒｊは，上式の理論値と観測値間の差の二乗和を評価して求められる．このとき誤差評価関数は，全ての観測データが等精度の場合，式２３のようになる．
【数２３】
2014021956.tif 000035
2014021956.tif 000036
味する．記号Σは，解析対象時間内の総和である．【００２１】  αｒｏ，ｕｒｊを含む新たな未知変数ａｒｊ，Ｕｒｊは，式２４のεｊの最小化必要条件式より正規方程式を立て，これを解いて求められる．
【数２４】
2014021956.tif 000037
一方，固有値λｒは，固有ベクトルｕｒを既知として扱う．前述したように，固有値λｒは
2014021956.tif 000038
についてテーラー展開をして，その一次まで考慮し線形化を行う．即ち，観測点ｊの時間
2014021956.tif 000039
【数２５】
2014021956.tif 000040
2014021956.tif 000041
【００２２】  式２５が固有値計算時の回帰式である．未知変数Δλｒは，式２５の値と観測値間の差の二乗和を評価して求められる．このとき，誤差評価関数εは式２６のようになる．
【数２６】
2014021956.tif 000042
Δλｒは式２７のεの最小化必要条件式より正規方程式を立て，これを解いて求められる．そして，固有値λｒはΔλｒが加えられ更新される．刺激係数の最も大きいモードは，加速度応答の大部分を占めるので，式２７ではじめに対象となるモードは最大刺激係数を持つモードである．
【数２７】
2014021956.tif 000043
以上のように，固有値および固有ベクトルは個別に求め，式２６の誤差が所定の許容値に達するまで繰り返し行われる．【００２３】  時間領域のモード解析より求められるモード定数を使用して相対速度，相対変位が式２８，２９よりそれぞれ計算される．
【数２８】
2014021956.tif 000044
【数２９】
2014021956.tif 000045
また，対象系の伝達関数は，以下のように作成される．【００２４】  式９のベクトルＸ（ｔ），外力Ｆ（ｔ）のフーリエスペクトルをそれぞれＹ（ω），Ｇ（ω）とすれば，式９の両辺をフーリエ変換して次式を得る．
【数３０】
2014021956.tif 000046
ここに，ωは角振動数であり，ｉは虚数単位である．【００２５】  ここで式３１のように，上式のベクトルＹ（ω）の式１１の固有ベクトルＶｒによる線形展開式を考える．
【数３１】
2014021956.tif 000047
上式の展開係数γｒ（ω）は，上式を式３０へ代入して直交性式１４，１５の両式を考慮すれば求められ，Ｙ（ω），そして式１０，１３より，相対変位ｘ（ｔ）のフーリエスペクトルｙ（ω）は次のようになる．
【数３２】
2014021956.tif 000048
【数３３】
2014021956.tif 000049
ここに，ｇ（ω）は，外力ｆ（ｔ）のフーリエスペクトルである．ただし，前述したようにモード解析では固有ベクトルｕｒに換えてＵｒｊ＝ｕｒｊｕｒ（ｊは観測点を指す）を計算しており，実際には，式３４により，観測点ｊの伝達関数ｙｊ（ω）が求められる．
【数３４】
2014021956.tif 000050
伝達関数は，加速度入力の場合と加振力入力で異なり，それぞれ次のようになる．【００２６】  地震波のような加速度入力の場合は，伝達関数は基準点の絶対値に対する対象点の絶対値（加速度，速度あるいは変位の各絶対値）の比で作成される．従って，基準点の絶対変位ｘｂ（ｔ）のフーリエスペクトルをｙｂ（ω丿とすれば，伝達関数Ｔｊ（ω）は式３５のようになる．
【数３５】
2014021956.tif 000051
2014021956.tif 000052
【数３６】
2014021956.tif 000053
ここに，Ｉは単位ベクトルである．これより式３５の伝達関数Ｔｊ（ω）は次のように表される．
【数３７】
2014021956.tif 000054
ところで，刺激係数ｕｒＴＭＩは系の定数であり，Ｍ単独では任意の値を採り得る．つまり，本解析から算出されるＵｒｊは適当なＭに応じた値となるが，刺激関数ＵｒｊＴＭＩ（＝ｕｒＴＭＩｕｒｊ）はＭに依らず一定である．【００２７】  一方，単位力加振の場合，対象系ｊ点の伝達関数Ｔｊ（ω）は，式３８のようになる．
【数３８】
2014021956.tif 000055
ここに，ｅｋは観測点ｋに単位力が働いた場合，第ｋ要素に１．０，他の要素に０．０を持つベクトルである．加振点が複数の場合は，対応する要素に１．０を与えればよい．
【発明の効果】
【００２８】  対象系の運動が運動方程式に従うことを仮定して最小二乗法に従って定式化した本発明によれば，系本来の固有値および固有ベクトルを精度良く算出可能な時間領域のモード解析手法，時間領域のモード解析プログラムおよび時間領域のモード解析プログラムを記録したコンピュータ読み取り可能な記録媒体を提供することができる．【００２９】  本発明によれば，加速度応答に対してそれぞれ非線形，線形の関係にある固有値および固有ベクトルを個別に計算して安定した高い収束性を確保し，精度良く算出可能な時間領域のモード解析手法，時間領域のモード解析プログラムおよび時間領域のモード解析プログラムを記録したコンピュータ読み取り可能な記録媒体を提供することができる．【００３０】  本発明によれば，算出された固有値および固有ベクトルを使用して，対象系の変位，速度および周波数特性を示す伝達関数を精度良く計算可能な時間領域のモード解析手法，時間領域のモード解析プログラムおよび時間領域のモード解析プログラムを記録したコンピュータ読み取り可能な記録媒体を提供することができる．【００３１】  時間領域で解析を行う本発明によれば，対象系にとって重要な，刺激係数のより大きいモードの固有値と固有ベクトルから順に精度良く計算可能な時間領域のモード解析手法，時間領域のモード解析プログラムおよび時間領域のモード解析プログラムを記録したコンピュータ読み取り可能な記録媒体を提供することができる．
【図面の簡単な説明】
【００３２】
【図１】図１は，本発明の実施形態に係る時間領域のモード解析手法を説明する流れ図である．
【図２】図２は，本発明の実施形態に係る時間領域のモード解析プログラム（不図示），時間領域のモード解析プログラムを記録したコンピュータ読み取り可能な記録媒体および時間領域のモード解析プログラムをメモリーに読み込み解析を実行するコンピュータの外観図である．
【発明を実施するための形態】
【００３３】  本発明の実施形態について添付の図面を参照しながら説明する．具体的には，（１）時間領域のモード解析手法，（２）時間領域のモード解析プログラムおよび時間領域のモード解析プログラムを記録したコンピュータ読み取り可能な記録媒体について説明する．
【００３４】（１）時間領域のモード解析手法の説明
【００３５】  図１は本発明の実施形態に係る時間領域のモード解析手法を説明する流れ図である．はじめに，予めコンピュータのハードディスクやコンパクトディスクＣＤ等の記録媒体に記録されている観測記録を読み込む（ステップＳ１）．次に，次数ｒに１を与え（ステップＳ２），刺激係数の最も大きいモードの固有値および固有ベクトルを解析対象とする．そして，非線形パラメータ固有値λｒの初期値を読み込み（ステップＳ３），その固有値に属する線形パラメータ固有ベクトルａｒｊ，Ｕｒｊの計算を行う（ステップＳ４）．求められたこれらモード定数を使用して加速度の理論値を求め（ステップＳ５），観測記録と理論値間の誤差の二乗和ε１を計算する（ステップＳ６）．
【００３６】  繰り返し演算に入り，繰り返し回数Ｌｏｏｐを１にセットする（ステップＳ７）．最初に固有ベクトルａｒｊ，Ｕｒｊを固定して，固有値の微小修正量Δλｒを計算する（ステップＳ８）．次に，その微小修正量より固有値を更新して（ステップＳ９），この固有値に属する固有ベクトル
2014021956.tif 000056
二乗和ε２を求め（ステップＳ１２），ε２のε１に対する相対誤差｜ε１−ε２｜／ε１を計算して誤差の許容値Ｅｒｒｏｒと比較する（ステップＳ１３）．誤差が収束して相対誤差がＥｒｒｏｒ以下であると判定された場合（Ｙｅｓ）は，次数ｒを１上げる（ステップＳ１４）．その次数が所定のモード数Ｍｏｄｅ以下と判定された場合（ステップＳ１５でＹｅｓ）はステップＳ３へ戻り，それまでの収束した各次の固有値と新たに加えた次数に対する固有値の初期値を読み込む（ステップＳ３）．次数が変わればモード定数にわずかながら変化が見られるので，前の次数で収束した固有値および固有ベクトルも改めてパラメータとして処理する．モード次数ｒがＭｏｄｅより大きい場合（ステップＳ１５でＮｏ）は，求められた全次数の固有値と固有ベクトルを使用して対象系の相対速度，相対変位（ステップＳ１６），ならびに伝達関数を計算し（ステップＳ１７），解析結果を図示等の後処理のためにハードディスクあるいはＣＤ等の記録媒体に出力して（ステップＳ１８）解析を終了する．ステップＳ１３で相対誤差がＥｒｒｏｒより大きいと判定された場合（Ｎｏ）は，繰り返し回数Ｌｏｏｐを１増やし（ステップＳ１９），それが予め設定した反復回数ＮＬｏｏｐと比較して（ステップＳ２０）小さい場合（Ｙｅｓ）は，相対誤差の基準値ε１をε２に換えて（ステップＳ２１）ステップＳ８へ進み，収束演算を繰り返す．大きい場合（Ｎｏ）は，一旦演算をストップして（ステップＳ２２），固有値の初期値や観測記録を見直して再検討する．
【００３７】（２）時間領域のモード解析プログラムおよび時間領域のモード解析プログラムを記録したコンピュータ読み取り可能な記録媒体の説明
【００３８】  図２は本発明の実施形態に係る時間領域のモード解析プログラムおよび時間領域のモード解析プログラムを記録したコンピュータ読み取り可能な記録媒体を説明する図である．コンピュータ３１０は，一般に処理部のＣＰＵ，メモリーやハードディスクのような記憶部，入力部，出力部および表示部等で構成されている．また，コンパクトディスクＣＤ３００のような記録媒体に対して書き込みおよび読み出しを可能にするドライバが備えられている．
【００３９】  図１を用いて説明した時間領域のモード解析は処理部ＣＰＵにおいて実行される．このとき，時間領域のモード解析プログラム（不図示）は，予め記録されているＣＤ３００等の記録媒体からメモリーに読み込まれて実行される．あるいは，予めハードディスクに記憶されている同プログラムをメモリーに読み込み実行されることもある．また，ネットワークを介して同プログラムを取得して実行されることも考えられる．
【符号の説明】  ３００−−−コンパクトディスクＣＤ等の記録媒体，３１０−−−コンピュータ
【図１】
2014021956.tif 000122
【図２】
2014021956.tif 000123