物体の振動

振動の原理
1. 振動現象
  バネやゴムひもなどの弾性体につながったおもりを引っ張って離すと、振動をします。おもりが弾性体の中心から離れると、中心方向に引っ張られます。中心まで戻ると、引っ張る力はなくなりますが、それまでの勢いで反対方向に行きすぎてしまいます。行きすぎると、また、元に戻す力が働き減速します。
  　振動を繰り返すと、摩擦や空気抵抗でエネルギを失うため、次第に振動は徐々に小さくなります。この振動を数式で表現し、これを数値解析して動きをシミュレーションします。
2. 力学モデル
  バネの中心を原点にします。おもりの重さをm、位置をxで表します。バネがおもりを引っ張る力Fはバネの長さに比例し、方向は逆ですから
```
    　F=-q・x
```
  となります。ここで、qはバネの強さです。おもりの加速度をαとすると、ニュートンの法則から
  　m・（dα/dt）　=　F = -q・x
  となります。（dα/dt）は時間に対する加速度の時間に対する変化率を意味しています。加速度は速度の変化率ですから　α=dv/dt　で、さらに速度は距離xの変化率ですから　v=dx/dt　です。従って、
  
  　　　m・(dx²/dt²) + q・x = 0　
  
  となります。この式は解析的に求められます。
  
  　　x(t)= A cos(ω・ｔ)　、w²=k/m
  　ここで、t=0で α=0、x=A としています。
  
  この動きを図にすると次のようになります。横軸は時間で、縦軸が中心からの距離(x)になります。
  q が大きくなると、戻りが速くなるため周期が短くなります。
  
  　　　　　　　
  　　　　　　　　　　　　q = 20　　　　　　　　　　　　　　　　　　　　　q = 10
3. 摩擦を考える
  この式では、エネルギーの損失項がないため、永久に振動を続けることになります。振動に対する摩擦などの項は、多くの場合速度に比例します。この摩擦に対する係数をpとしますと、摩擦の力は　
  　　　p・dx/dt　
  ですから、摩擦までを含めたｘを定める式（運動方程式）は次のようになります。
  
  　　m・(dx²/dt²) +p ・(dx/dt) + q・x = 0
  
  この式も、解析的に解くことができますが、かなり複雑な式になります。
  　しかし、時刻0におけるxの値とそのときの速度、加速度から、微小時間経過後の位置を予測できれば、t>0に対するxを順次求めることができます。計算で求めた、摩擦により減衰しながら振動する様子を下の図に示します。
  
  　　　　　　　　
  　　　　　　　　　p=0.2　　　　　　　　　　　　　　　　　　　　　　　p=0.6
運動方程式を計算する
1. オイラー法
  　時刻tにおける位置、速度、加速度、を x(t),v(t),a(t) とします。微小時間 h が経過後の位置、速度、加速度は
  hが十分小さければ、
```
　x(t+h) = x(t) + v(t)*h
　v(t+h) = v(t) + a(t)*h
　a(t+h) = a(t)-(q/m)x(t) + (p/m)・v(t)
```
  で求められます。この計算方法をオイラー法と呼びます。オイラー法は、h の間では「加速度は一定」であることを仮定しています。実際には、この仮定は正しくないので、h が大きくなると誤差が発生します。この誤差は、h と同じオーダーとなります。実際に計算すると、p=0 の場合でも定常的な振動せずに、振動が大きくなってしまいます。
2. ルンゲクッタ法
  　オイラー法よりも高い計算精度が約束される方法にルンゲクッタ法があります。この方法は、時刻　t+h の間の中間時刻の位置を予測し、この点を含めて高次の近似を行います。
  　ルンゲ・クッタ法は、次の時刻の値を求めるとき、これまでの位置、速度、から定まる直線の延長上ではなく、曲線上に予測します。このため、振動のような変化の多い場合でも比較的正確な計算ができます。理論的には、時間間隔hに対して、オイラー法が　o(h²) の誤差であるのにたいし、ルンゲクッタ法では、o(h⁵) と見積もられます。
  　ルンゲクッタ法では次のような計算をします。詳細は、数値計算法のテキストを参照してください。Xが位置、Yが速度、Hが時間の刻み幅となります。Dt の値を大きくすると、計算する時間の間隔が短くなり計算値も正確になります。Dt の値を10(h=0.1) 以上にすると、正しい結果が得られません。
```
  Y=0; // 初期速度
  X=1; //初期位置
  double Dt=40.0;H = 1 /( Dt);//時間の刻み幅
  for (int i = T0; i < TF; i++ ) {
    g.drawOval(i+xofs, (int)(hh-X*dh),1,1 );
    k1 = H * fnk( Y );
    l1 = H * fnl(Y, X);
    k2 = H * fnk(Y + l1 /2);
    l2 = H * fnl((Y + l1 /2), (X + k1 /2));
    k3 = H * fnk(Y + l2 /2);
    l3 = H * fnl((Y + l2 /2), (X + k2 /2));
    k4 = H * fnk(Y + l3);
    l4 = H * fnl((Y + l3), (X + k3));
    X += (k1 + 2 * k2 + 2 * k3 + k4)/6;
    Y += (l1 + 2 * l2 + 2 * l3 + l4)/6;
  }

public double fnk(double y) {
 return ( y );
}

public double fnl(double y, double x) {
  return ( -P*y-Q*x );
  }
```
プロジェクト
1. レイアウト
  p,qの値を設定するスライダを用意します。下の実行画面を参考にして下さい。
2. メソッド
  paint()の中で、計算をします。
```
int w = getSize().width;
　int h = getSize().height;
```
  で、このウインドウのサイズを取得します。
```
　P = scrP.getValue()/50.0;
　Q = scrQ.getValue()/2.0;
```
  で、パラメータを取り込み、数値の範囲を適正にするため、わり算をしています。
  　fillOval();は楕円を、drawspring()でバネを描画しています。rkgDeq() でルンゲクッタ法で　dh時間後の位置(rX)、速度(rY)、を計算しています。
3. 実行
  　開始ボタンを押すと、指定のバネ強度と摩擦係数でのシミュレーションを開始します。上段はおもりの運動を、下段は時刻に対するバネの長さをグラフで表示します。バネ強度と摩擦係数はスライダで変更可能です。クリアボタンで、下のグラフをクリアします。
ダウンロード
　このアプレットのソースはここにあります。また、このプロジェクト（Jbuilder）をダウンロードできます。次の行をクリックして、Bane.exeファイルを適当なフォルダに保存します。
　ダウンロード開始
このファイルは自己解凍型の圧縮ファイルです。このファイルを実行すると指定したフォルダに必要なファイルが生成されます。

振動の原理

振動現象

力学モデル

摩擦を考える

運動方程式を計算する

オイラー法

ルンゲクッタ法

プロジェクト

レイアウト

メソッド

実行

ダウンロード