VisualC++ と OpenGL を利用した仮想物理実験室
ラグランジュ運動方程式１：極座標を用いた単振子

拘束条件を取り扱う場合，一般化座標を導入したラグランジュ運動方程式を解くことで取り扱うことができます。 ラグランジュ運動方程式の一つ目の例として極座標を用いた単振子のシミュレーションを行います。

一般化座標θを用いて， 位置rを時刻 t で微分するとことで、速度位置vが得られます。

運動エネルギー T とポテンシャルエネルギーV 、ラグラジアンはそれぞれ次のように得られます。

これにより，ラグランジュ運動方程式から，一般化座標θの満たすべき方程式が導かれます。

シミュレーションを行うため，数値計算法としてルンゲクッタ法を利用します。 ルンゲクッタ法を用いるために，次のように変形します．

極座標を用いた単振子＋ルンゲクッタ法のアニメーション

極座標を用いた単振子＋ルンゲクッタ法のプログラム

【0-2-3】仮想物理実験室の構築 （ver1.2）を参照ください

構造体の定義

// ボールの定義
double ball_r = 4.0;//ボールの半径
struct BALL {
  double theta , theta_d;
  double l;
  double x, y, z;
};
const int N = 3;

BALL ball[N];
void SetUp(void){
  for(int i=0; i<N; i++ ){
    ball[i].l      = L/double(N) * double(N-i);
    ball[i].theta  = 30.0/180.0 * PI * (i+1);    //
    ball[i].theta_d  = 0,0;
  }
}

ルンゲクッタ法を利用した数値計算

/////////////////////////////////////////////
// ルンゲクッタ
double f1(double t,double theta,double theta_d, double l){
  return theta_d;
}
double f2(double t,double theta,double theta_d, double l){
  return -g/l*sin(theta);
}
void RungeKutta4(BALL *ball){
  double k1[2],k2[2],k3[2],k4[2];
  double theta   = ball->theta;
  double theta_d = ball->theta_d;
  double l       = ball->l;
  k1[0]=dt*f1(t,theta,theta_d, l);
  k1[1]=dt*f2(t,theta,theta_d, l);
  k2[0]=dt*f1(t+dt/2.0,theta+k1[0]/2.0,theta_d+k1[1]/2.0, l);
  k2[1]=dt*f2(t+dt/2.0,theta+k1[0]/2.0,theta_d+k1[1]/2.0, l);
  k3[0]=dt*f1(t+dt/2.0,theta+k2[0]/2.0,theta_d+k2[1]/2.0, l);
  k3[1]=dt*f2(t+dt/2.0,theta+k2[0]/2.0,theta_d+k2[1]/2.0, l);
  k4[0]=dt*f1(t+dt,theta+k3[0],theta_d+k3[1], l);
  k4[1]=dt*f2(t+dt,theta+k3[0],theta_d+k3[1], l);
  ball->theta    = ball-> theta   + (k1[0]+2.0*k2[0]+2.0*k3[0]+k4[0])/6.0;
  ball->theta_d  = ball-> theta_d + (k1[1]+2.0*k2[1]+2.0*k3[1]+k4[1])/6.0;

}