力学演習II 数値計算(3): 強制振動

2000-11-22


注意! 演習室のプリンタの用紙やトナーには限りがあります。 必要以上にたくさん印刷することは避けてください。 特に、WWWの文書(この文書を含む)は、 なるべくオンラインで見るようにしてもらえると助かります。


数値計算プログラムの高精度化

数値計算にEuler法をそのまま使うと、 dtをものすごく細かくしない限り、 大きな誤差を生じてしまう。 しかし、数値計算法を工夫すれば、精度は劇的に向上する。 たとえば、単振動の方程式

(d/dt)(p,q) =(-q, p)   (1)

2次の方法(改良Euler法またはHeun法)で解く場合、 dt=0.1とか dt=0.05程度の刻みで 10周期分計算しても、0.1%程度の誤差しか生じない。 計算結果を位相空間内の軌跡として図示すると、 ほとんど完全に閉じた曲線が得られる。

強制振動子の運動方程式

質量 m の質点が、 次の3種類の力を受けて、1次元的に運動しているとする:

この場合、質点の運動方程式は

m (d^2(x)/dt^2) = -kx - gamma dx/dt + F_ex

すなわち
m (d^2(x)/dt^2) + gamma dx/dt + kx = F_ex   (2)

となる。 方程式(2)を 適当に無次元化すると、 次のような形に書きなおすことができる:
(d^2(x)/dt^2) + beta(dx/dt) + x = cos(Omega t)   (3)

この方程式は、 betaOmegaという 2つの無次元パラメータを含んでいる。
実験や数値計算の結果を比較する場合には、 無次元パラメータ(今の例では betaOmega)が 同じであれば、 たとえ mkが 異なったとしても ``同じ実験条件''と見なすことができる(したがって、たとえば、 減衰が大きいか小さいかを論じる際には、 gammaの大小よりも betaの大小で考えるのが正しい)。 なお、本当は初期条件まで合わせてはじめて``同じ実験条件''となるのだが、 初期条件の影響が時間とともに減衰して消え去るような場合には、 初期条件の違いは重要ではない。

練習

以上の課題を、あさって(24日)の夕方までに提出せよ。

12月の課題

  1. 外力の振動数 Omegaの値を、 ちょうどbeta=0で 共振が起きるような値に設定する。 微小な減衰(beta=0.1)が存在する場合について、 系が定常状態に達するまで数値計算をおこない、 解を tの関数として図示せよ。 この結果を、beta=0の場合と比較せよ。 なお、初期条件は好きなように選んでよい。
  2. 減衰の大きさを beta = 0.2とする。 外力の振動数 Omegaの値を さまざまに変えて数値計算をおこない、 横軸にOmega, 縦軸に定常応答の振幅をとってグラフにまとめよ。 また、定常振幅と Omegaのあいだの関係式を理論的に求め、 数値計算の結果と一致することを確認せよ。
  3. 減衰を beta=1.5として、課題[2]と同様に さまざまなOmegaに対する定常応答振幅を求め、 理論と比較せよ。

これらの課題をレポートにまとめ、 12月18日の朝までに 事務室前のレポート箱に提出すること。


[力学演習のページに戻る]