A..1 単振り子 vs 単振動

振幅 $\theta_0$ を指定して、単振り子の運動 $\theta(t)$ と、 同じ初期条件を課した単振動の運動 $\theta_{\mathrm{harmonic}}(t)$ を計算する。 Runge-Kutta 法のようなアルゴリズムを用いても良いが、 ここは GSL (GNU Scientific Library) のテストをかねて、 Jacobi の楕円関数を使って計算した。

/*
 * furiko.c --- 単振り子の運動
 *   θ=2 arcsin (k sn(√g/l t,k))
 */

#include <stdio.h>
#include <math.h>
#include <gsl/gsl_specfunc.h>

int main()
{
  double g,l,omega,theta0,k;
  int i,n;
  double t,dt,sn,cn,dn;

  g = 9.8;
  l = 1.0;
  omega = sqrt(g / l);
  scanf("%lf", &theta0);
  k = sin(theta0 / 2);

  dt = 0.01;
  for (i = 0; i <= 1000; i++) {
    t = i * dt;
    gsl_sf_elljac_e(omega * t, k*k, &sn, &cn, &dn);
    printf("%f %f %f\n", t, 2 * asin(k * sn), 2 * k * sin(omega * t));
  }
  return 0;
}

MacPorts で GSL をインストールした場合のコンパイル

gcc -I/opt/local/include furiko.c -L/opt/local/lib -lgsl

Mathematica ではこんな感じ

g=9.8; l=1; omega=Sqrt[g/l];
theta=1; k=Sin[theta/2];
Plot[{2k Sin[omega t], 2ArcSin[k JacobiSN[omega t,k^2]]}, {t,0,10}]

注意: これまで GSL では、

gsl_sf_elljac_e(omega * t, k, &sn, &cn, &dn)     (間違い！)

Mathematica では

JacobiSN[omega t,k]     (これも間違い！)

と書いて来たのですが、 GSL, Mathematica の楕円関数、楕円積分関係の関数は、 引数として $k$ でなく $m (=k^2)$ を使うようになっていることを教えていただいたので訂正しました (ご指摘に感謝します。2012/3/2)。