単振り子の話
桂田祐史

2007 年 10 月 8 日, 2012/3/2, 2015/11/12 訂正

http://nalab.mind.meiji.ac.jp/~mk/labo/text/furiko/

はじめに1

たんしんし

紐の長さが ℓ である単振り子 (単振子, simple pendulum) の運動方程式は、m を重りの質量、g を重力加速度としてmℓθ′′ (t) = mg sin θ(t)となる。これから√(1)′′2θ (t) + ω sin θ(t) = 0,

ω :=gℓという微分方程式が得られる (ここに m は出て来ないので、単振り子の運動は質量には依らない、ということはすぐ分かる)。
微小振動 (|θ| ∼ 0) の場合は、sin θ は θ に近いので、単振り子の運動はθ′′ (t) + ω 2 θ(t) = 0
に従う単振動で良く近似できる、とされる (高校物理の相場)。単振動は (大学の理系の学科では)「常識的」で1 、一般解はθ(t) = C1 cos ωt + C2 sin ωt (C1 , C2 は任意定数)
であり、これは周期√T =2π= 2πωℓgを持つ周期関数である、と明解に解ける。この周期が g と ℓ だけで決まることから、「振り子の等時性が成り立つ (周期は振幅によらない)」ということにもなる。
しかし、単振動で近似するのはあくまでも近似である。どうして近似を用いた説明が多いのかというと、近似しないで解くのが難しいからである (後で見るように、解や周期を表すのに、楕円関数、楕円積分というものが必要になる2 )。この文書の目的は、本来の振り子の運動 (微分方程式 (1)の解) がどうなるか (近似を用いずに)、きちんと考えてみよう、ということである。
1

単 振 動 に つ い て は 、『 単 振 動 の 数 学 的 解 析 』(http://nalab.mind.meiji.ac.jp/~mk/keisanki2/
simple-harmonic-motion/) という解説を書いたことがある。
2
楕円積分、楕円関数については、内輪向けだが、『πノート (http://nalab.mind.meiji.ac.jp/~mk/labo/text/pi.pdf) というのを書いた。

1

単振り子は、単振動ではないとすると、本当に周期運動をするか？と心配となるが、次節でも証明するエネルギー保存則
1 ′ 2
θ (t) − ω 2 cos θ(t) = 定数
2
より大丈夫であることが分かる (戻ってきた時、前と同じ速さであり、減衰したりはしない3 。
210-1-2-40-224

図 1: 横軸 θ, 縦軸 θ′ . 閉曲線は、行ったり来たり。勢いが良いとぐるぐる (振り子じゃない！)。
ガリレオが発見したという『振り子の等時性』は、実は厳密には成立しない性質であるが (つまり周期は本当は振幅による)、どの程度成り立っているものなのかも、明らかにしよう (目で見て分かるようにする— 後の図 5)。

微分方程式の解析

2

前節で述べたことの繰り返しになるが、長さ ℓ の単振り子の微分方程式は

θ′′ (t) = −ω 2 sin θ(t).

(2)
ただし

√ω :=

2.1g.
ℓ

エネルギー保存則(2) に θ′ (t) をかけると、
θ′′ (t)θ′ (t) + ω 2 θ′ (t) sin θ(t) = 0.
すなわち[]
d 1 ′ 22θ (t) − ω cos θ(t) = 0.
dt 2

ゆえに

1 ′ 2
θ (t) − ω 2 cos θ(t) = E (E はある定数).
2
これは物理的にはエネルギー保存則に相当する (mℓE が全エネルギーを表す)。

(3)

3

厳密に証明するには、途中で止まらないことを言う必要があるけれど、それは明らかでしょう。注意が必要なのは、90 以下の触れの範囲の (0◦ < θ0 ≤ 90◦ として、−θ0 ≤ θ ≤ θ0 の範囲を行ったり来たりする) 「普通の」振り子以外に、90◦ より大きな触れの範囲の振り子 (90◦ < θ0 ≤ 180◦ — あまり現実的ではない？) や、鉄棒の大車輪のようなぐるぐる回る振り子 (これはリアリティがあるけれど、「振り子」という名前に相応しくない？) などの現象が方程式に含まれていることである。
◦2

予備的考察後のために、E を振幅 (θ の最大値) θ0 で書き直しておこう4 。θ = θ0 のとき θ′ = 0 であるから、

1 2· 0 − ω 2 cos θ0 = E.
2
これから

E = −ω 2 cos θ0 .
ゆえに (3) は次のように書き換えられる。

1 ′ 2
θ (t) − ω 2 (cos θ(t) − cos θ0 ) = 0.
2(4)ある時刻に θ = 0 となったとすると、(θ′ )2 = 2ω 2 (1 − cos θ0 ) = 4ω 2 sin2

θ′ = ±2ω sin

θ0であるから、2θ02である。つまりつりあいの位置にあるときの速度が振幅で表せたわけである。

2.2厳密解結果の式が簡単になるように、適当な初期条件を課した初期値問題の解を求めよう。
θ0 を 0 < θ0 < π なる定数として、(5)k := sin

θ02とおき、初期条件を(6)

θ(0) = 0,θ′ (0) = 2kωとする。前項の議論から振幅が θ0 であることが分かる。
この初期値問題の解は
θ(t) = 2 sin−1 (k sn (ωt, k))であることを示そう。
まずエネルギー保存則 (4) からdθ √ √= 2ω cos θ − cos θ0 .
dt(± をつけるべきかもしれないが、θ′ (0) > 0 であるから、t が小さいうちはこの微分方程式に従うはずである。)
これは変数分離形の微分方程式である。天下りであるが、(7)kz = sin4

θ2

ぐるぐる回る大車輪似の「振り子らしくない」振り子では、θ の最大値が存在しない。そういう場合は考えない、ということにしておく。この文書の議論は、色々な本に載っているものを適当にブレンドしたものではあるけれど、なかなかすっきりしないですねぇ。

3

と変数変換する。

(cos θ − cos θ0 =1 − 2 sin)2 θ2(− 1 − 2 sin

2 θ0)(2 θ0= 2 sin

222 θ− sin2)= 2(k 2 − k 2 z 2 )= 2k 2 (1 − z 2 ).
またk dz =1
θcos · dθ
22より

dθ =

2k dz
2k dz
=√.
θcos 2
1 − k2z 2(多分、(4) を、z(t) に関する微分方程式に変換する、とやればもっとすっきりした議論になるのでしょう。いつか書き直そう。)
これから∫ z∫ t∫ θ
1
1
12k ds1√ ·√√ √dθ = √·√t=dt =2cos θ − cos θ0
1 − k 2 s22ω2ω 0
2k 1 − s00
∫ z
ds
1
√
=.
2
ω 0
(1 − s )(1 − k 2 s2 )
すなわち

∫ z
ωt =

√0ds
(1 − s2 )(1 − k 2 s2 ).

これは、Jacobi の楕円関数 sn を用いると (ωt = sn−1 (z, k) であるから)

sn(ωt, k) = z =1
θsink2と解ける。ゆえに

θ(t) = 2 sin−1 (k sn(ωt, k)) .

(8)

これが初期値問題 (2), (6) の厳密解である。

単振動との比較3

前節の問題の単振動 (調和振動) 近似

θ′′ (t) = −ω 2 θ(t),θ(0) = 0,

θ′ (0) = 2kωの解は明らかに(9)

θ(t) = 2k sin (ωt) .

θ0 = 0.1, 0.2, 1 の場合に、単振り子の運動と、単振動を比較してみる。以下の図は横軸が t, 縦軸が θ で、単振動の解を赤、振り子の解を緑で描いたものである。

4

0.1
harmonic
theta=0.1

0.05

0-0.05-0.10246810

図 2: θ0 = 0.1: 振幅が小さいと単振動とほとんど同じ (上書きされて赤線が見えない)0.3
harmonic
theta=0.3

0.2

0.1

0-0.1-0.2-0.3
0246810

図 3: θ0 = 0.3: 振幅が大きくなると単振動からずれてくる

1
harmonic
theta=1
0.8
0.6
0.4
0.2
0-0.2
-0.4
-0.6
-0.8
-102468

図 4: θ0 = 1 (≒ 57.3◦ ) となるとかなりずれる510

「等時性」が成り立っているかを調べる

4

振幅が θ0 である単振り子の周期は (s が 0 から 1 に変化するのに要する時間の 4 倍だから)√∫
1 1dsℓθ0√
T =4×=4K(k), k := sin .
ω 0g2(1 − s2 )(1 − k 2 s2 )ただし K(k) は k を母数とする第 1 種完全楕円積分である:∫ π/2K(k) :=√0

これと単振動の周期

dt
1 − k 2 sin2 t.

√Th := 2πℓgとの比

2T= K(k)
Thπを具体的に計算してみよう。

T (θ0 = 1)
= 1.06633424557996 · · · ,ThT (θ0 = 45◦ )= 1.039973343196804 · · · ,Th
T (θ0 = 90◦ )= 1.180340599016096 · · · .
Thつまり θ0 = 45◦ でほぼ 4% しかずれていない。さすがに θ0 = 90◦ になると 18% ほどずれている。
[]( )2(()2)2
1
1·31·3·5π2461+k +k +k + ···
K(k) =
222·42·4·6であるから(10)T
2= K(k) = 1 +Th
π((( )2)2)2
1·3
1·3·5k2
124k +k +k6 + · · · ≒ 1 + .
22·42·4·6
4θ0 = 45◦ のとき、k = sin π8 = 0.38268 · · · であるから、1 +かな？(もう 1 項くらい取れば…)。

6k2= 1.0366 · · · . あまり良い近似でない
4θ0 (度)0102030405060708090100110120130140150160170179

T /Th
1.000000000000000
1.001907188143217
1.007669025791545
1.017408797595956
1.031340519130037
1.049782960623032
1.073182007149365
1.102144909639270
1.137492559923922
1.180340599016096
1.232229196737520
1.295339998876550
1.372880500618350
1.469819325894477
1.594446101177228
1.762203729503756
2.007507401244124
2.439362719673885
3.901065160389150

図 5 は、これをグラフにしたものである。90◦ < θ0 < 180◦ はあまり現実的な振り子の運動と言えないから、0◦ < θ0 < 90◦ の場合のみ注目すると、振り子の等時性はほどほどの精度で成立する、と言えそうである。
(振り子の等時性が成り立たないことは実験してみれば分かる、という人がいるけれど、かなり真剣に取り組まないとはっきりしないのではないか？と思う。)

A

プログラム

この文書に示した数値計算結果を得るのに用いたプログラムを紹介する。
趣味色が強い。情報処理教室で再現してもらうためには、Mathematica でやるのがよいか？と考えて試しに書いてみたら、1 行で済んだりして…

A.1

単振り子 vs 単振動

振幅 θ0 を指定して、単振り子の運動 θ(t) と、同じ初期条件を課した単振動の運動 θharmonic (t) を計算する。Runge-Kutta 法のようなアルゴリズムを用いても良いが、ここは GSL (GNU ScientificLibrary) のテストをかねて、Jacobi の楕円関数を使って計算した。
/** furiko.c --- 単振り子の運動*
θ=2 arcsin (k sn(√ g/l t,k))*/#include <stdio.h>#include <math.h>#include <gsl/gsl_specfunc.h>7

4"toujisei.txt"
3.5

32.5

2

1.5

10.5

0
02040

60

80

100120140160180

図 5: 振り子の等時性: 振幅で周期がどう変るか
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 では、8gsl sf elljac e(omega * t, k, &sn, &cn, &dn)

(間違い！)

Mathematica では
JacobiSN[omega t,k]

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

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;

9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}]
BRunge-Kutta 法で解く

今さら数値計算する気がなくなる…まあ、気が向いたら。と書いたのですが、歳は取りたくないもので、ずっと昔、授業 (『振り子の運動』5 ) でやったことがありました。その計算でも、確かに振幅 1 ラジアンで 6.63% 周期が長くなっていました。問題として、完全楕円積分を級数展開で計算しろ、なんて書いてある (笑)。ともあれ、全然違う方法で計算して得た数値が良く一致したので、とりあえず安心。

C

周期を数値積分で求める

ははみまかりける年某月某日、学生 Y の襲撃に会う。何でも振り子の周期を求めるために、端点に特異性のある定積分を計算するように指示されたけど、自分でネットを徘徊して、変数変換によって特異性を除去してから台形則で計算したら (おおっ)、粗い分割で非常な高精度を達成できた、という。実質的に第 2 種完全楕円積分∫ π/2K(k) =√0dt
1 − k 2 sin2 tの計算に帰着させたらしい。
5

http://nalab.mind.meiji.ac.jp/~mk/syori2-1992/syo2-92-06.pdf

10台形則で高精度…あれだな、きっと。ちょっとやってみよう。まずは Mathematica で。
k = 1/2 のときの K(k) の値

f[t_,k_]:=1/Sqrt[1-k^2 Sin[t]^2]
k=1/2
g=Plot[f[t,k],{t,-Pi,Pi}]
被積分関数の周期は π か。

Integrate[f[t,k],{t,0,Pi/2}]→ EllipticK[1/4] と表示されるけれど、これは Mathematica 方言で、本当は K(1/2) ということだな。

N[%,50]→ 1.6857503548125960428712036577990769895008008941411 となる。
ちなみに Mathematica で数値積分をするのは、NIntegrate[] だ。

NIntegrate[f[t, k], {t, 0, Pi/2}, AccuracyGoal -> 50, WorkingPrecision -> 60]
さて、それで台形則をしてみる。

TR[N_] := Block[{i, h, S}, h = (Pi/2)/N; S = 0;h (f[0, k]/2 + Sum[f[i*h, k], {i, 1, N - 1}] + f[Pi/2, k]/2)]
N[Table[TR[n], {n, 2, 8}], 16]%-EllipticK[k^2]{0.000024522126040, 1.04290737*10^-7, 4.68015*10^-10, 2.165*10^-12,
1.0*10^-14, 0.*10^-16, 0.*10^-16}なるほど、これは凄い。
11実質的に、有名な「解析的周期関数の 1 周期区間の積分は、台形則で非常に高精度に計算できる」という数値積分の常識の実例に遭遇したわけだ。
C のプログラムを示せば、

12
period.c#include <stdio.h>#include <math.h>double k = 0.5;
double sqr(double x) { return x * x; }double f(double x){
return 1 / sqrt(1- sqr(k * sin(x)));}int main(void){int i, n;double h,s,pi;
pi = 4 * atan(1.0);
for (n = 2; n <= 8; n++) {h = pi / 2 / n;s = (f(0) + f(pi/2)) / 2;
for (i = 1; i < n; i++)s += f(i*h);s *= h;printf("%2d %20.15f\n", n, s);
}
}





$ gcc period.c$ ./a.out
2
1.685774876938636
3
1.685750459103333
4
1.685750355280611
51.685750354814761
61.685750354812606
7
1.685750354812596
8
1.685750354812596$

ふむふむ (区間を 7 等分で倍精度一杯の精度だ)。AGM とどっちが速いんだろう。
それにしても、自力でそういう工夫をして、高精度な結果が得られたことを不思議に思って追求して質問に来るとは、大したものだと思う。

13