(57)【特許請求の範囲】
【請求項１】
時刻ｔn のときの加速度ａn を用いて，時刻ｔn から時間τだけ経過した時刻ｔs ＝ｔn ＋τのときの加速度ａs を近似することにより，外力ｆが与えられたときの次式（１）で表される多質点線形振動系の時刻歴応答をコンピュータにより演算する振動解析方法において，
上記加速度ａs を，次式（２）の通り，上記時間τの２次式で近似してなることを特徴とする振動解析方法。
Ｍａ＋Ｃｖ＋Ｋｄ＝ｆ                                            （１）
ａs ＝ａn ＋ｂτ＋ｃτ2                                         （２）
但し，Ｍ，Ｃ，Ｋはそれぞれ質量，減衰，剛性マトリクス，ａは加速度ベクトル，ｖは速度ベクトル，ｄは変位ベクトル，ｆは外力ベクトルである。また，状態量ｂ，ｃはそれぞれ加加速度，及びその微分量を表すベクトルである。
【請求項２】
上記時刻ｔs のときの速度ｖs ，変位ｄs を，時刻ｔn のときの速度ｖn ，変位ｄn を用いて次式（３），（４）によりそれぞれ近似してなる請求項１に記載の振動解析方法。
ｖs ＝ｖn ＋ａn τ＋０．５ｂτ2 ＋ｅ3 ｃτ3                     （３）
ｄs ＝ｄn ＋ｖn τ＋０．５ａn τ2 ＋ｅ6 ｂτ3 ＋ｅ12ｃτ4       （４）
但し，ｅ3 ，ｅ6 ，ｅ12は係数である。
【請求項３】
時刻ｔn から演算周期ｈだけ経過した時刻ｔn+1 ＝ｔn ＋ｈにおける運動方程式（５），及び上記時間τが演算周期ｈのθ（θ＞１）倍の時間であるときの上記時刻ｔs ＝ｔn ＋θｈにおける運動方程式（６）を満たすように，上記状態量ｂ，ｃを演算してなる請求項２に記載の振動解析方法。
Ｍａn+1 ＋Ｃｖn+1 ＋Ｋｄn+1 ＝ｆn+1                             （５）
Ｍａs   ＋Ｃｖs   ＋Ｋｄs   ＝ｆs                               （６）
【請求項４】
上記状態量ｂ，ｃを次式（７），及び（７ａ）乃至（７ｆ）に従って演算してなる請求項３に記載の振動解析方法。
【数１】
0003640843.tif 000012
【請求項５】
外力ｆが与えられたときの次式（１）で表される多質点線形振動系の時刻歴応答を，時刻ｔn のときの加速度ａn を用いた次式（２）によって，時刻ｔn から時間τだけ経過した時刻ｔs ＝ｔn ＋τのときの加速度ａs を近似することにより，コンピュータにより演算する振動解析方法であって，
初期変位ｄ0 ，初期速度ｖ0 を含む初期条件，及び次式（８）に基づき初期加速度ａ0 を算出する初期加速度算出工程と，
上記初期条件，上記初期加速度算出工程により算出された初期加速度，及び次式（７ａ）〜（７ｆ）に従って，行列要素Ａ11，Ａ12，Ａ21，Ａ22，ｐ1 ，ｐ2 を算出するマトリクス算出工程と，
上記マトリクス算出工程により算出された上記行列要素Ａ11，Ａ12，Ａ21，Ａ22，ｐ1 ，ｐ2 ，及び次式（７）に従って，加加速度，及び加加速度の微分量をそれぞれ表す状態量ｂ，ｃを算出する状態量算出工程と，
上記状態量算出工程により算出された上記状態量ｂ，ｃ，及び次式（２）〜（４）に従って，時刻ｔn から演算周期ｈだけ時間が経過した時刻ｔn+1 における加速度ａn+1 ，速度ｖn+1 ，変位ｄn+1 を含む次時刻状態量を算出する次時刻状態量算出工程とを具備し，
上記次時刻状態量算出工程により算出された上記次時刻状態量を用いて，上記行列要素Ａ11，Ａ12，Ａ21，Ａ22，ｐ1 ，ｐ2 を更新し，上記状態量算出工程，上記次時刻状態量算出工程を順次行う工程を繰り返すことにより，上記時刻歴応答を求めてなる振動解析方法。
Ｍａ＋Ｃｖ＋Ｋｄ＝ｆ                                            （１）
ａs ＝ａn ＋ｂτ＋ｃτ2                                         （２）
ｖs ＝ｖn ＋ａn τ＋０．５ｂτ2 ＋ｅ3 ｃτ3                     （３）
ｄs ＝ｄn ＋ｖn τ＋０．５ａn τ2 ＋ｅ6 ｂτ3 ＋ｅ12ｃτ4       （４）
Ｍａn+1 ＋Ｃｖn+1 ＋Ｋｄn+1 ＝ｆn+1                             （５）
Ｍａs   ＋Ｃｖs   ＋Ｋｄs   ＝ｆs                               （６）
【数２】
0003640843.tif 000013
ａ0 ＝Ｍ（−Ｃｖ0 −Ｋｄ0 ＋ｆ0 ）                              （８）
但し，Ｍ，Ｃ，Ｋはそれぞれ質量，減衰，剛性マトリクス，ａは加速度ベクトル，ｖは速度ベクトル，ｄは変位ベクトル，ｆは外力ベクトル，ａn+1 ，ｖn+1 ，ｄn+1 ，ｆn+1 はそれぞれ時刻ｔn から演算周期ｈだけ経過した時刻ｔn+1 における加速度ベクトル，速度ベクトル，変位ベクトル，外力ベクトル，ａs ，ｖs ，ｄs ，ｆs はそれぞれ時刻ｔn から時間τ（＝θｈ（θ＞１））だけ経過した時刻ｔs における加速度ベクトル，速度ベクトル，変位ベクトル，外力ベクトルである。また，状態量ｂ，ｃはそれぞれ加加速度，及びその微分量を表すベクトルであり，ｅ3 ，ｅ6 ，ｅ12は係数である。
【請求項６】
θ＝１．２５，ｅ3 ＝１／３，ｅ6 ＝１／６，ｅ12＝１／１２である請求項３〜５のいずれか１項に記載の振動解析方法。
【請求項７】
時刻ｔn のときの加速度ａn を用いて，時刻ｔn から時間τだけ経過した時刻ｔs ＝ｔn ＋τのときの加速度ａs を近似することにより，外力ｆが与えられたときの次式（１）で表される多質点線形振動系の時刻歴応答を演算する振動解析装置において，
上記加速度ａs を，次式（２）の通り，上記時間τの２次式で近似してなることを特徴とする振動解析装置。
Ｍａ＋Ｃｖ＋Ｋｄ＝ｆ                                            （１）
ａs ＝ａn ＋ｂτ＋ｃτ2                                         （２）
但し，Ｍ，Ｃ，Ｋはそれぞれ質量，減衰，剛性マトリクス，ａは加速度ベクトル，ｖは速度ベクトル，ｄは変位ベクトル，ｆは外力ベクトルである。また，状態量ｂ，ｃはそれぞれ加加速度，及びその微分量を表すベクトルである。
【請求項８】
上記時刻ｔs のときの速度ｖs ，変位ｄs を，時刻ｔn のときの速度ｖn ，変位ｄn を用いて次式（３），（４）によりそれぞれ近似してなる請求項７に記載の振動解析装置。
ｖs ＝ｖn ＋ａn τ＋０．５ｂτ2 ＋ｅ3 ｃτ3                     （３）
ｄs ＝ｄn ＋ｖn τ＋０．５ａn τ2 ＋ｅ6 ｂτ3 ＋ｅ12ｃτ4       （４）
但し，ｅ3 ，ｅ6 ，ｅ12は係数である。
【請求項９】
時刻ｔn から演算周期ｈだけ経過した時刻ｔn+1 ＝ｔn ＋ｈにおける運動方程式（５），及び上記時間τが演算周期ｈのθ（θ＞１）倍の時間であるときの上記時刻ｔs ＝ｔn ＋θｈにおける運動方程式（６）を満たすように，上記状態量ｂ，ｃを演算してなる請求項８に記載の振動解析装置。
Ｍａn+1 ＋Ｃｖn+1 ＋Ｋｄn+1 ＝ｆn+1                             （５）
Ｍａs   ＋Ｃｖs   ＋Ｋｄs   ＝ｆs                               （６）
【請求項１０】
上記状態量ｂ，ｃを次式（７），及び（７ａ）乃至（７ｆ）に従って演算してなる請求項９に記載の振動解析装置。
【数３】
0003640843.tif 000014
【請求項１１】
外力ｆが与えられたときの次式（１）で表される多質点線形振動系の時刻歴応答を，時刻ｔn のときの加速度ａn を用いた次式（２）によって，時刻ｔn から時間τだけ経過した時刻ｔs ＝ｔn ＋τのときの加速度ａs を近似することにより，演算する振動解析装置であって，
初期変位ｄ0 ，初期速度ｖ0 を含む初期条件，及び次式（８）に基づき初期加速度ａ0 を算出する初期加速度算出手段と，
上記初期条件，上記初期加速度算出手段により算出された初期加速度，及び次式（７ａ）〜（７ｆ）に従って，行列要素Ａ11，Ａ12，Ａ21，Ａ22，ｐ1 ，ｐ2 を算出するマトリクス算出手段と，
上記マトリクス算出手段により算出された上記行列要素Ａ11，Ａ12，Ａ21，Ａ22，ｐ1 ，ｐ2 ，及び次式（７）に従って，加加速度，及び加加速度の微分量をそれぞれ表す状態量ｂ，ｃを算出する状態量算出手段と，
上記状態量算出手段により算出された上記状態量ｂ，ｃ，及び次式（２）〜（４）に従って，時刻ｔn から演算周期ｈだけ時間が経過した時刻ｔn+1 における加速度ａn+1 ，速度ｖn+1 ，変位ｄn+1 を含む次時刻状態量を算出する次時刻状態量算出手段とを具備し，
上記次時刻状態量算出手段により算出された上記次時刻状態量を用いて，上記行列要素Ａ11，Ａ12，Ａ21，Ａ22，ｐ1 ，ｐ2 を更新し，上記状態量算出手段による上記状態量の演算，上記次時刻状態量算出手段による上記次時刻状態量の演算を繰り返し，上記時刻歴応答を求めてなる振動解析装置。
Ｍａ＋Ｃｖ＋Ｋｄ＝ｆ                                            （１）
ａs ＝ａn ＋ｂτ＋ｃτ2                                         （２）
ｖs ＝ｖn ＋ａn τ＋０．５ｂτ2 ＋ｅ3 ｃτ3                     （３）
ｄs ＝ｄn ＋ｖn τ＋０．５ａn τ2 ＋ｅ6 ｂτ3 ＋ｅ12ｃτ4       （４）
Ｍａn+1 ＋Ｃｖn+1 ＋Ｋｄn+1 ＝ｆn+1                             （５）
Ｍａs   ＋Ｃｖs   ＋Ｋｄs   ＝ｆs                               （６）
【数４】
0003640843.tif 000015
ａ0 ＝Ｍ（−Ｃｖ0 −Ｋｄ0 ＋ｆ0 ）                              （８）
但し，Ｍ，Ｃ，Ｋはそれぞれ質量，減衰，剛性マトリクス，ａは加速度ベクトル，ｖは速度ベクトル，ｄは変位ベクトル，ｆは外力ベクトル，ａn+1 ，ｖn+1 ，ｄn+1 ，ｆn+1 はそれぞれ時刻ｔn から演算周期ｈだけ経過した時刻ｔn+1 における加速度ベクトル，速度ベクトル，変位ベクトル，外力ベクトル，ａs ，ｖs ，ｄs ，ｆs はそれぞれ時刻ｔn から時間τ（＝θｈ（θ＞１））だけ経過した時刻ｔs における加速度ベクトル，速度ベクトル，変位ベクトル，外力ベクトルである。また，状態量ｂ，ｃはそれぞれ加加速度，及びその微分量を表すベクトルであり，ｅ3 ，ｅ6 ，ｅ12は係数である。
【請求項１２】
θ＝１．２５，ｅ3 ＝１／３，ｅ6 ＝１／６，ｅ12＝１／１２である請求項９〜１１のいずれか１項に記載の振動解析装置。
【請求項１３】
時刻ｔn のときの加速度ａn を用いて，時刻ｔn から時間τだけ経過した時刻ｔs ＝ｔn ＋τのときの加速度ａs を近似することにより，外力ｆが与えられたときの次式（１）で表される多質点線形振動系の時刻歴応答をコンピュータに演算させる振動解析プログラムを記録したコンピュータ読み取り可能な記録媒体において，
上記加速度ａs を，次式（２）の通り，上記時間τの２次式で近似してなることを特徴とする振動解析プログラムを記録したコンピュータ読み取り可能な記録媒体。
Ｍａ＋Ｃｖ＋Ｋｄ＝ｆ                                            （１）
ａs ＝ａn ＋ｂτ＋ｃτ2                                         （２）
但し，Ｍ，Ｃ，Ｋはそれぞれ質量，減衰，剛性マトリクス，ａは加速度ベクトル，ｖは速度ベクトル，ｄは変位ベクトル，ｆは外力ベクトルである。また，状態量ｂ，ｃはそれぞれ加加速度，及びその微分量を表すベクトルである。
【請求項１４】
上記時刻ｔs のときの速度ｖs ，変位ｄs を，時刻ｔn のときの速度ｖn ，変位ｄn を用いて次式（３），（４）によりそれぞれ近似してなる請求項１３に記載の振動解析プログラムを記録したコンピュータ読み取り可能な記録媒体。
ｖs ＝ｖn ＋ａn τ＋０．５ｂτ2 ＋ｅ3 ｃτ3                     （３）
ｄs ＝ｄn ＋ｖn τ＋０．５ａn τ2 ＋ｅ6 ｂτ3 ＋ｅ12ｃτ4       （４）
但し，ｅ3 ，ｅ6 ，ｅ12は係数である。
【請求項１５】
時刻ｔn から演算周期ｈだけ経過した時刻ｔn+1 ＝ｔn ＋ｈにおける運動方程式（５），及び上記時間τが演算周期ｈのθ（θ＞１）倍の時間であるときの上記時刻ｔs ＝ｔn ＋θｈにおける運動方程式（６）を満たすように，上記状態量ｂ，ｃを演算してなる請求項１４に記載の振動解析プログラムを記録したコンピュータ読み取り可能な記録媒体。
Ｍａn+1 ＋Ｃｖn+1 ＋Ｋｄn+1 ＝ｆn+1                             （５）
Ｍａs   ＋Ｃｖs   ＋Ｋｄs   ＝ｆs                               （６）
【請求項１６】
上記状態量ｂ，ｃを次式（７），及び（７ａ）乃至（７ｆ）に従って演算してなる請求項１５に記載の振動解析プログラムを記録したコンピュータ読み取り可能な記録媒体。
【数５】
0003640843.tif 000016
【請求項１７】
外力ｆが与えられたときの次式（１）で表される多質点線形振動系の時刻歴応答を，時刻ｔn のときの加速度ａn を用いた次式（２）によって，時刻ｔn から時間τだけ経過した時刻ｔs ＝ｔn ＋τのときの加速度ａs を近似することにより，コンピュータに演算させる振動解析プログラムを記録したコンピュータ読み取り可能な記録媒体であって，
コンピュータを，初期変位ｄ0 ，初期速度ｖ0 を含む初期条件，及び次式（８）に基づき初期加速度ａ0 を算出する初期加速度算出手段，
上記初期条件，上記初期加速度算出手段により算出された初期加速度，及び次式（７ａ）〜（７ｆ）に従って，行列要素Ａ11，Ａ12，Ａ21，Ａ22，ｐ1 ，ｐ2 を算出するマトリクス算出手段，
上記マトリクス算出手段により算出された上記行列要素Ａ11，Ａ12，Ａ21，Ａ22，ｐ1 ，ｐ2 ，及び次式（７）に従って，加加速度，及び加加速度の微分量をそれぞれ表す状態量ｂ，ｃを算出する状態量算出手段，
上記状態量算出手段により算出された上記状態量ｂ，ｃ，及び次式（２）〜（４）に従って，時刻ｔn から演算周期ｈだけ時間が経過した時刻ｔn+1 における加速度ａn+1 ，速度ｖn+1 ，変位ｄn+1 を含む次時刻状態量を算出する次時刻状態量算出手段として機能させ，
コンピュータに，上記次時刻状態量算出手段により算出された上記次時刻状態量を用いて，上記行列要素Ａ11，Ａ12，Ａ21，Ａ22，ｐ1 ，ｐ2 を更新させ，上記状態量算出手段による上記状態量の演算，上記次時刻状態量算出手段による上記次時刻状態量の演算を繰り返させ，上記時刻歴応答を求めさせる振動解析プログラムを記録したコンピュータ読み取り可能な記録媒体。
Ｍａ＋Ｃｖ＋Ｋｄ＝ｆ                                            （１）
ａs ＝ａn ＋ｂτ＋ｃτ2                                         （２）
ｖs ＝ｖn ＋ａn τ＋０．５ｂτ2 ＋ｅ3 ｃτ3                     （３）
ｄs ＝ｄn ＋ｖn τ＋０．５ａn τ2 ＋ｅ6 ｂτ3 ＋ｅ12ｃτ4       （４）
Ｍａn+1 ＋Ｃｖn+1 ＋Ｋｄn+1 ＝ｆn+1                             （５）
Ｍａs   ＋Ｃｖs   ＋Ｋｄs   ＝ｆs                               （６）
【数６】
0003640843.tif 000017
ａ0 ＝Ｍ（−Ｃｖ0 −Ｋｄ0 ＋ｆ0 ）                              （８）
但し，Ｍ，Ｃ，Ｋはそれぞれ質量，減衰，剛性マトリクス，ａは加速度ベクトル，ｖは速度ベクトル，ｄは変位ベクトル，ｆは外力ベクトル，ａn+1 ，ｖn+1 ，ｄn+1 ，ｆn+1 はそれぞれ時刻ｔn から演算周期ｈだけ経過した時刻ｔn+1 における加速度ベクトル，速度ベクトル，変位ベクトル，外力ベクトル，ａs ，ｖs ，ｄs ，ｆs はそれぞれ時刻ｔn から時間τ（＝θｈ（θ＞１））だけ経過した時刻ｔs における加速度ベクトル，速度ベクトル，変位ベクトル，外力ベクトルである。また，状態量ｂ，ｃはそれぞれ加加速度，及びその微分量を表すベクトルであり，ｅ3 ，ｅ6 ，ｅ12は係数である。
【請求項１８】
θ＝１．２５，ｅ3 ＝１／３，ｅ6 ＝１／６，ｅ12＝１／１２である請求項１５〜１７のいずれか１項に記載の振動解析プログラムを記録したコンピュータ読み取り可能な記録媒体。
【発明の詳細な説明】
【０００１】
【発明の属する技術分野】
本発明は，振動解析方法，振動解析装置，及び振動解析プログラムを記録したコンピュータ読み取り可能な記録媒体に係り，詳しくは，時刻ｔn のときの加速度ａn を用いて，時刻ｔn から時間τだけ経過した時刻ｔs ＝ｔn ＋τのときの加速度ａs を近似することにより，外力が与えられたときの多質点線形振動系の時刻歴応答を演算する振動解析方法，振動解析装置，及び振動解析プログラムを記録したコンピュータ読み取り可能な記録媒体に関するものである。
【０００２】
【従来の技術】
多質点線形振動系の運動方程式は次のＭＣＫ形の方程式（１）で表される。
Ｍａ＋Ｃｖ＋Ｋｄ＝ｆ                                            （１）
ここで，Ｍ，Ｃ，Ｋはそれぞれ質量，減衰，剛性マトリクス，ａは加速度ベクトル，ｖは速度ベクトル，ｄは変位ベクトル，ｆは外力ベクトルである。
外力ｆが与えられたときの時刻歴応答は，式（１）を時間積分することにより求められる。（数値）積分法としては，種々の方法があるが，構造振動系の問題には，一般的にNewmark 法が用いられることが多い。
上記Newmark 法は，加速度を一次式で取り扱う線形加速度法の一種であり，時刻ｔから微小時間Δｔが経過した時刻ｔ＋Δｔにおける変位ベクトルｄ（ｔ＋Δｔ），速度ベクトルｖ（ｔ＋Δｔ）は，時刻ｔにおける変位ベクトルｄ（ｔ），速度ベクトルｖ（ｔ）に対して，それぞれ次式（９），（１０）により表される。
0003640843.tif 000002
上記Newmark 法は，線形加速度法における係数（１／３！）をパラメータβに置き換えたものであり，β＝１／６とする他，β＝１／４として用いられる場合が多い。
【０００３】
【発明が解決しようとする課題】
例えばβ＝１／４としたとき上記Newmark 法は，無条件安定で減衰精度も比較的良いが，図１０に示す如く，周期精度は優れず，また硬い系では高次モードの不要な振動が現れてしまう。尚，図１０におけるｈは演算周期であり，Ｔは振動周期であり，図１０（ａ）は振動に対して演算が比較的細かい場合，図１０（ｂ）は振動に対して演算が比較的粗い場合を示す。上記Newmark 法の他，数値積分法としては，高次モードの不要な振動を抑える数値減衰機能を有したWilson法やHoubout 法なども用いられるが，これらは精度に問題がある。
上記構造振動系の問題に適した解法としては，少なくとも安定性，精度，数値減衰機能の３つの条件を満足することが望まれるが，従来の数値積分法では，調整し得るパラメータも少なく，これら３つの条件を全て満足させることができなかった。
このため，従来の数値積分法を用いて構造振動系の時刻歴応答を定める振動解析装置や振動解析方法により計算を行う場合，また，コンピュータ読み取り可能なＣＤ−ＲＯＭ等の記録媒体に記録された振動解析プログラムを用いてコンピュータに計算を行わせる場合にも，上記３つの条件のうちのいずれか又は複数に問題が生じていた。
本発明は，このような従来の技術における課題を解決するために，振動解析方法，振動解析装置，及び振動解析プログラムを改良し，安定性を確保しながら，必要とする周波数帯域では極めて精度がよく，しかも不必要な高周波モードを減衰させることができる振動解析方法，振動解析装置，及び振動解析プログラムを記録したコンピュータ読み取り可能な記録媒体を提供することを目的とするものである。
【０００４】
【課題を解決するための手段】
上記目的を達成するために，請求項１に係る発明は，時刻ｔn のときの加速度ａn を用いて，時刻ｔn から時間τだけ経過した時刻ｔs ＝ｔn ＋τのときの加速度ａs を近似することにより，外力ｆが与えられたときの上式（１）で表される多質点線形振動系の時刻歴応答をコンピュータにより演算する振動解析方法において，上記加速度ａs を，次式（２）の通り，上記時間τの２次式で近似してなることを特徴とする振動解析方法として構成されている。
ａs ＝ａn ＋ｂτ＋ｃτ2                                         （２）
但し，Ｍ，Ｃ，Ｋはそれぞれ質量，減衰，剛性マトリクス，ａは加速度ベクトル，ｖは速度ベクトル，ｄは変位ベクトル，ｆは外力ベクトルである。また，状態量ｂ，ｃはそれぞれ加加速度，及びその微分量を表すベクトルである。
また，請求項２に係る発明は，上記請求項１に記載の振動解析方法において，上記時刻ｔs のときの速度ｖs ，変位ｄs を，時刻ｔn のときの速度ｖn ，変位ｄn を用いて次式（３），（４）によりそれぞれ近似してなることをその要旨とする。
ｖs ＝ｖn ＋ａn τ＋０．５ｂτ2 ＋ｅ3 ｃτ3                     （３）
ｄs ＝ｄn ＋ｖn τ＋０．５ａn τ2 ＋ｅ6 ｂτ3 ＋ｅ12ｃτ4       （４）
但し，ｅ3 ，ｅ6 ，ｅ12は係数である。
また，請求項３に係る発明は，上記請求項２に記載の振動解析方法において，時刻ｔn から演算周期ｈだけ経過した時刻ｔn+1 ＝ｔn ＋ｈにおける運動方程式（５），及び上記時間τが演算周期ｈのθ（θ＞１）倍の時間であるときの上記時刻ｔs ＝ｔn ＋θｈにおける運動方程式（６）を満たすように，上記状態量ｂ，ｃを演算してなることをその要旨とする。
Ｍａn+1 ＋Ｃｖn+1 ＋Ｋｄn+1 ＝ｆn+1                             （５）
Ｍａs   ＋Ｃｖs   ＋Ｋｄs   ＝ｆs                               （６）
【０００５】
また，請求項４に係る発明は，上記請求項３に記載の振動解析方法において，上記状態量ｂ，ｃを次式（７），及び（７ａ）乃至（７ｆ）に従って演算してなることをその要旨とする。
【数７】
0003640843.tif 000003
また，請求項５に係る発明は，外力ｆが与えられたときの上式（１）で表される多質点線形振動系の時刻歴応答を，時刻ｔn のときの加速度ａn を用いた上式（２）によって，時刻ｔn から時間τだけ経過した時刻ｔs ＝ｔn ＋τのときの加速度ａs を近似することにより，コンピュータにより演算する振動解析方法であって，初期変位ｄ0 ，初期速度ｖ0 を含む初期条件，及び次式（８）に基づき初期加速度ａ0 を算出する初期加速度算出工程と，上記初期条件，上記初期加速度算出工程により算出された初期加速度，及び上式（７ａ）〜（７ｆ）に従って，行列要素Ａ11，Ａ12，Ａ21，Ａ22，ｐ1 ，ｐ2 を算出するマトリクス算出工程と，上記マトリクス算出工程により算出された上記行列要素Ａ11，Ａ12，Ａ21，Ａ22，ｐ1 ，ｐ2 ，及び上式（７）に従って，加加速度，及び加加速度の微分量をそれぞれ表す状態量ｂ，ｃを算出する状態量算出工程と，上記状態量算出工程により算出された上記状態量ｂ，ｃ，及び上式（２）〜（４）に従って，時刻ｔn から演算周期ｈだけ時間が経過した時刻ｔn+1 における加速度ａn+1 ，速度ｖn+1 ，変位ｄn+1 を含む次時刻状態量を算出する次時刻状態量算出工程とを具備し，上記次時刻状態量算出工程により算出された上記次時刻状態量を用いて，上記行列要素Ａ11，Ａ12，Ａ21，Ａ22，ｐ1 ，ｐ2 を更新し，上記状態量算出工程，上記次時刻状態量算出工程を順次行う工程を繰り返すことにより，上記時刻歴応答を求めてなる振動解析方法である。
ａ0 ＝Ｍ（−Ｃｖ0 −Ｋｄ0 ＋ｆ0 ）                              （８）
また，請求項６に係る発明は，上記請求項３〜５のいずれか１項に記載の振動解析方法において，θ＝１．２５，ｅ3 ＝１／３，ｅ6 ＝１／６，ｅ12＝１／１２であることをその要旨とする。
【０００６】
上記請求項１〜４のいずれか１項に記載の発明によれば，外力ｆが与えられたときの上式（１）で表される多質点線形振動系の時刻歴応答を演算する際に，時刻ｔn における加速度ａn を用いて，時刻ｔn から時間τだけ経過した時刻ｔs における加速度ａs が，上式（２）の通り，時間τの２次式にて近似されるため，調整し得るパラメータ（θ，ｅ3 ，ｅ6 ，ｅ12）を増やすことができ，その結果安定性を確保しながら，必要とする周波数帯域では極めて精度がよく，しかも不必要な高周波モードを減衰させることができる振動解析方法を提供することができる。
調整し得るパラメータが増加するのと同時に未知数も増加するから，それを同定するための方程式が必要になるが，例えば時刻ｔn+1 ＝ｔn ＋ｈ（演算周期）における上記運動方程式（５），及び時刻ｔs ＝ｔn ＋θｈ（θ＞１）における上記運動方程式（６）を用いればよい。
また，上記請求項５に記載の発明によれば，外力ｆが与えられたときの上式（１）で表される多質点線形振動系の時刻歴応答を演算する際に，時刻ｔn における加速度ａn を用いて，時刻ｔn から時間τだけ経過した時刻ｔs における加速度ａs が，上式（２）の通り，時間τの２次式にて近似され，時刻ｔn における速度ｖn を用いて，上式（３）の通り，時刻ｔs における速度ｖs が近似され，時刻ｔn における変位ｄn を用いて，上式（４）の通り，時刻ｔs における変位ｄs が近似され，これらの式（２）〜（４），及び上記運動方程式（５），（６）から構成される連立方程式（７），（７ａ）〜（７ｆ）から状態量ｂ，ｃが導き出され，次の演算周期の加速度ａn+1 ，ｖn+1 ，ｄn+1 が順次求められるため，調整し得るパラメータ（θ，ｅ3 ，ｅ6 ，ｅ12）を増やすことができ，その結果安定性を確保しながら，必要とする周波数帯域では極めて精度がよく，しかも不必要な高周波モードを減衰させることができる振動解析方法を提供することができる。
また，上記請求項１〜５のいずれか１項に記載の振動解析方法では，上記パラメータθ，ｅ3 ，ｅ6 ，ｅ12は，安定性，精度，数値減衰機能の観点から，θ＝１．２５，ｅ3 ＝１／３，ｅ6 ＝１／６，ｅ12＝１／１２とするのが好ましい。
【０００７】
また，請求項７に係る発明は，時刻ｔn のときの加速度ａn を用いて，時刻ｔn から時間τだけ経過した時刻ｔs ＝ｔn ＋τのときの加速度ａs を近似することにより，外力ｆが与えられたときの上式（１）で表される多質点線形振動系の時刻歴応答を演算する振動解析装置において，上記加速度ａs を，上式（２）の通り，上記時間τの２次式で近似してなることを特徴とする振動解析装置である。
また，請求項８に係る発明は，上記請求項７に記載の振動解析装置において，上記時刻ｔs のときの速度ｖs ，変位ｄs を，時刻ｔn のときの速度ｖn ，変位ｄn を用いて上式（３），（４）によりそれぞれ近似してなることをその要旨とする。
また，請求項９に係る発明は，上記請求項８に記載の振動解析装置において，時刻ｔn から演算周期ｈだけ経過した時刻ｔn+1 ＝ｔn ＋ｈにおける上記運動方程式（５），及び上記時間τが演算周期ｈのθ（θ＞１）倍の時間であるときの上記時刻ｔs ＝ｔn ＋θｈにおける上記運動方程式（６）を満たすように，上記状態量ｂ，ｃを演算してなることをその要旨とする。
また，請求項１０に係る発明は，上記請求項９に記載の振動解析装置において，上記状態量ｂ，ｃを上式（７），及び（７ａ）乃至（７ｆ）に従って演算してなることをその要旨とする。
また，請求項１１に係る発明は，外力ｆが与えられたときの上式（１）で表される多質点線形振動系の時刻歴応答を，時刻ｔn のときの加速度ａn を用いた上式（２）によって，時刻ｔn から時間τだけ経過した時刻ｔs ＝ｔn ＋τのときの加速度ａs を近似することにより，演算する振動解析装置であって，初期変位ｄ0 ，初期速度ｖ0 を含む初期条件，及び上式（８）に基づき初期加速度ａ0 を算出する初期加速度算出手段と，上記初期条件，上記初期加速度算出手段により算出された初期加速度，及び上式（７ａ）〜（７ｆ）に従って，行列要素Ａ11，Ａ12，Ａ21，Ａ22，ｐ1 ，ｐ2 を算出するマトリクス算出手段と，上記マトリクス算出手段により算出された上記行列要素Ａ11，Ａ12，Ａ21，Ａ22，ｐ1 ，ｐ2 ，及び上式（７）に従って，加加速度，及び加加速度の微分量をそれぞれ表す状態量ｂ，ｃを算出する状態量算出手段と，上記状態量算出手段により算出された上記状態量ｂ，ｃ，及び上式（２）〜（４）に従って，時刻ｔn から演算周期ｈだけ時間が経過した時刻ｔn+1 における加速度ａn+1 ，速度ｖn+1 ，変位ｄn+1 を含む次時刻状態量を算出する次時刻状態量算出手段とを具備し，上記次時刻状態量算出手段により算出された上記次時刻状態量を用いて，上記行列要素Ａ11，Ａ12，Ａ21，Ａ22，ｐ1 ，ｐ2 を更新し，上記状態量算出手段による上記状態量の演算，上記次時刻状態量算出手段による上記次時刻状態量の演算を繰り返し，上記時刻歴応答を求めてなる振動解析装置である。
また，請求項１２に係る発明は，上記請求項９〜１１のいずれか１項に記載の振動解析装置において，θ＝１．２５，ｅ3 ＝１／３，ｅ6 ＝１／６，ｅ12＝１／１２であることをその要旨とする。
【０００８】
上記請求項７〜１０のいずれか１項に記載の発明によれば，外力ｆが与えられたときの上式（１）で表される多質点線形振動系の時刻歴応答を演算する際に，時刻ｔn における加速度ａn を用いて，時刻ｔn から時間τだけ経過した時刻ｔs における加速度ａs が，上式（２）の通り，時間τの２次式にて近似されるため，調整し得るパラメータ（θ，ｅ3 ，ｅ6 ，ｅ12）を増やすことができ，その結果安定性を確保しながら，必要とする周波数帯域では極めて精度がよく，しかも不必要な高周波モードを減衰させることができる振動解析装置を提供することができる。
調整し得るパラメータが増加するのと同時に未知数も増加するから，それを同定するための方程式が必要になるが，例えば時刻ｔn+1 ＝ｔn ＋ｈ（演算周期）における上記運動方程式（５），及び時刻ｔs ＝ｔn ＋θｈ（θ＞１）における上記運動方程式（６）を用いればよい。
また，上記請求項１１に記載の発明によれば，外力ｆが与えられたときの上式（１）で表される多質点線形振動系の時刻歴応答を演算する際に，時刻ｔn における加速度ａn を用いて，時刻ｔn から時間τだけ経過した時刻ｔs における加速度ａs が，上式（２）の通り，時間τの２次式にて近似され，時刻ｔn における速度ｖn を用いて，上式（３）の通り，時刻ｔs における速度ｖs が近似され，時刻ｔn における変位ｄn を用いて，上式（４）の通り，時刻ｔs における変位ｄs が近似され，これらの式（２）〜（４），及び上記運動方程式（５），（６）から構成される連立方程式（７），（７ａ）〜（７ｆ）から状態量ｂ，ｃが導き出され，次の演算周期の加速度ａn+1 ，ｖn+1 ，ｄn+1 が順次求められるため，調整し得るパラメータ（θ，ｅ3 ，ｅ6 ，ｅ12）を増やすことができ，その結果安定性を確保しながら，必要とする周波数帯域では極めて精度がよく，しかも不必要な高周波モードを減衰させることができる振動解析装置を提供することができる。
また，上記請求項７〜１１のいずれか１項に記載の振動解析装置では，上記パラメータθ，ｅ3 ，ｅ6 ，ｅ12は，安定性，精度，数値減衰機能の観点から，θ＝１．２５，ｅ3 ＝１／３，ｅ6 ＝１／６，ｅ12＝１／１２とするのが好ましい。
【０００９】
また，請求項１３に係る発明は，時刻ｔn のときの加速度ａn を用いて，時刻ｔn から時間τだけ経過した時刻ｔs ＝ｔn ＋τのときの加速度ａs を近似することにより，外力ｆが与えられたときの上式（１）で表される多質点線形振動系の時刻歴応答をコンピュータに演算させる振動解析プログラムを記録したコンピュータ読み取り可能な記録媒体において，上記加速度ａs を，上式（２）の通り，上記時間τの２次式で近似してなることを特徴とする振動解析プログラムを記録したコンピュータ読み取り可能な記録媒体である。
また，請求項１４に係る発明は，上記請求項１３に記載の振動解析プログラムを記録したコンピュータ読み取り可能な記録媒体において，上記時刻ｔs のときの速度ｖs ，変位ｄs を，時刻ｔn のときの速度ｖn ，変位ｄn を用いて上式（３），（４）によりそれぞれ近似してなることをその要旨とする。
また，請求項１５に係る発明は，上記請求項１４に記載の振動解析プログラムを記録したコンピュータ読み取り可能な記録媒体において，時刻ｔn から演算周期ｈだけ経過した時刻ｔn+1 ＝ｔn ＋ｈにおける上記運動方程式（５），及び上記時間τが演算周期ｈのθ（θ＞１）倍の時間であるときの上記時刻ｔs ＝ｔn ＋θｈにおける上記運動方程式（６）を満たすように，上記状態量ｂ，ｃを演算してなることをその要旨とする。
また，請求項１６に係る発明は，上記請求項１５に記載の振動解析プログラムを記録したコンピュータ読み取り可能な記録媒体において，上記状態量ｂ，ｃを上式（７），及び（７ａ）乃至（７ｆ）に従って演算してなることをその要旨とする。
また，請求項１７に係る発明は，外力ｆが与えられたときの次式（１）で表される多質点線形振動系の時刻歴応答を，時刻ｔn のときの加速度ａn を用いた次式（２）によって，時刻ｔn から時間τだけ経過した時刻ｔs ＝ｔn ＋τのときの加速度ａs を近似することにより，コンピュータに演算させる振動解析プログラムを記録したコンピュータ読み取り可能な記録媒体であって，コンピュータを，初期変位ｄ0 ，初期速度ｖ0 を含む初期条件，及び次式（８）に基づき初期加速度ａ0 を算出する初期加速度算出手段，上記初期条件，上記初期加速度算出手段により算出された初期加速度，及び上式（７ａ）〜（７ｆ）に従って，行列要素Ａ11，Ａ12，Ａ21，Ａ22，ｐ1 ，ｐ2 を算出するマトリクス算出手段，上記マトリクス算出手段により算出された上記行列要素Ａ11，Ａ12，Ａ21，Ａ22，ｐ1 ，ｐ2 ，及び上式（７）に従って，加加速度，及び加加速度の微分量をそれぞれ表す状態量ｂ，ｃを算出する状態量算出手段，上記状態量算出手段により算出された上記状態量ｂ，ｃ，及び次式（２）〜（４）に従って，時刻ｔn から演算周期ｈだけ時間が経過した時刻ｔn+1 における加速度ａn+1 ，速度ｖn+1 ，変位ｄn+1 を含む次時刻状態量を算出する次時刻状態量算出手段として機能させ，コンピュータに，上記次時刻状態量算出手段により算出された上記次時刻状態量を用いて，上記行列要素Ａ11，Ａ12，Ａ21，Ａ22，ｐ1 ，ｐ2 を更新させ，上記状態量算出手段による上記状態量の演算，上記次時刻状態量算出手段による上記次時刻状態量の演算を繰り返させ，上記時刻歴応答を求めさせる振動解析プログラムを記録したコンピュータ読み取り可能な記録媒体である。
また，請求項１８に係る発明は，上記請求項１５〜１７のいずれか１項に記載の振動解析プログラムを記録したコンピュータ読み取り可能な記録媒体において，θ＝１．２５，ｅ3 ＝１／３，ｅ6 ＝１／６，ｅ12＝１／１２であることをその要旨とする。
【００１０】
上記請求項１３〜１６のいずれか１項に記載の発明によれば，外力ｆが与えられたときの上式（１）で表される多質点線形振動系の時刻歴応答を演算する際に，時刻ｔn における加速度ａn を用いて，時刻ｔn から時間τだけ経過した時刻ｔs における加速度ａs が，上式（２）の通り，時間τの２次式にて近似されるため，調整し得るパラメータ（θ，ｅ3 ，ｅ6 ，ｅ12）を増やすことができ，その結果安定性を確保しながら，必要とする周波数帯域では極めて精度がよく，しかも不必要な高周波モードを減衰させることができる振動解析プログラムを記録したコンピュータ読み取り可能な記録媒体を提供することができる。
調整し得るパラメータが増加するのと同時に未知数も増加するから，それを同定するための方程式が必要になるが，例えば時刻ｔn+1 ＝ｔn ＋ｈ（演算周期）における上記運動方程式（５），及び時刻ｔs ＝ｔn ＋θｈ（θ＞１）における上記運動方程式（６）を用いればよい。
また，上記請求項１７に記載の発明によれば，外力ｆが与えられたときの上式（１）で表される多質点線形振動系の時刻歴応答を演算する際に，時刻ｔn における加速度ａn を用いて，時刻ｔn から時間τだけ経過した時刻ｔs における加速度ａs が，上式（２）の通り，時間τの２次式にて近似され，時刻ｔn における速度ｖn を用いて，上式（３）の通り，時刻ｔs における速度ｖs が近似され，時刻ｔn における変位ｄn を用いて，上式（４）の通り，時刻ｔs における変位ｄs が近似され，これらの式（２）〜（４），及び上記運動方程式（５），（６）から構成される連立方程式（７），（７ａ）〜（７ｆ）から状態量ｂ，ｃが導き出され，次の演算周期の加速度ａn+1 ，ｖn+1 ，ｄn+1 が順次求められるため，調整し得るパラメータ（θ，ｅ3 ，ｅ6 ，ｅ12）を増やすことができ，その結果安定性を確保しながら，必要とする周波数帯域では極めて精度がよく，しかも不必要な高周波モードを減衰させることができる振動解析プログラムを記録したコンピュータ読み取り可能な記録媒体を提供することができる。
また，上記請求項１３〜１７のいずれか１項に記載の振動解析プログラムを記録したコンピュータ読み取り可能な記録媒体では，上記パラメータθ，ｅ3 ，ｅ6 ，ｅ12は，安定性，精度，数値減衰機能の観点から，θ＝１．２５，ｅ3 ＝１／３，ｅ6 ＝１／６，ｅ12＝１／１２とするのが好ましい。
【００１１】
【発明の実施の形態】
以下，添付図面を参照して，本発明の実施の形態につき説明し，本発明の理解に供する。尚，以下の実施の形態は，本発明の具体的な一例であって，本発明の技術的範囲を限定する性格のものではない。ここに，図１は本発明の実施の形態に係る振動解析方法を説明するためのフローチャート，図２は本発明の実施の形態に係る振動解析装置の概略構成を示す機能ブロック図，図３は上記振動解析装置として機能させるコンピュータの概略構成を示す機能ブロック図である。
まず，本発明の実施の形態に係る振動解析方法は，図１に示す如く，外力ｆが与えられたときの上式（１）で表される多質点線形振動系の時刻歴応答を，時刻ｔn のときの加速度ａn を用いた上式（２）によって，時刻ｔn から時間τだけ経過した時刻ｔs ＝ｔn ＋τのときの加速度ａs を近似することにより，演算する振動解析方法であって，上式（１）の運動方程式を導出する工程Ｓ１と，パラメータθ，ｅ6 ，ｅ12，演算周期ｈ，初期速度ｖ0 ，初期変位ｄ0 を含む初期条件を設定する工程Ｓ２と，上記初期条件ｖ0 ，ｄ0 ，及び上式（８）に基づき初期加速度ａ0 を算出する初期加速度算出工程Ｓ３と，上記初期条件，上記初期加速度算出工程Ｓ３により算出された初期加速度，及び上式（７ａ）〜（７ｆ）に従って，行列要素Ａ11，Ａ12，Ａ21，Ａ22，ｐ1 ，ｐ2 を算出するマトリクス算出工程Ｓ４と，上記マトリクス算出工程Ｓ４により算出された上記行列要素Ａ11，Ａ12，Ａ21，Ａ22，ｐ1 ，ｐ2 ，及び上式（７）に従って，加加速度，及び加加速度の微分量をそれぞれ表す状態量ｂ，ｃを算出する状態量算出工程Ｓ５と，上記状態量算出工程Ｓ５により算出された上記状態量ｂ，ｃ，及び上式（２）〜（４）に従って，時刻ｔn から演算周期ｈだけ時間が経過した時刻ｔn+1 における加速度ａn+1 ，速度ｖn+1 ，変位ｄn+1 を含む次時刻状態量を算出する次時刻状態量算出工程Ｓ６とを具備し，上記次時刻状態量算出工程Ｓ６により算出された上記次時刻状態量を用いて，上記行列要素Ａ11，Ａ12，Ａ21，Ａ22，ｐ1 ，ｐ2 を更新し，上記状態量算出工程Ｓ５，上記次時刻状態量算出工程Ｓ６を順次行う工程を繰り返すことにより，上記時刻歴応答を求めるものである。
【００１２】
また，本発明の実施の形態に係る振動解析方法を実施するのに好適な（本発明の実施の形態に係る）振動解析装置は，図２に示す如く，外力ｆが与えられたときの上式（１）で表される多質点線形振動系の時刻歴応答を，時刻ｔn のときの加速度ａn を用いた上式（２）によって，時刻ｔn から時間τだけ経過した時刻ｔs ＝ｔn ＋τのときの加速度ａs を近似することにより，演算する振動解析装置であって，初期変位ｄ0 ，初期速度ｖ0 を含む初期条件，及び上式（８）に基づき初期加速度ａ0 を算出する初期加速度算出手段１と，上記初期条件，上記初期加速度算出手段１により算出された初期加速度，及び上式（７ａ）〜（７ｆ）に従って，行列要素Ａ11，Ａ12，Ａ21，Ａ22，ｐ1 ，ｐ2 を算出するマトリクス算出手段２と，上記マトリクス算出手段２により算出された上記行列要素Ａ11，Ａ12，Ａ21，Ａ22，ｐ1 ，ｐ2 ，及び上式（７）に従って，加加速度，及び加加速度の微分量をそれぞれ表す状態量ｂ，ｃを算出する状態量算出手段３と，上記状態量算出手段３により算出された上記状態量ｂ，ｃ，及び上式（２）〜（４）に従って，時刻ｔn から演算周期ｈだけ時間が経過した時刻ｔn+1 における加速度ａn+1 ，速度ｖn+1 ，変位ｄn+1 を含む次時刻状態量を算出する次時刻状態量算出手段４とを具備し，上記次時刻状態量算出手段４により算出された上記次時刻状態量を用いて，上記行列要素Ａ11，Ａ12，Ａ21，Ａ22，ｐ1 ，ｐ2 を更新し，上記状態量算出手段３による上記状態量の演算，上記次時刻状態量算出手段４による上記次時刻状態量の演算を繰り返し，上記時刻歴応答を求めるものである。
上記振動解析装置は，例えば上記各手段をメモリやハードディスクなどの記憶装置上に生成し演算装置に上記時刻歴応答を演算させる振動解析プログラムにより機能させられたコンピュータとして具体化される。
上記コンピュータ３０には，図３に示す如く，マウス，キーボードなどの入力装置３１，ディスプレイなどの表示装置３２，演算装置３３，メモリ３４，ハードディスクなどの記憶装置３５，プログラム読み取り用のＣＤ−ＲＯＭドライブなどの記録媒体読取装置３６等の構成を備えた一般的なものを用いることができる。
上記振動解析プログラムは，ＣＤ−ＲＯＭ等のコンピュータ読み取り可能な記録媒体３７に無圧縮，又は圧縮状態で記録されるものであり，必要に応じてコンピュータ３０の記憶装置３５に実行可能な状態で転送記憶される。上記コンピュータ３０は，入力装置３１を介して使用者から上記振動解析プログラムの実行の指示があると，上記演算装置３３により上記記録媒体読取装置３６又は上記記憶装置３５から上記振動解析プログラムを読み出し，上記振動解析プログラムの指令に従って，上記メモリ３４や上記記憶装置３５上に上記振動解析装置の各手段を生成させる。
【００１３】
次に，本発明に係る解析理論，上記振動解析装置，及び上記振動解析プログラムの詳細について説明する。
時刻ｔn における加速度，速度，変位をａn ，ｖn ，ｄn と表し，時刻ｔn から時間τ秒間経過した時刻ｔs ＝ｔn ＋τにおける加速度ａs ，速度ｖs ，変位ｄs と表すものとすると，本発明に係る解析理論では，上記時刻ｔs における加速度ａs ，速度ｖs ，変位ｄs は，上記時刻ｔn における加速度ａn ，速度ｖn ，変位ｄn を用いて，それぞれ次式（２）〜（４）により表される。
ａs ＝ａn ＋ｂτ＋ｃτ2                                         （２）
ｖs ＝ｖn ＋ａn τ＋０．５ｂτ2 ＋ｅ3 ｃτ3                     （３）
ｄs ＝ｄn ＋ｖn τ＋０．５ａn τ2 ＋ｅ6 ｂτ3 ＋ｅ12ｃτ4       （４）
但し，状態量ｂ，ｃはそれぞれ加加速度，及びその微分量を表すベクトルであり，ｅ3 ，ｅ6 ，ｅ12は係数である。
これらは，加速度ａs が２次式で変化すると仮定し，それを積分して速度ｖs ，変位ｄs を表したものである。単に，上式（２）の加速度ａs を積分しただけでは，上式（３）及び（４）のｅ3 ，ｅ6 ，ｅ12は，それぞれ１／３，１／６，１／１２となるが，安定性や精度を調整するために，上記係数ｅ3 ，ｅ6 ，ｅ12をパラメータとして用いる。
従来のNewmark 法等では，加速度ａs は１次式で変化すると仮定しているため，本発明による数値積分の方が精度が向上する。しかしながら，加速度ａs が２次式で変化すると仮定することにより，求めるべき未知数が増加する。このために，上記時間τとして演算周期ｈ，及びθｈ（θ＞１）の２つの時刻を考え，それぞれの時刻に対して，上式（１）の運動方程式を用いる。即ち，次式（１ａ），（１ｂ）を用いる。
Ｍａn+1 ＋Ｃｖn+1 ＋Ｋｄn+1 ＝ｆn+1                           （１ａ）
Ｍａs   ＋Ｃｖs   ＋Ｋｄs   ＝ｆs                             （１ｂ）
まず，時刻ｔn+1 ＝ｔn ＋ｈにおける，釣り合いを考える。上式（２）〜（４）のτをｔn+1 に置き換えて，上式（１ａ）に代入すると，次式（１１）が得られる。
0003640843.tif 000004
次に，時刻ｔs ＝ｔn ＋θｈにおける，釣り合いを考える。上式（２）〜（４）を上式（１ｂ）に代入すると，次式（１２）が得られる。
0003640843.tif 000005
上式（１１），（１２）をまとめ，行列で表示すると，次式（７），（７ａ）〜（７ｆ）のようになる。
【数８】
0003640843.tif 000006
時刻ｔn における状態量ａn ，ｖn ，ｄn が既知であれば，上式（７）の右辺が算定でき，この１次連立方程式を解くことにより，上記状態量ｂ，ｃを求めることができる。上記状態量ｂ，ｃが求められたら，時刻ｔn から演算周期ｈが経過した次のステップ時刻ｔn+1 における加速度ａn+1 ，速度ｖn+1 ，変位ｄn+1 を次式（２ａ），（３ａ），（４ａ）に従って求めることができる。
ａn+1 ＝ａn ＋ｂｈ＋ｃｈ2                                     （２ａ）
ｖn+1 ＝ｖn ＋ａn ｈ＋０．５ｂｈ2 ＋ｅ3 ｃｈ3                 （３ａ）
ｄn+1 ＝ｄn ＋ｖn ｈ＋０．５ａn ｈ2 ＋ｅ6 ｂｈ3 ＋ｅ12ｃｈ4   （４ａ）
そして，次時刻ｔn+1 における状態量ａn+1 ，速度ｖn+1 ，変位ｄn+1 が求まったら，これを現時刻ｔn として，上記の演算を繰り返し，順次次時刻ｔn+1 における上記状態量を求めれば，上式（１）の多質点振動系の時刻歴応答を得ることができる。
尚，最初のステップでは，初期条件（初期速度ｖ0 ，初期変位ｄ0 ）に基づき，初期加速度ａ0 を求めておく必要がある。
また，必要に応じて上式（４ａ）の係数ｅ6 を係数ｅ6 ’に代えた次式（４ａ’）を次のステップ時刻ｔn+1 における変位ｄn+1 の演算に用い，上式（４ａ）を出力用の式として用いるようにしてもよい。
0003640843.tif 000007
【００１４】
上記の解析理論に従って演算を行うために，上記振動解析装置は，コンピュータ３０のメモリ３４又は記憶装置３５上に生成される上記した各手段を具備する。上記振動解析装置において，上記初期加速度算出手段１は，使用者により入力装置３１を介して指示があると，演算装置３３を用いて，上記初期条件ｖ0 ，ｄ0 ，及び上式（１）を変形した次式（８）に従って，初期加速度ａ0 を求めるものである。
ａ0 ＝Ｍ（−Ｃｖ0 −Ｋｄ0 ＋ｆ0 ）                              （８）
また，上記マトリクス算出手段２は，演算しようとする時刻ｔn+1 の前の時刻ｔn における加速度ａn ，速度ｖ0 ，変位ｄ0 ，また最初のステップでは，上記初期条件ｖ0 ，ｄ0 ，及び上記初期加速度算出手段１により算出された初期加速度ａ0 が与えられると，設定された演算周期ｈやパラメータθ，ｅ6 ，ｅ12等を参照しながら，上式（７ａ）〜（７ｆ）に従って，上記演算装置３３により行列要素Ａ11，Ａ12，Ａ21，Ａ22，ｐ1 ，ｐ2 を算出させるものである。
また，上記状態量算出手段３は，上記マトリクス算出手段２により上記行列要素Ａ11，Ａ12，Ａ21，Ａ22，ｐ1 ，ｐ2 が算出されると，上式（７）に従って，上記演算装置３３により，加加速度，及び加加速度の微分量をそれぞれ表す状態量ｂ，ｃを算出させるものである。
さらに，上記次時刻状態量算出手段４は，上記状態量算出手段３により上記状態量ｂ，ｃが算出されると，上式（２）〜（４）に従って，時刻ｔn から演算周期ｈだけ時間が経過した時刻ｔn+1 における加速度ａn+1 ，速度ｖn+1 ，変位ｄn+1 を含む次時刻状態量を上記演算装置３３を用いて算出するものである。
上記振動解析プログラムは，上記振動解析装置に係る各手段を上記コンピュータ３０上のメモリ３４や記憶装置３５上に生成させる他，表示装置３２にパラメータθ，ｅ6 ，ｅ12などを入力する欄を表示させたり，演算された時刻歴応答を表示させるなど，使用者とのインターフェースを提供する。また，設定されたパラメータや演算された状態量などのデータは，必要に応じてメモリ３４や記憶装置３５に記憶され参照されるが，このデータの管理も行う。
【００１５】
次に，上記振動解析装置の動作の詳細について説明する（上記振動解析方法の説明も兼ねる）。
はじめに，上式（１）の運動方程式が定められる（Ｓ１）。即ち，質量マトリクスＭ，減衰マトリクスＣ，剛性マトリクスＫ，外力ベクトルｆが定められる。次に，演算周期ｈ，パラメータθ，ｅ6 ，ｅ12，初期条件ｖ0 ，ｄ0 の設定が行われる（Ｓ２）。この設定は，例えば，上記振動解析プログラムの指令に従って上記表示装置３２に表示されるパラメータ入力欄に，使用者が，キーボードやマウスなどの入力装置３１を用いて所望の数値の入力を行うことにより行われる。
使用者により入力された演算周期ｈや，パラメータθ，ｅ6 ，ｅ12，初期条件ｖ0 ，ｄ0 などは，上記初期加速度算出手段１に供給され，上式（８）に従って，初期加速度ａ0 が算出される（Ｓ３）。
また，上記演算周期ｈや，パラメータθ，ｅ6 ，ｅ12，初期条件ｖ0 ，ｄ0 などは，上記初期加速度算出手段１により算出された初期加速度ａ0 とともに，上記マトリクス算出手段２に供給され，上式（７ａ）〜（７ｆ）に従って，行列要素Ａ11，Ａ12，Ａ21，Ａ22，ｐ1 ，ｐ2 が算出される（Ｓ４）。
上記行列要素Ａ11，Ａ12，Ａ21，Ａ22，ｐ1 ，ｐ2 は，上記状態量算出手段３に供給され，上式（７）に従って，上記状態量ｂ，ｃが算出される（Ｓ５）。
上記状態量算出手段３により算出された上記状態量ｂ，ｃは，次時刻状態量算出手段４に供給され，上式（２ａ），（３ａ），（４ａ）（又は（４ａ’））に従って，初期状態から演算周期ｈに相当する時間だけが経過した時刻ｔ1 における加速度ａ1 ，速度ｖ1 ，変位ｄ1 が求められる（Ｓ６）。
次に，上記時刻ｔ1 が，予め設定した演算終了時刻ｔend よりも大きいか否かの判別が行われる（Ｓ７）。
上記時刻ｔ1 が，上記演算終了時刻ｔend よりも小さい場合には，演算対象となる時刻が，時刻ｔ1 から演算周期ｈに相当する時間だけが経過した時刻ｔ2 に移行される。また，処理が上記工程Ｓ４に戻され，上記初期条件ｖ0 ，ｄ0 ，及び初期加速度ａ0 の代わりに，先のステップで演算された上記時刻ｔ1 における速度ｖ1 ，変位ｄ1 ，加速度ａ1 を用いて，上記マトリクス算出手段２により，上記行列要素Ａ11，Ａ12，Ａ21，Ａ22，ｐ1 ，ｐ2 が算出される。そして，この更新された行列要素Ａ11，Ａ12，Ａ21，Ａ22，ｐ1 ，ｐ2 を用いて，上記工程Ｓ５，Ｓ６が行われ，上記時刻ｔ2 における加速度ａ2 ，速度ｖ2 ，変位ｄ2 が求められる。
上記時刻ｔ2 についての演算が終了すると，再び上記工程Ｓ７において，上記時刻ｔ2 が上記演算終了時刻ｔend よりも大きいか否かの判別が行われ，上記演算終了時刻ｔend よりも演算対象となる時刻ｔi が大きくなるまで，各時刻ｔ3 ，…，ｔi における加速度ａi ，速度ｖi ，変位ｄi の演算が繰り返される。
そして，上記時刻ｔi が，上記演算終了時刻ｔend よりも大きくなったとき，外力ベクトルｆに対する多質点振動系の時刻歴応答の演算が完了する。
尚，上記時刻歴応答の演算では，上記状態量ｂ，ｃを用いたが，ｚb ＝ｂｈ，ｚC ＝ｃｈ2 とおいて，上式（７），（７ａ）〜（７ｆ）を書き直し，未知数を加速度の次元に揃えておく方が演算上都合がよい。
【００１６】
次に，上記振動解析方法，振動解析装置，及び振動解析プログラムによる時刻歴応答演算に関する性能の検討と最適係数の探索を行う。
上記時刻歴応答の演算に関する，安定性，並びに，周期精度及び減衰精度を検討するにあたっては，例えば文献「清水ほか，機論Ｃ（４６−４０１）（昭５５），２７」の手法を用いることができる。
テスト用の方程式として，次式（１３）で表される１自由度系の自由振動方程式を用いる。
ｘ”＋２ζｘ′＋ω2 ｘ＝０                                    （１３）
但し，ｘは変位（ｄに対応），ｘ”はｘの２階の時間微分（ａに対応），ｘ′はｘの１階の時間微分（ｖに対応），ζは減衰係数である。
上式（１３）に対して，上式（１１），（１２）の関係を用いて数値的に次式（１４）の状態推移マトリクスを導く。
【数９】
0003640843.tif 000008
但し，Ｈは演算周期ｈ，係数ｅ3 ，ｅ6 ，ｅ12の関数であり，式（１１），（１２）において，Ｍ＝１，Ｃ＝２ζω，Ｋ＝ω2 とおくことにより，関数Ｈの内容を記述することができる。数値的に関数Ｈを求めるには，（ｄn ，ｖn ，ａn ）T ＝（１，０，０）T ，（０，１，０）T ，（０，０，１）T として，各初期ベクトルに対する応答計算を３回行い，その結果を関数T の各列成分として求めればよい。
このようにして定めた上式（１４）の関数Ｈの固有値Λi とすると，スペクトル半径ρ，減衰誤差Δζ，振動周期Ｔについての周期誤差ΔＴ／Ｔは，それぞれ次式（１５ａ），（１５ｂ），（１５ｃ）に従って求めることができる。
ρ＝max ｜Λi ｜                                            （１５ａ）
Δζ＝ζ’−ζ                                              （１５ｂ）
ΔＴ／Ｔ＝（Ｔ’−Ｔ）／Ｔ                                  （１５ｃ）
但し，絶対値最大の固有値max Λi の実部，虚部をそれぞれσ’，ω’としたとき，
ζ’＝−σ’／（σ'2＋ω'2）0.5                             （１５ｄ）
Ｔ’＝２π／ω’                                            （１５ｅ）
である。
上記スペクトル半径ρ，減衰誤差Δζ，周期誤差ΔＴ／Ｔを用いて評価を行いながら，上記パラメータθ，ｅ3 ，ｅ6 ，ｅ6 ’，ｅ12の最適な値を次の観点から探索する。
１）無条件安定である（０＜ｈ／Ｔ＜∞においてρ≦１）
２）ｈ／Ｔ＜０．１（１周期１０分割）程度まで精度良好
３）ｈ／Ｔ＞１程度以上では人口減衰が作用する（高次モード抑制）
そこで，上式（１３）において減衰係数ζ＝０とした非減衰自由振動方程式を取り上げ，上記パラメータを幾つか変更して広範囲のｈ／Ｔに対するスペクトル半径と精度を求めた。
そのうちの好ましい３候補の係数と性能を表１，２に示す。
【表１】
0003640843.tif 000009
【表２】
0003640843.tif 000010
但し，表１，２におけるρ∞はｈ／Ｔ≫１におけるスペクトル半径，Δζ，ΔＴ／Ｔはそれぞれ減衰精度，周期精度であり，正値が減衰増，周期伸長を表す。
表１，２に示されるPlan1 ，Plan2 は，減衰精度が良いが，ｅ6 ≠ｅ6 ’の場合，Wilson法などと同じくオーバーシュートの問題が生じるため，Fujinew と称された案が好ましい。Fujinew 案のスペクトル半径と精度を図４に示す。
図４では，減衰係数ζを０，０．２，１．０にした場合の結果が示されており，Fujinew 案は無条件安定であることが理解される。
【００１７】
次に，上記Fujinew 案を適用した振動解析方法，振動解析装置，及び振動解析プログラムに係る３つの具体例（第一の具体例，第二の具体例，第三の具体例）について説明する。
（第一の具体例）
例えば上記時刻歴応答の演算を，１自由度非減衰振動系に適用する。
上式（１３）の方程式においてζ＝０とし，初期条件ｄ0 ＝１，ｖ0 ＝０後の応答を演算周期ｈを変えて解析する。Newmark 法（β＝１／４）による結果，理論解（破線）とともに，解析結果を図５に示す。尚，図５（ａ）はｈ／Ｔ＝０．１の設定で上記時刻歴応答の演算結果と理論解を，図５（ｂ）はｈ／Ｔ＝０．１の設定でNewmark 法による解析結果と理論解を，図５（ｃ）はｈ／Ｔ＝０．５の設定でｈ／Ｔ＝０．１の設定で上記時刻歴応答の演算結果と理論解を，図５（ｄ）はNewmark 法による解析結果と理論解をそれぞれ示す図である。
上記図５から理解される通り，上記時刻歴応答の演算では，１周期１０分割程度まではきわめて精度がよく，さらにそれより粗くなると強力な人口減衰が作用して高次モードを抑制することができる。
（第二の具体例）
また，上記時刻歴応答の演算を，図６に示すような２自由度振動系に適用すると，ｍ1 ＝ｍ2 ＝１，ｋ1 ＝ｋ2 ＝１０3 ，ｋ12＝１０6 とした硬い振動系（固有振動数は５．０３２２５Ｈｚであり，ｍ，ｋはそれぞれ質量マトリクスＭ，剛性マトリクスＫの行列要素に対応）において，初期条件ｘ1 ′＝１０，ｘ1 ＝ｘ2 ′＝ｘ2 ＝０を与えたときの，演算周期ｈ＝０．０２秒での過渡応答は，図７に示す通りである。ここで，図７（ａ），（ｂ）は，それぞれ同条件におけるNewmark 法によるｍ1 の加速度，変位であり，図７（ｃ），（ｄ）は，それぞれ本案（Fujinew 案）による加速度，変位である。特に加速度において，上記Newmark 法で生じるハンチングが効果的に抑えられている。
（第三の具体例）
また，上記時刻歴応答の演算を建築構造物にエルセントロ(El Centro) 地震波が入力された場合に適用する。ここで，上記建築構造物が，図８に示すような，集中質量とせん断ばね要素で表現される６質点系モデルにモデル化可能であるとする。固有振動数は低い方から，３．１５，６，９４，９．０７，１３．３，３２．９，１９５Ｈｚである。Ｍ，Ｋの作成の後，Ｃ＝０．００２Ｋとし，エルセントロ地震波３４２Gaに対する応答を計算した計算結果を図９に示す。尚，図９（ａ），（ｂ），（ｃ）は，それぞれ質点ｍ1 ，ｍ3 ，ｍ5 の加速度であり，図９（ｄ）は入力されたエルセントロ地震波である。参考のため，上記第三の具体例について，本案（Fujinew 案）の他，Newmark 法，Wilson法（θ＝１．４２），Collocation 法（θ＝１．２１５８，β＝０．１９）を用いて解析した最大加速度応答を表３に示す。
【表３】
0003640843.tif 000011
このように，本実施の形態に係る振動解析方法，振動解析装置，及び振動解析プログラムを記録したコンピュータ読み取り可能な記録媒体によれば，外力ｆが与えられたときの上式（１）で表される多質点線形振動系の時刻歴応答を演算する際に，時刻ｔn における加速度ａn を用いて，時刻ｔn から時間τだけ経過した時刻ｔs における加速度ａs が，上式（２）の通り，時間τの２次式にて近似されるため，調整し得るパラメータ（θ，ｅ3 ，ｅ6 ，ｅ12）を増やすことができ，その結果安定性を確保しながら，必要とする周波数帯域では極めて精度がよく，しかも不必要な高周波モードを減衰させることができる振動解析プログラムを記録したコンピュータ読み取り可能な記録媒体を提供することができる。
調整し得るパラメータが増加するのと同時に未知数も増加するから，それを同定するための方程式が必要になるが，例えば時刻ｔn+1 ＝ｔn ＋ｈ（演算周期）における上記運動方程式（５），及び時刻ｔs ＝ｔn ＋θｈ（θ＞１）における上記運動方程式（６）を用いればよい。
さらに，上記パラメータθ，ｅ3 ，ｅ6 ，ｅ12は，θ＝１．２５，ｅ3 ＝１／３，ｅ6 ＝１／６，ｅ12＝１／１２とすることによって，好ましい安定性，精度，数値減衰機能を得ることができる。
【００１８】
【発明の効果】
以上説明した通り，上記請求項１〜４のいずれか１項に記載の発明によれば，外力ｆが与えられたときの上式（１）で表される多質点線形振動系の時刻歴応答を演算する際に，時刻ｔn における加速度ａn を用いて，時刻ｔn から時間τだけ経過した時刻ｔs における加速度ａs が，上式（２）の通り，時間τの２次式にて近似されるため，調整し得るパラメータ（θ，ｅ3 ，ｅ6 ，ｅ12）を増やすことができ，その結果安定性を確保しながら，必要とする周波数帯域では極めて精度がよく，しかも不必要な高周波モードを減衰させることができる振動解析方法を提供することができる。
また，上記請求項５に記載の発明によれば，外力ｆが与えられたときの上式（１）で表される多質点線形振動系の時刻歴応答を演算する際に，時刻ｔn における加速度ａn を用いて，時刻ｔn から時間τだけ経過した時刻ｔs における加速度ａs が，上式（２）の通り，時間τの２次式にて近似され，時刻ｔn における速度ｖn を用いて，上式（３）の通り，時刻ｔs における速度ｖs が近似され，時刻ｔn における変位ｄn を用いて，上式（４）の通り，時刻ｔs における変位ｄs が近似され，これらの式（２）〜（４），及び上記運動方程式（５），（６）から構成される連立方程式（７），（７ａ）〜（７ｆ）から状態量ｂ，ｃが導き出され，次の演算周期の加速度ａn+1 ，ｖn+1 ，ｄn+1 が順次求められるため，調整し得るパラメータ（θ，ｅ3 ，ｅ6 ，ｅ12）を増やすことができ，その結果安定性を確保しながら，必要とする周波数帯域では極めて精度がよく，しかも不必要な高周波モードを減衰させることができる振動解析方法を提供することができる。
さらに，上記請求項１〜５のいずれか１項に記載の振動解析方法では，上記パラメータθ，ｅ3 ，ｅ6 ，ｅ12を，θ＝１．２５，ｅ3 ＝１／３，ｅ6 ＝１／６，ｅ12＝１／１２とすることによって，好ましい安定性，精度，数値減衰機能を得ることができる。
また，上記請求項７〜１０のいずれか１項に記載の発明によれば，外力ｆが与えられたときの上式（１）で表される多質点線形振動系の時刻歴応答を演算する際に，時刻ｔn における加速度ａn を用いて，時刻ｔn から時間τだけ経過した時刻ｔs における加速度ａs が，上式（２）の通り，時間τの２次式にて近似されるため，調整し得るパラメータ（θ，ｅ3 ，ｅ6 ，ｅ12）を増やすことができ，その結果安定性を確保しながら，必要とする周波数帯域では極めて精度がよく，しかも不必要な高周波モードを減衰させることができる振動解析装置を提供することができる。
【００１９】
また，上記請求項１１に記載の発明によれば，外力ｆが与えられたときの上式（１）で表される多質点線形振動系の時刻歴応答を演算する際に，時刻ｔn における加速度ａn を用いて，時刻ｔn から時間τだけ経過した時刻ｔs における加速度ａs が，上式（２）の通り，時間τの２次式にて近似され，時刻ｔn における速度ｖn を用いて，上式（３）の通り，時刻ｔs における速度ｖs が近似され，時刻ｔn における変位ｄn を用いて，上式（４）の通り，時刻ｔs における変位ｄs が近似され，これらの式（２）〜（４），及び上記運動方程式（５），（６）から構成される連立方程式（７），（７ａ）〜（７ｆ）から状態量ｂ，ｃが導き出され，次の演算周期の加速度ａn+1 ，ｖn+1 ，ｄn+1 が順次求められるため，調整し得るパラメータ（θ，ｅ3 ，ｅ6 ，ｅ12）を増やすことができ，その結果安定性を確保しながら，必要とする周波数帯域では極めて精度がよく，しかも不必要な高周波モードを減衰させることができる振動解析装置を提供することができる。
さらに，上記請求項７〜１１のいずれか１項に記載の振動解析装置では，上記パラメータθ，ｅ3 ，ｅ6 ，ｅ12をθ＝１．２５，ｅ3 ＝１／３，ｅ6 ＝１／６，ｅ12＝１／１２とすることによって，好ましい安定性，精度，数値減衰機能を得ることができる。
また，上記請求項１３〜１６のいずれか１項に記載の発明によれば，外力ｆが与えられたときの上式（１）で表される多質点線形振動系の時刻歴応答を演算する際に，時刻ｔn における加速度ａn を用いて，時刻ｔn から時間τだけ経過した時刻ｔs における加速度ａs が，上式（２）の通り，時間τの２次式にて近似されるため，調整し得るパラメータ（θ，ｅ3 ，ｅ6 ，ｅ12）を増やすことができ，その結果安定性を確保しながら，必要とする周波数帯域では極めて精度がよく，しかも不必要な高周波モードを減衰させることができる振動解析プログラムを記録したコンピュータ読み取り可能な記録媒体を提供することができる。
調整し得るパラメータが増加するのと同時に未知数も増加するから，それを同定するための方程式が必要になるが，例えば時刻ｔn+1 ＝ｔn ＋ｈ（演算周期）における上記運動方程式（５），及び時刻ｔs ＝ｔn ＋θｈ（θ＞１）における上記運動方程式（６）を用いればよい。
また，上記請求項１７に記載の発明によれば，外力ｆが与えられたときの上式（１）で表される多質点線形振動系の時刻歴応答を演算する際に，時刻ｔn における加速度ａn を用いて，時刻ｔn から時間τだけ経過した時刻ｔs における加速度ａs が，上式（２）の通り，時間τの２次式にて近似され，時刻ｔn における速度ｖn を用いて，上式（３）の通り，時刻ｔs における速度ｖs が近似され，時刻ｔn における変位ｄn を用いて，上式（４）の通り，時刻ｔs における変位ｄs が近似され，これらの式（２）〜（４），及び上記運動方程式（５），（６）から構成される連立方程式（７），（７ａ）〜（７ｆ）から状態量ｂ，ｃが導き出され，次の演算周期の加速度ａn+1 ，ｖn+1 ，ｄn+1 が順次求められるため，調整し得るパラメータ（θ，ｅ3 ，ｅ6 ，ｅ12）を増やすことができ，その結果安定性を確保しながら，必要とする周波数帯域では極めて精度がよく，しかも不必要な高周波モードを減衰させることができる振動解析プログラムを記録したコンピュータ読み取り可能な記録媒体を提供することができる。
さらに，上記請求項１３〜１７のいずれか１項に記載の振動解析プログラムを記録したコンピュータ読み取り可能な記録媒体では，上記パラメータθ，ｅ3 ，ｅ6 ，ｅ12をθ＝１．２５，ｅ3 ＝１／３，ｅ6 ＝１／６，ｅ12＝１／１２とすることによって，好ましい安定性，精度，数値減衰機能を得ることができる。
【図面の簡単な説明】
【図１】  本発明の実施の形態に係る振動解析方法を説明するためのフローチャート。
【図２】  本発明の実施の形態に係る振動解析装置の概略構成を示す機能ブロック図。
【図３】  上記振動解析装置を具体化するために用いられるコンピュータの概略構成を示す機能ブロック図。
【図４】  本発明による時刻歴応答演算のスペクトル半径と精度を説明するための図。
【図５】  本発明を適用した第一の具体例の演算結果を示す図。
【図６】  本発明を適用した第二の具体例において，対象となる振動系を説明するための図。
【図７】  本発明を適用した第二の具体例の演算結果を示す図。
【図８】  本発明を適用した第三の具体例において，対象となる振動系を説明するための図。
【図９】  本発明を適用した第三の具体例の演算結果を示す図。
【図１０】  従来の数値積分法の精度を説明するための図。
【符号の説明】
１…初期加速度算出手段
２…マトリクス算出手段
３…状態量算出手段
４…次時刻状態量算出手段
３０…コンピュータ
３７…記録媒体
Ｓ３…初期加速度算出工程
Ｓ４…マトリクス算出工程
Ｓ５…状態量算出工程
Ｓ６…次時刻状態量算出工程
【図１】
0003640843.tif 000018
【図２】
0003640843.tif 000019
【図３】
0003640843.tif 000020
【図４】
0003640843.tif 000021
【図５】
0003640843.tif 000022
【図６】
0003640843.tif 000023
【図７】
0003640843.tif 000024
【図８】
0003640843.tif 000025
【図９】
0003640843.tif 000026
【図１０】
0003640843.tif 000027