随分と以前に自分のメモとして紙にまとめていたもの。
いつかは電子的に残しておきたいと思っていたので書きます。
さて、地動を受ける１質点系の運動方程式は、
と現す事がができます。ここに、
ある時刻 t における応答変位が x(t) であったとします。
すると Δt 秒後の変位は下図の台形を利用して、
となります。これを速度の計算にもあてはめると。
この（３）式に（２）式を代入すると、
となります。更に展開を進めると、
ここで第４項の 1/4 を β に置き換えたものを「ニューマークの β 法」と呼びます。
β は加速度変化の仮定をあらわすもので、
β = 0 ：衝撃加速度法
= 1/8：段階加速度法
= 1/6：線形加速度法
= 1/4：平均加速度法
としてあらわされます。
ところで（１）式より、
が成り立ちます。
（３）（４）（６）式で以下の３つは分かっているもの。
そして、以下の３つが次に知りたいものです。
式が３つで未知数が３つだから（３）（４）（６）式を連立方程式と考えて解けば未知数であるところの時刻 t+Δt の加速度、速度、変位が求まります。
あとは気合いの問題なのですが、式の代入まで書いているとわけが分からなくなるので、ここでは端折って、最後に求められる式を記載すると。
この t+Δt の時の加速度を（３）式に代入して t+Δt の時の速度を求め、次に（４）式から変位を算出すれば OK です。
初期条件として時刻 0 の加速度、速度、変位があれば次々と各時刻の応答が求まるのです。
しかし、（７）式は複雑ですねぇ。
ところがプログラムにするとそうでもないのです。
fortran で書いたサンプルはこんな感じ。
program dof1 implicit real*8(a-h,o-z) real*8 eq(10000) tau = 0.02 h = 0.05 w = 100.0 bk = 1000.0 ac = 0.0 ve = 0.0 di = 0.0 write(*,*) 'Number of step' read(*,*) nstep read(*,*) (eq(i),i=1,nstep) e1 = tau / 2.0d0 e2 = e1 * e1 g1 = ve + e1 * ac g2 = di + tau * ve + e1 * ac am = w / 980.665d0 cv = 2.0d0 * h * sqrt(am) * sqrt(bk) do i = 1 , nstep r = am + e1 * cv + e2 * bk ff = am * eq(i) + g1 * cv + g2 * bk ac = -ff / r ve = g1 + e1 * ac di = g2 + e2 * ac g1 = ve + e1 * ac g2 = di + tau * ve +e2 * ac t = float(i) * tau write(*,*) i,t,ac,ve,di end do stop end