A..2 振り子の等時性

0 度から $179$ 度まで$1$ 度刻みの振幅に対して、 振り子の周期を計算する。 完全楕円積分であるから、 算術幾何平均(AGM)アルゴリズムを使って計算した (本当は GSL を呼び出せばよいわけで、半分道楽…いや、 GSL を使わず標準的な関数だけですむのは多少意味があるか？ 冪級数展開 (10) を使って計算する方が良いかな？)。

/*
 * toujisei.c --- 単振り子の等時性
 *   T=4√l/g K(k), k=sin θ0/2
 *   Th=2π√l/g
 *   http://www.math.meiji.ac.jp/~mk/labo/text/prog/toujisei.c
 */

#include <stdio.h>
#include <math.h>

double pi;

/*
 * 算術幾何平均アルゴリズムで完全楕円積分 K(k), E(k) を求める
 */
int kanzen(double k, double *K, double *E)
{
  int i;
  double a,b,an,bn,cn,anp1,bnp1,cnp1;
  double power2, s, AGM;

  if (k < 0 || k >= 1) {
    fprintf(stderr, "k must be in [0,1)\n");
    *K = *E = 0.0;
    return 1;
  }
  a = 1.0; b = sqrt(1.0 - k * k);

  an = a; bn = b; cn = sqrt(an * an - bn * bn);

  /* s = 2^(n-1) * cn * cn の和 (n=0,1,2,...) */
  power2 = 0.5;
  s = power2 * cn * cn;

  for (i = 0; i <= 10; i++) {
    anp1 = (an + bn) / 2; bnp1 = sqrt(an * bn); cnp1 = (an - bn) / 2;
    an = anp1; bn = bnp1; cn = cnp1;
/*    printf("%20.15f, %20.15f, %20.15f\n", an, bn, cn); */

    power2 *= 2;
    s += power2 * cn * cn;

    if (fabs(an-bn) < 2e-16)
      break;
  }
  if (fabs(an-bn)>1e-15) {
    fprintf(stderr,"k=%f: does not converge. i=%d, an=%f, bn=%f, |an-bn|=%e\n",
	    k, i, an, bn, fabs(an-bn));
    return 1;
  }
  AGM = an;
  *K = pi / 2 / AGM;
  *E = (a - s) * (*K);
  return 0;
}

int main()
{
  double theta,k,K,E;
  int deg;

  pi = 4 * atan(1.0);
  for (deg = 0; deg < 180; deg++) {
    theta = deg * pi / 180.0;
    k = sin(theta / 2);
    kanzen(k, &K, &E);
    printf("%2d %20.15f\n", deg, 2 * K / pi);
  }
  return 0;
}

Mathematica ではこんな感じ

  Plot[2 EllipticK[Sin[t Degree/2]^2] / Pi, {t,0,179}, PlotRange->{0,4}]