対話・非線形振動（その１）
K E N Z OU
2008 年 9 月 28 日
♣ まだ残暑の厳しさが残る初秋のある日，コニーが少し日焼けした顔で，お気に入りのスニーカーを履いてＫ氏を訪ねて
きた。コンニチワ∼Ｋさん，まだ残暑厳しいわねェ∼
• Ｋ氏：Ｏｈ！
！コニー久しぶりだねぇ∼。お彼岸も過ぎたけど，残暑はまだ厳しいようだね。 デートでいろいろ忙し
そうじゃないか，この前キャサリンが言っていたよ。大いに遊んで大いに学ぶ，結構なことだ。
• コニー：ハイ，おかげさまで。しかし，暑さは異常だわね。体温以上の日が続くなんて地球温暖化の影響ね。そう
いえば最新オゾンホールの話題がホットでなくなったけど紫外線対策はバッチリしたつもりよ。
• Ｋ氏：だろうね，シミがでるとまずいからね，
，
，オッと失礼。
• コニー：いいえ，私は色白だから，ことのほかシミ対策には気を配っているわ。ところで今日お伺いしたのは，最
近，非線形振動というものに興味が湧いて，正確に言うと非線形波動なんだけど，そのとっかかりとしてます非線
形振動あたりのお話がいろいろ聞ければなぁと思ってきたの。
• Ｋ氏：非線形振動か。非線形という奴は数学的にうるさくなるからね，しかしいろいろ面白いことも確かだね。僕
も詳しくは知らないけど，一緒に勉強するつもりで聞いてくれるのなら話をしてもいいと思うけど。。。
• コニー：日頃のＫさんに似合わず少し弱気ね。暑さのせいかしら。数学的にややこしいところは適当に飛ばしても
らっていいのよ。そういうところは物理的な意味がつかめればいいの。一から十まで聞けるとは思っていないわ。
まぁ，いろいろ聞きたいことはあるのだけど，まずは非線形振動のイントロ的なところのお話が聞ければいいと思っ
ているの。
• Ｋ氏：そうなんだ。例えば単振り子で振幅が大きくなるとこれは線形振動ではなくなるけど，そのあたりの話を今
日はしようか。こいつは楕円積分というのが登場してくるね。大抵のテキストはそこまで踏み込んで書かれていな
いけど，格段難しいものでもないよね。初等関数で表すことができないから数値計算で値を求めることになるんだ
けど，それが障害になっていりのかなぁ。
• コニー：楕円積分については詳しくは知らないけど，わかりやすく説明していただければいいわ。非線形波動にい
くと楕円関数というのがドシドシでてくるとかいくことだから，そんなことで腰を引くわけにはいかないわね。
• Ｋ氏：そうだね。その意気込みが大事だね。それじゃ，いきなり非線形振動というのもなんだから，急がば回れじゃ
ないけど，先ずは振動現象の一般的なおさらいをしてからということでどうだい。
• コニー：結構よ。振動現象の基本的なところの復習をしてからっということね。よろしくお願いします。
• Ｋ氏：振動に関するテキストは本屋にいろいろでていると思うけど，僕が参考にするのは有山正孝「振動・波動」，
吉岡大二郎「振動と波動」といったところだ。他にもいろいろあるから気に入ったテキストで僕の話の不足分は補っ
てね。それでははじめようか。
1
単振動とエネルギー保存則
1.1
単振動と等時性
• Ｋ氏：最初は単振動の話だけど，これはよく物体が原点からの距離に比例する力 F = −kx を受けている
場合の運動だ。ここで距離 x の 2 乗や 3 乗などでなく，素直に距離の１乗に比例している力ということが
大事なんだ。例えばバネを自然の長さから伸ばしたり縮めたりすると、元の長さに戻ろうとするいわゆる
復元力がはたらくけど，この 力と変形量との間には比例関係がある ということをイギリスの物理学者ロ
バート・フックが 1678 年に発見し，これをフックの法則とよんでいることはよく知っていると思う。
1
この系の運動方程式は次の 2 階線形微分方程式で表され，
m¨
x = −kx,
(k > 0)
(1.1)
これは 1 次元調和振動子と呼ばれる。この方程式を解いて変位 x の時間変化は
x(t) = A cos(ωt + α)
(1.2)
と求まるんだけど，ここで微分方程式の線形性がどのように効いているのかみるために，具体的に（1.1）
を解いてみよう。解き方はいろいろあるけど，ここでは因数分解を使ったやり方をやってみる。
d/dx = D という微分演算子を使って（1.1）を因数分解すると
d2
k
x+ x=
dt2
m
D2 +
k
m
x=0
.
..
D+i
...
D = ±i
...
x = A1 exp(iωt),
k
m
x=
or
k
m
D+i
k
m
D−i
k
dx
−→
= ±iωx ,
m
dt
(ω =
k
m
D−i
x=0
x=0
k/m)
x = A2 exp(−iωt)
となって基本解が求まる。A1 , A2 は任意の定数だ。一般解はこの解の線形和（重ね合わせ）で表されるから
x = A1 exp(iωt) + A2 exp(−iωt)
(1.3)
とおける。（1.3）はオイラーの公式（eiθ = cosθ + i sin θ）を使うと
x = (A1 + A2 )cosωt + i(A1 − A2 )sinωt
(1.4)
となる。左辺の x は実数だから，右辺も実数でないとまずい。つまり，(A1 + A2 ), i(A1 − A2 ) がともに実
数であれば等号が成り立つ。そのためには A1 , A2 が互いに共役な複素数（A1 = A¯2 , A2 = A¯1 ）であれば
よい。そこで ξ, η を実数として
A1 =
1
(ξ + iη) ,
2
A2 =
1
(ξ − iη)
2
(1.5)
とおいてやると,（1.4）は
x = ξ cos ωt − η sin ωt =
ξ
ξ2 + η2
ξ2 + η2
cos ωt −
η
ξ2 + η2
sin ωt
= A(cos α cos ωt − sin α sin ωt)
ξ 2 + η 2 , α = tan−1 (η/ξ))
= A cos(ωt + α) (ただし，A =
となって（1.2）がでてくる。
ここで用語を整理しておこう。
A：
ω：
α：
ωt + α：
f：
T：
角振動数 [rad/sec]・
・
・単位時間に何回振動するかを表す。
初期位相 [rad]
振幅・
・
・定数となっていることに留意。
位相 [rad]
周波数 [Hz]（f = ω/2π)
周期 [sec]（T = 1/f = 2π/ω)
(1.6)
上の話はバネの場合たけど，壁に吊り下げた振り子の場合も，振れ角度 θ が小さい場合は（1.1）と同じ微
分方程式
θ¨ = −ω 2 θ,
(ω =
g/l, g : 重力加速度, l：振り子の長さ)
2
(1.7)
で表され，同様に解が求まることはよく知っていることと思う。
この系で重要なポイントは，周期 T には振幅 A が含まれていない という点だ。これが単振動の大きな特
長で，
「振り子の振幅 A がどんなに大きくなっても，小さくなっても，振り子の周期 T には無関係」とい
ういわゆるガリレオ1 が発見した等時性をあらわしているということだね。
• コニー：ピサ大聖堂の天井から釣り下がっているランプの揺れを見て，振り子の周期は振幅には無関係で
振り子の長さに比例する (T = 2π l/g) という「振り子の等時性」を発見したという有名なお話ね。
ところで，位相が無次元量というのはわかるけど，無次元量に rad という ”単位 ”が付くのはなにかおか
しな感じね。
• Ｋ氏：そうだね。え∼っと，rad というのは，[rad] ＝「弧の長さ÷半径」で定義され，長さを長さで割っ
ているから無次元数なんだね。そして rad を無次元の組立単位とするとなっているんだね。
• コニー：そうなの。無次元量といっても大きさを持っているものね。
1.2
エネルギー保存則
• Ｋ氏：単振動の系のエネルギーについて調べてみよう。運動エネルギー Ek とポテンシャルエネルギー（弾
性エネルギー）はそれぞれ
1
1
Ek = mx˙ 2 , U = k(x − x
¯)2
(1.8)
2
2
と表される。ポテンシャルエネルギーの基準点を x = 0 におくと U = 21 kx2 となるので，この系の全エネ
ルギー W は
1
1
mx˙ 2 + kx2
2
2
と表される。（1.1）の両辺に x˙ をかけて変形していくと
W = Ek + U =
mx¨
˙ x − k x˙ =
...
(1.9)
1 d
mx˙ 2 + kx2 = 0
2 dt
1
1
mx˙ 2 + kx2 = W = const
2
2
(1.10)
となって，全エネルギーが保存されることがわかる。W は（1.9）より
1
1
W = mω 2 A2 sin2 (ωt + α) + kA2 cos2 (ωt + α)
2
2
1
2 2
= mω A
2
(1.11)
で，振幅 A の 2 乗に比例する。さて，運動エネルギー，ポテンシャルエネルギーの 1 周期にわたっての時
間平均を求めてごらん。
• コニー：はい，運動エネルギー，ポテンシャルエネルギーの平均値をそれぞれ < Ek >, < U > とすると
< Ek >=
< U >=
T
1
T
1
T
0
T
0
T
1
1
mx˙ 2 dt =
2
T
1 2
1
kx dt =
2
T
... < Ek >=< U >
0
T
0
1
1
1
mω 2 A2 sin2 (ωt + α)dt = mω 2 A = kA2
2
4
4
1 2
1
1
kA cos2 (ωt + α)dt = mω 2 A = kA2
2
4
4
運動エネルギーの平均値とポテンシャルエネルギーの平均値は等しいということになるわね。
1
Galileo Galilei 1564-1642
3
1.3
任意のポテンシャル下での微小振動
• 任意のポテンシャルの下での極小値付近における微小振動の一般論を展開しておこう。ポテンシャルを
U (x) とし，x = x0 で極小値（平衡点）を持つとする。振幅の変位 x は小さいので U (x) を x0 の周りでテ
イラー展開してやると
U (x) = U (x0 ) + (x − x0 )U (x0 ) +
(x − x0 )2
(x − x0 )3
U (x0 ) +
U (x0 ) + · · ·
2!
3!
(1.12)
x = x0 でポテンシャルは極小値をとるので U (x0 ) = 0 かつ U (x0 ) > 0 だね。次に平衡点からのズレ
(x − x0 ) は微小として 3 次以降の項を無視すると
U (x)
U (x0 ) +
U (x0 )
(x − x0 )2
2
(1.13)
と書ける。力は −U (x) から導かれるので運動方程式は
m¨
x=−
d
U (x)
dx
(1.14)
となる。U (x) に（1.13）を入れると2
m¨
x = −U (x0 )(x − x0 )
(1.15)
が得られるね。ここで ω 2 = U (x0 )/m, x − x0 = x とおいてやると
x
¨ = −ω 2 x
(1.16)
これは単振動の方程式となる。つまり，極小値をもつ任意のポテンシャルの下での微小な振動運動は単振
動（調和振動子）で表されるということなんだ。
• コニー：なるほど，調和振動子が物理で大事にされる理由がよくわかるわ。ところで任意のポテンシャル
といわれたけど，力 F が保存力で F = −∇U の関係にあるポテンシャルでなければならないのね。
• Ｋ氏：うん，いわゆる保存力場ということだね。それじゃ一つ演習問題として，ポテンシャルの形状がワ
イングラスの底のような形をした U = k cosh x を考えてみよう。Fig.1 のイメージだね。
y
F ig.1
y = cosh x = (ex + e−x )/2
ex = 1 + x + (1/2)x2 + · · ·
e−x = 1 − x + (1/2)x2 + · · ·
0
x
運動方程式は（1.14）から
m¨
x=−
d
d
k d x
k
U (x) = −k
cosh x = −
(e + e−x ) = − (ex − e−x )
dx
dx
2 dx
2
(1.17)
だね。変位 x は微小とするから，ex ，e−x を平衡点 x = 0 の周りでテイラー展開（マクローリン展開）す
ると
ex = 1 + x + (1/2)x2 + · · · ,
2
e−x = 1 − x + (1/2)x2 + · · ·
U (x0 ), U (x0 ) は値の定まった定数となっていることに注意。
4
したがって平衡点まわりの微小振動の運動方程式は
k
m¨
x = − (1 + x + · · · − 1 + x − · · · ) = −kx
2
となる。この方程式は角振動数 ω =
2
(1.18)
k/m の単振動を表しているね。
周期とエネルギーの関係
• Ｋ氏：単振動の場合，周期 T は角振動数より T = 2π/ω と求まったが，ここではエネルギー保存則を使っ
て周期用いて T を求めてみよう。全エネルギー W は
W =
1
m
2
dx
dt
2
+ U (x)
これから
dx
=
dt
...
2{W − U (x)}/m
x
m
2
t=
dx
W − U (x)
(2.1)
+ const
が得られる。このやり方をエネルギー積分方式と呼ぶことにしよう3 。
（2.1）の右辺の W − U (x) は運動エ
ネルギーを表しているから，仮に x = A1 と x = A2 の間で振動しているとすると，その両端では速度がゼ
ロとなるから，V (x) を x の関数として
W − U (x) = (x − A1 )(A2 − x)V (x),
(A2 > A1 )
(2.2)
と書ける。質点は x = A1 から x = A2 へ動くのに要する時間は半周期 T /2 の時間だから，
m
2
T =2×
A2
dx
W − U (x)
A1
√
＝ 2m
A2
dx
(x − A1 )(A2 − x)V (x)
A1
(2.3)
ということになる。そして
x=
1
1
(A1 + A2 ) − (A2 − A1 ) cos φ ,
2
2
(0 ≤ φ ≤ π)
と変数変換すると
dx = (1/2)(A2 − A1 ) sin φdφ
1
W − U (x) = (A2 − A1 )V (x) sin φ ,
2
(A2 > A1 )
となるので（2.3）は
T =2×
m
2
A2
A1
dx
W − U (x)
=
√
π
2m
0
dφ
V (x)
(2.4)
となる。
3
サイクロイド振子
• Ｋ氏：振幅が小さい単振り子の場合は等時性が成立すること見てきたけど，振幅が大きくなれば等時性は
成立しなく。しかし，サイクロイド振子とかホイヘンス（Huygens) 振子と呼ばれている，鉛直面内のサイ
クロイド曲線上で振動する質点は振幅が大きくなっても一定の周期の単振動をするんだ。これを完全等時
性と呼んでいる。
3
ここだけの方言ですからご注意あれ！
5
y
F ig.2
dx2 + dy 2 = 2a cos(θ/2)dθ
ds =
2a
ds
−πa
mgy
θ
0
s=
θ
0
y=
1
2a
√
ds = 4a sin(θ/2) = 2 2ay
s
2
2
x
πa
振動の周期 T は（2.1）より求められる。ポイントは質点の速度 s˙ を求めることだね。
s
dt =
θ
ds
m/2
=
W − U (s)
m/2
2a cos(θ/2)dθ
W − U (s)
θ
=
m/2
2a cos(θ/2)dθ
s˙
(3.1)
全エネルギー W は最上点でのポテンシャルエネルギー Umax に等しいので
W = Umax = 2mga
途中点 x でのポテンシャルエンルギー U (s) は
U (s) = mgy = mg
1
2a
s
2
2
運動エネルギー Ek はしたがって
Ek =
1
m
2
ds
dt
2
= W − U (x) = 2mga − mg
1
2a
s
2
2
これから質点の速度 s˙ を求めると
s˙ = 4ga −
g 2
s
4a
1/2
= 4ga 1 − sin2
θ
√
= 2 ga cos
2
θ
2
(3.2)
となる。これを（3.1）に入れ，積分範囲を 0 ≤ θ ≤ 2π にとると T /2 が求まる。積分を実行すると
2π
T /2 =
0
.
..
T = 4π
2a cos(θ/2)
dθ =
√
2 ga cos(θ/2)
a
g
2π
dθ = 2π
0
a
g
a
g
(3.3)
が得られる。周期には振幅が含まれていないので 完全等時性が成立する ことがわかる。
また，角振動数 ω は
ω=
2π
1
=
T
2
g
a
(3.4)
となる。尚，ここのサイト http://members.jcom.home.ne.jp/sheer-heart-attack/に「サイクロイ
ド曲線に束縛された質点の運動」というレポートが掲載されているので，是非見ておけばいいと思うよ。
• コニー：面白そうね。
4
非線形振動
4.1
振幅の大きい振り子
• Ｋ氏：天井からつり下げた振り子の運動方程式は
g
θ¨ = − sin θ = −ω 2 sin θ,
l
6
(ω =
g/l)
(4.1)
となる。振幅が小さい場合
θ3
θ5
+ −
(4.2)
3!
5!
と展開して 3 次以上の項を無視すると，単振動の方程式（θ¨ = −ω 2 θ）となるよね。ここではこの近似が成
sin θ = θ −
立しない振幅の大な場合を取り扱う。
F ig.3
y
運動方程式
運動エネルギー
2
l
θ
2
Ek = (1/2)(x˙ + y˙ )
m¨
x = mg − T cos θ
m¨
y = −T sin θ
ポテンシャルエネルギー
T
U = −mgl cos θ
=⇒ θ¨ = −(g/l) sin θ
束縛条件
¨ = −l(φ¨ sin φ + φ˙ 2 cos φ)
x = l cos θ → x
y = l sin θ → y¨ = l(φ¨ cos φ − φ˙ 2 sin φ)
x
mg
系の全エネルギー W は W = Ek + U だから
W =
1 2 ˙2
ml θ − mgl cos θ = ml2
2
1 ˙2
θ − ω 2 cos θ
2
(4.3)
運動方程式（4.1）の両辺に θ˙ をかけて整理すると
d
θ˙θ¨ = −ω 2 sin θθ˙ −→
dt
1 ˙2
θ − ω 2 cos θ
2
=0
(4.4)
となって，この系は保存量 E
1 ˙2
θ − ω 2 cos θ
2
をもつことがわかる。E に ml2 をかけたものが先ほどの（4.3）で，エネルギー保存則が成り立っている。
E=
（4.4）の微分方程式を t = 0 で θ = θ0 , θ˙ = 0 という初期条件で解くと
θ˙2 = 2ω 2 (cos θ − cos θ0 )
... θ˙ = 2ω
sin2 (θ0 /2) − sin2 (θ/2)
(4.5)
が得られる。ただし最後の変形では三角関数の公式 cos φ = 1 − 2 sin2 (φ/2) を使った。θ = 0 から θ = θ0
となる時間を計算しよう。ここで θ0 は t = 0 での振れ角だったから最大振れ角ということになるね。だか
ら所要時間は T /4 となる。積分すると
t=
1
T
=
4
2ω
θ0
0
1
sin2 (θ0 /2) − sin2 (θ/2)
dθ =
1
ω
π/2
0
1
1 − a2 sin2 φ
dφ
(4.6)
1 − a2 sin2 φ )dφ だね。φ の
積分範囲は，θ と φ の関係式より θ = 0 のときは φ = 0，θ = θ0 のときは φ = π/2 となる。したがって周
期 T は次式で与えられるね。
となる。ここで sin(θ0 /2) = a, sin(θ/2) = a sin φ とおいた。また dθ = (2/
T =
4
ω
π/2
1
2
1 − a2 sin φ
0
π/2
ただし K(a) =
dφ =
1
1 − a2 sin2 φ
0
7
dφ
4
K(a)
ω
(4.7)
(4.8)
（4.8）の積分はいわゆる第 1 種の完全楕円積分（complete elliptic integral）4 と呼ばれるモノで，残念な
がらこの不定積分は初等関数で表すことはできない。そこで被積分関数をテイラー展開して積分するとい
う手をとる。テイラー展開すると
∞
31 4 4
531 6 6
(2n − 1)!! 2n 2n
1
= 1 + a2 sin2 φ +
a sin φ +
a sin φ + · · · =
a sin φ
2
2
2
4
2
6
4
2
(2n)!!
1 − a sin φ
n=0
1
となる。右辺最後の式で !! と 2 重階乗がでてきたけどこれは n から一つおきに下った数の掛け算で例えば
6!! = 6 × 4 × 2 = 48 という調子だね。積分公式
π/2
π (2n − 1) · (2n − 3) · · · 3 · 1
π (2n − 1)!!
=
2
n · (n − 2) · · · 4 · 2
2
n!!
sin2n θdθ =
0
を使うとこの積分は容易にできて
π/2
1
K(k) =
1 − a2 sin2 φ
0
dφ =
1
92 4
25 6
1 + a2 +
a +
a + ···
4
64
256
π
2
これから周期 T は
T =
4
2π
1
92 4
25 6
K(k) =
1 + a2 +
a +
a + ··· ,
ω
ω
4
64
256
a = sin
θ0
2
(4.9)
が得られる。これを眺めると周期 T は最大振れ角 θ0 に依存することがわかる。つまり，振幅が大きくな
れば a の高次の寄与が効いてくるので 等時性は成立しない ということだね。もっとも a
1 という条件
下では a の 2 次以上の項が無視できて等時性が成立する
振幅が大きくなった場合，周期 T と単振動の周期 Th の比をを見積もってみると
Th = 2π/ω, T = (4/ω)K(k) −→ T /Th = (2/π)K(k)
下表の左側は近似式を使わずに求めた値で，右側は（4.9）の近似式の第 3 項目までをとって求めた値だ。
近似式の第 3 項までとれば θ0 が 50 度くらいまではほぼ同じ値を示すことがわかるよね。
θ0
T /Th
θ0
T /Th
θ0
0
1.01348 30
1.01738
0
1.00000 30
1.01741
5
T /Th
θ0
T /Th
1.00048 40
1.03134
5
1.00048 40
1.03117
10 1.00191 50
1.04978
10
1.00191 50
1.04914
15 1.00430 60
1.07318
15
1.00430 60
1.07129
20 1.00767 80
1.13749
20
1.00767 80
1.12730
25 1.01203 90
1.18034
25
1.01202 90
1.16016
図は振幅 θ0 と T /T h の関係で θ を 179 °までとったものだ。これらの結果から θ0 ≤ 20°までは等時性が
成り立つとみてよいと思うけど，どう思う。
Π2
T Th
1
2Π
0
T Th
1
Φ, a sin
a2 sin Φ
2
Θ0
2
4
3.5
3
2.5
2
1.5
0.5
4
1
1.5
2
2.5
3
Θ0
Karl Gustav Jakob Jacobi(1804 年―1851 年）ドイツの数学者。1829 年に楕円関数論を発表。楕円の弧長を求める時に発見したと
かいわれている。
8
• コニー：そうねぇ，テーブルから判断して θ0 ≤ 5°がまず無難なところね。あと Ｋ氏 さんが言われるよう
に短時間の範囲であれば 20 °未満もいけそうといったとかしら。
• Ｋ氏：そうだね。ところで，いま見たように振れ角や振幅が大きくなると変位の 2 乗や 3 乗という高次の
項を無視できなくなり，振動は単振動（調和振動）ではなくなるんだね。いわゆる非線形振動（非調和振
動）ということになり，このような場合は厳密には等時性は成り立たないし，解の重ね合わせも成立しな
い。単振動のところでやったように解 x はというわけにはなかなかいかないんだね。
参考：Mathematica での計算例
For[θ=0Degree, θ<=90Degree, θ=θ+5Degree,a=Sin[θ/2]^2;Ans=N[(2/Pi)*ElliptickK[a]];
Print[{θ,Ans}]]
Plot[EllipticK[Sin[θ/2]*(2/Pi),{θ,0,179Degree},PlotStyle->Thickness[0.004],
Ticks->{Table[θ Degree,{θ Degree,0Degree,180Degree,20Degree}],Automatic},
AxesLabel1->{"θ0 ",T/Th},PlotLavel->StyleForm["T/Th=(2/π)K(a),a=sin(θ0 /2),
(π/2)
√ 1 2 dφ”,
K(a)= 0
2
1−a sin φ
FontFamily->"Times",FontSize->14,FontWeight->"Bold"]]
4.2
非線形のバネ
バネの伸びと力との関係がフックの法則に従わない場合，具体的にはバネの力が
F = m¨
x = −(αx + βx3 ), （α > 0, β = 0）
(4.10)
で表される場合を考えよう。この方程式はダフィング方程式（Duffing equation）と呼ばれている5 。β は
非線形度を表すパラメーターで β > 0 のものをハードニング型（漸硬バネ），β < 0 のものをソフトニン
グ型（漸軟バネ）と呼ばれる。この非線形の復元力を図示すると Fig.4 のようなイメージだね。
■定性的把握
ところで，非線形バネの物理的イメージをまず掴むために，まず質点がポテンシャル U (x, β) = (1/2)αx2 +
(1/4)βx4 のもとで運動している場合を定性的に考えてみよう。質点の運動方程式は次のようになる。
F = m¨
x=−
F ig.4
d
U = −αx − βx3
dx
(4.11)
復元力 F
β>0
β=0
F = m¨
x = −(αx + βx3 ),
(α > 0)
W = Ek + U
β<0
変位 x
0
U = (1/2)αx2 + (1/4)βx4
Z x
p
dx
p
t = ± m/2
W − U (x)
振幅が小さい場合は単振動となるが，振幅の大きくなるとこの第 2 項は無視できない。質点は x = −A から x =
5
A の間をいったりきたり振幅 A の振動運動をしており，周期 T が非線形パラメーター β の大きさによっ
てどのように変わるか調べてみることにしよう。
（4.1
）で sinθ = θ − θ 3 /3! + θ 5 /5! − · · · 展開し第 2 項まで取ると θ¨ = −ω 2 (θ − θ 3 /6) となって同様の方程式がでてくる。
9
F ig.5
U
U = (1/2)αx2 + (1/4)βx4 ,
(β > 0)
Umax = (1/2)αA2 + (1/4)βA4
(β < 0)
−A
x
v(x) A
0
変位位置 x における質点の速度を v(x) とすると距離 dx 進むに要する時間 dt は dt = dx/v となる。した
がって x = −A から x = A まで進むに要した時間を t とすると，これは周期 T の 1/2 だから
A
t = T /2 =
−A
dx
v(x)
(4.12)
となる。右辺分母の v(x) はエネルギー保存則より求める。系の全エネルギー W は最大振幅 x = A のとき
のポテンシャルエネルギーに等しいので
Umax (A, β) = W =
1
1
αA2 + βA4
2
4
(4.13)
変位 x での運動エネルギー Ek (x, β) は全エネルギー W からその位置でのポテンシャルエネルギー U (x, β)
を差し引いたものだから
Ek (x, β) = W − U (x, β) =
1
1
1
1
αA2 + βA4 − αx2 + βx4
2
4
2
4
1
1
= α(A2 − x2 ) + β(A4 − x4 )
2
4
1
1
2
2
α + β(A2 + x2 )
= (A − x )
2
4
(4.14)
(4.15)
となる。Ek (x, β) = (1/2)mv 2 なので，これから
v(x) =
2Ek (x, β)/m
(4.16)
と求まる。これを（4.12）に入れると
A
T
=
2
dx
= m1/2
v(x)
−A
A
−A
1
2Ek (x, β)
dx
(4.17)
が得られる。これから周期 T と β の定性的な関係はわかる。
β > 0 の場合：β → 大 =⇒ T → 短 / A → 大 =⇒ T → 短
復元力が強くなると質点の速度が速くなり周期は短くなる。また，振幅が大きくなると質点の速度が
速くなり周期は短くなるということだね。
β < 0 の場合：| β |→ 大 =⇒ T → 長 /
A → 大 =⇒ T → 長
• コニー：非線形振動では周期 T を求めるだけでも大変なのね。ところで今の議論を整理すると，全エネル
ギー保存則を用いて
W =
.
..
1
m
2
dx
dt
x
dt = ±
2
+ U (x) −→
dx
=±
dt
dx
2{W − U (x)}/m
ということから周期 T を求めていこうというわけね。
10
2{W − U (x)}/m
(4.18)
• Ｋ氏：そうなんだ。§２でやったよね。結局（4.18）は
x
φ
dx
dt = ±
2{W − U (x)}
dφ
=±
2V (x)/m
(4.19)
となるんだったね。
さて，いよいよこの微分方程式を解いていこう。ここでは 3 つの方法＋αで解いてみることにする。
• コニー；その＋αとはなにかしら？
• Ｋ氏：気を持たすようないい方をしたけど，解析力学でいわゆる変分原理というのがあるだろう。作用積
分を極値するような運動が実現する運動というやつだね。この変分原理（Hamilton の原理）を使って解を
見つけようとする方法だね。ここでの話は解析力学からのアプローチはしていないから＋αとして紹介し
ようと思ったのだけど。
• コニー：ランダウの力学のテキストではラグランアンを使っていろいろ振動の問題を議論しているわね。
面白そうだから是非よろしくお願いするわ。
（１）エネルギー積分方式
• Ｋ氏：エネルギー積分方式といっているけどこのネーミングはあくまでここでの方言だからね。それはと
もかく，バネは x = −a, a (a > 0) の間を振動しているとすると
W − U (x) = (a2 − x2 )V (x)
(4.20)
と書ける。変位 x の最大，最小値は E − U (x) = 0 を満たす値だから，その値を | x |= a, (a > 0) とおく
と，振幅 a は次式を満たす。
1
1
W − αa2 − βa4 = 0,
2
4
... W =
1 2 1 4
αa + βa
2
4
(4.21)
次に，V (x) = Ax2 + Bx + C とおいて
1
1
1
1
W − U (x) = αa2 + βa4 − αx2 − βx4
2
4
2
4
= (a2 − x2 )(Ax2 + Bx + C)
(4.22)
この式は恒等式なので x の同次係数を等しいとおくことから係数 A, B, C を求めると
A=
1
1
1
β, B = 0, C = α + βa2
4
2
4
... V (x) =
1
{2α + β(x2 + a2 )}
4
(4.23)
を得る。これで V (x) が具体的に求められたので W − U (x) をキチンと書くと
1
1
1
W − U (x) = W − αx2 − βx4 = (2a2 α + a4 β − 2αx2 − βx4 )
2
4
4
1
= (a2 − x2 ){2α + β(a2 + x2 )}
4
となる。これをエネルギー積分方式の式に入れて
t=±
m
2
√
=± m
=±
x
dx
W − U (x)
ϕ
√
= ± 2m
x
dx
(a2
a sin ϕdϕ
a sin ϕ
m
α + βa2
ϕ
(α + βa2
1 − k 2 sin ϕ
dϕ
1 − k 2 sin ϕ
11
−
x2 ){2α
+ β(a2 + x2 )}
が得られる。右辺は§4.1「振幅の大きい振り子」のところででてきた楕円積分だね。したがって周期 T は
m
α + βa2
T =4
π/2
dϕ
1 − k 2 sin ϕ
0
=4
m
K(k)
α + βa2
(4.24)
となる。ここで母数 k は
k2 =
βa2
2(1 + a2 )
ちなみに β → 0 とすると（線形）単振動の周期 T0 を得る。
m
K(0) = 2π
α
T0 = 4
m
2π
=
,
α
ω0
( ω02 = α/m )
(4.25)
（２）変数変換方式（スケール変換）
• Ｋ氏：簡単化のために時間の尺度を変えて
α/m t を改めて t と書き， β/α x を改めて x と書く。
s=
α
t→t
m
y=
β/α x → x
そうすると dx/ds = x , d2 x/ds2 = x として
x(t) = x
x
¨(t) = −
m/α s ,
α
m
d2 y
d
=
ds2
ds
x˙ =
m/α x
m/α s ,
x
¨(t) =
m
x
α
m/α s
α
(y + y 3 ) · · （
· 4.10）より
β
dy dx
dx ds
β d2 x
=
α ds2
=
βm
x
¨=−
αα
βmα
ααm
α
(y + y 3 ) = −(y + y 3 )
β
となり，変数変換後の元の方程式のスタイルは
x
¨ = −(x + x3 )
(4.26)
と極めてサッパリした形に書ける。（4.26）の両辺に (dx/dt) をかけると
dx d2
dx
1 d
= −(x + x3 )
−→
2
dt dt
dt
2 dt
これを t で積分すると
1
2
dx
dt
2
2
dx
dt
=−
1 2 1 4
x + x
2
4
=W −
d
dt
1 2 1 4
x + x
2
4
= W − U (x)
(4.27)
を得る。W は積分定数で全エネルギーだ。右辺第 2 項はポテンシャルエネルギーを意味しているよね。
U (x) =
1 2 1 4
x + x
2
4
(4.28)
（4.27）を書き換えて
dx
= ± 2(W − U (x))
dt
1
t = ±√
2
x
dx
W − U (x)
(4.29)
これからの計算は１．のエネルギー積分方式でやったのと同じなので省略してもいいけど，ついでだから
フォローしておくと。。。
E − U (x) の中が正の範囲でおこなわれ，変位 x の最大，最小値は E −
U (x) = 0 を満たす値だね。その値を | x |= a, (a > 0) とおくと
質点の運動は上の被積分関数の
√
1
1
E − a2 − a4 = 0 −→ a2 = −1 + 1 + 4E
2
4
12
(4.30)
と求まる。したがって
1
1
E − U (x) = E − x2 − x4 = (a2 − x2 )V (x)
2
4
V (x) は x の 2 次式になるから V (x) = Ax2 + Bx + C とおいて x の同次数の係数比較より
1
1
1
a2 + 2
E − x2 − x4 = (a2 − x2 )(Ax2 + Bx + C) −→ A = , B = 0, C =
2
4
4
4
(4.31)
(4.32)
を得る。したがって（4.29）は x = a cos ϕ とおくと dx = −a sin ϕdϕ なので
√
t=± 2
x
dx
(a2
−
x2 )(2
+
a2
+
= ±√
x2 )
ϕ
1
1 + a2
dϕ
1 − k 2 sin2 ϕ
(4.33)
となる。
k2 =
a2
2(1 + a2 )
とした。k は母数と呼ばれるパラメーターだ。振動の周期 T はしたがって
√
π/2
4
4 2
dϕ
T = 4t = √
kK(k)
=
a
1 + a2 0
1 − k 2 sin2 ϕ
(4.34)
となる。T はここで使った単位で測った周期だね。元の単位に戻すと
1
1
1
W − U (x) = W − αx2 − βx4 = (2a2 α + a4 β − 2αx2 − βx4 )
2
4
4
1
= (a2 − x2 ){2α + β(a2 + x2 )}
4
.
. . t = ±
m
2
√
=± m
=±
x
dx
W − U (x)
ϕ
√
= ± 2m
x
dx
(a2
−
x2 ){2α
+ β(a2 + x2 )}
a sin ϕdϕ
(α + βa2
a sin ϕ
ϕ
m
α + βa2
1 − k 2 sin ϕ
dϕ
1 − k 2 sin ϕ
となる。したがって周期 T は
T =4
m
α + βa2
π/2
dφ
1−
0
k2
sin ϕ
m
K(k)
α + βa2
=4
(4.35)
と得られる。ここで母数 k は
k2 =
βa2
2(1 + a2 )
（３）摂動計算方式（逐次近似法）
• Ｋ氏：系の運動方程式 m¨
x = −(αx + βx3 ) を m = 1 として簡略化して書くと
x
¨ + ω 2 x + βx3 = 0,
ω=
√
と書けるね6 。この方程式をここでは採用することにしよう。
非線形項の βx3 は次の条件を満たす摂動として捉えるわけだね。
| ω 2 x | | βx3 |
初期条件は
t = 0 で x = a,
6
α/m → α, m/β → β と置き換えてもよい。
13
x˙ = 0
α
(4.36)
としておこう。摂動計算だから逐次近似で解を求めていくことになるんだが，ここでは 3 次の近似解まで
求めていこう。
◆摂動項を無視できる場合には，方程式は x
¨ + ω 2 x = 0 となり，その解は x = a cos ωt であるので，
（4.36）
の第 1 次の近似解として ω とわずかに異なる角振動数 p を持つとして
x = a cos pt
ω2
(ω 2 − p2 )
(4.37)
(4.38)
とおく。第 2 次の近似解は第 1 次の近似解をベースとしてより近似を高めていくというやり方で求めてい
くわけだね。あとの計算の都合上，ω を次のように表記しておく。
ω 2 = p2 + (ω 2 − p2 ),
(4.39)
これを使って（4.36）を書きなおすと
x
¨ + p2 x = −(ω 2 − p2 )x − βx3
(4.40)
と書ける。この式の右辺の項はいずれも微少量であるので，右辺をゼロとした場合の解（4.37）は（4.40）
の真の解にかなり近いということがいえる。そこでこの第 1 次の近似解を（4.40）に入れ，三角関数の 3
倍角公式7 を使って整理すると
x
¨ + p2 x = −(ω 2 − p2 )a cos ωt − βa3 cos3 pt
1
3
= − a(ω 2 − p2 ) + βa3 cos pt − βa3 cos 3pt
4
4
(4.41)
が得られる。この方程式を解いて第 2 近似の解が求まるわけだね。ここで注意すべきポイントとして，こ
の方程式の右辺第 1 項までを考えると
x
¨ + p2 x = f cos pt
(4.42)
と書け，外部から強制的に周期的な力を加えた強制振動の方程式の形をしている点だ。強制振動の方程式
は
x
¨ + ω02 = F cos ωt
で，この方程式の一般解は
x = a cos(ω0 t + α) + F
1
cos ωt
ω02 − ω 2
そして ω0 = ω の場合に共鳴（共振）が起こることが知られている。このことを頭において（4.42）をみ
ると，この方程式の一般解は
x = a cos(pt + α) + f
1
cos ωt
p 2 − p2
(4.43)
で右辺第 2 項に共鳴項が現れ，解が無限大の振幅を持つことになる。近似解を求める方程式（4.41）は右辺第
1 項により振幅無限大を起こす事態を含むことになるが，もともと解は周期解となることがわかっている
ので，このような解を含むのはまずい。そこで cos pt の項がゼロになるように p を選択し，この事態を避
けるようにする。つまり p を
3
3
a(ω 2 − p2 ) + βa3 = 0 −→ p2 = ω 2 + a2 β
4
4
(4.44)
とするわけだね。そうすると（4.41）は
1
x
¨ + p2 x = − βa3 cos 3pt
4
7
cos3θ = 4cos3 θ − 3 cos θ
14
(4.45)
となって 2 階線形微分方程式となり，一般解は“ 同次方程式の解＋特解 ”という形で表される8 。特解を
x = Acos3pt
(4.46)
とおいて，上の方程式に代入して係数 A を求めると
βa3
1
(−9p2 A + p2 A) cos 3pt = − βa3 cos 3pt −→ A =
4
32p2
(4.47)
となるね。同次方程式（x
¨ + p2 x = 0）の解は c1 ，c2 を定数として c1 cos pt + c2 sin pt となるので，一般解は
x = c1 cos pt + c2 sin pt + Acos3pt = c1 cos pt + c2 sin pt +
βa3
cos 3pt
32p2
と求まる。係数 c1 , c2 は初期条件から求めていくわけだけど，まず t で微分して
x˙ = −pc1 sin pt + pc2 cos pt −
初期条件より c1 , c2 を求めると
t = 0：x = a
t = 0：x˙ = 0
(4.48)
βa3
βa3
−→ c1 = a −
2
32p
32p2
0 = pc2 −→ c2 = 0
a = c1 +
よって一般解として
x=
3βa3
sin 3pt
32p
a−
βa3
32p2
cos pt +
βa3
cos 3pt
32p2
(4.49)
が得られる。これは第 2 次の近似解だね。振動の周期は T = 2π/p だね。ちなみに p2 は
p2 = ω 2 + (3/4)a2 β
(4.50)
・ コニー：非線形項 β > 0 なら周期は短くなるし，逆に β < 0 なら周期は長くなるということね。ところ
で摂動計算で共鳴現象のような項がでてきたけど，もともと単振動の系に非線形項が摂動として加わった
のだから振幅が発散するようなことはないわね。しかし，摂動計算ではそのような項がでてくる。これは
近似のやり方がまずいということだと思うけど，条件を設けてその困難を回避したわけね。
・Ｋ氏：そうなんだ。昔，天体力学での摂動計算でこのような事態に遭遇し，摂動計算はだめだと深刻な
議論になったそうだ。ちなみに，この擬似共鳴項を永年項（secular term）と呼んでいるね。
・コニー：そうなの，いろいろ歴史があるのね。ところで，2 次の近似解で非線形パラメータ β が顔をだ
したわね。また cos 3pt の高調波も混ざってきたわ。3 次の近似解は 2 次の近似解を元の方程式に入れて求
めていくのかしら？
◆高次の近似解
・ Ｋ氏：そうだね。基本的にそうなんだが，実は摂動展開の手法でもっと見通しのよいやり方をするんだ。
具体的にいうと，第 2 次の近似解である（4.50）と（4.49）をみると
ω 2 = p2 + C1 β, （C1 : 定数）
x = ϕ0 (t) + βϕ1 (t)
(4.51)
と先ほど コニー が指摘したように非線形パラメータ β の 1 次がかかっている形をしているだろう。そう
すると高次の近似解は
x = ϕ0 + βϕ1 + β 2 ϕ2 + β 3 ϕ4 · · ·
となるだろうと想像がつく。そこで
ω 2 = p2 + C1 β + C2 β 2 + · · ·
x = ϕ0 + βϕ1 + β 2 ϕ2 + · · ·
8
適当な微分方程式のテキストを参照ください。
15
(4.52)
という β のベキ級数を考え，これが（4.40）を満たすように C1 , C2 , · · · , p と関数 ϕ0 , ϕ1 , ϕ2 を決めること
ができれば，x の高次の近似解が得られることになるね。そこで試みに（4.52）を（4.40）に代入して β 3
の項までとると ( ←第 4 次の近似解まで考慮）
ϕ¨0 + β ϕ¨1 + β 2 ϕ¨2 + β 3 ϕ¨3 + p2 (ϕ0 + βϕ1 + β 2 ϕ2 + β 3 ϕ3 )
= ϕ¨0 + p2 ϕ0 + β(ϕ¨1 + ϕ1 ) + β 2 (ϕ¨2 + ϕ2 ) + β 3 (ϕ¨3 + ϕ3 )
(4.53)
= −{βC1 ϕ0 + β 2 (C2 ϕ0 + C1 ϕ1 ) + β 3 (C3 ϕ0 + C2 ϕ1 + C1 ϕ2 )}
(4.54)
−{βϕ30 + β 2 (3ϕ20 ϕ1 ) + β 3 (3ϕ20 ϕ2 + 3ϕ0 ϕ21 )}
(4.55)
= −{β(C1 ϕ0 + ϕ30 ) + β 2 (C2 ϕ0 + C1 ϕ1 + 3ϕ20 ϕ1 )
(4.56)
+β 3 (C3 ϕ0 + C2 ϕ1 + C1 ϕ2 + 3ϕ20 ϕ2 + 3ϕ0 ϕ21 )}
(4.57)
これが β の任意の値に対して成り立つためには β の同次数の係数が等しくなくてはならないから
ϕ¨0 + p2 ϕ0 = 0
ϕ¨1 + p2 ϕ1 = −C1 ϕ0 − ϕ30
ϕ¨2 + p2 ϕ2 = −C2 ϕ0 − C1 ϕ1 − 3ϕ20 ϕ1
ϕ¨3 + p2 ϕ3 = −C3 ϕ0 − C2 ϕ1 − C1 ϕ2 − 3ϕ20 ϕ2 − 3ϕ0 ϕ21
(4.58)
また，初期条件 t = 0：x = a, x˙ = 0 より
ϕ0 (0) + βϕ1 (0) + β 2 ϕ1 (0) + β 2 ϕ3 (0) = a
(ϕ˙ 0 )t=0 + β(ϕ˙ 1 )t=0 + β 2 (ϕ˙ 2 )t=0 + β 3 (ϕ˙ 3 )t=0 = 0
(4.59)
これも任意の β に対して成りたたなければならないので
ϕ0 (0) = a,
ϕ1 (0) = 0,
ϕ2 (0) = 0,
ϕ3 (0) = 0,
(ϕ˙ 0 )t=0
(ϕ˙ 1 )t=0
(ϕ˙ 2 )t=0
(ϕ˙ 3 )t=0
(4.60)
（4.58））の第 1 式より
ϕ¨0 + p2 ϕ0 = 0 −→ ϕ0 (t) = a cos pt + const
（4.60）の第 1 式の初期条件を入れると const = 0 となり
ϕ0 (t) = a cos pt
(4.61)
を得る。これは第 1 次の近似解だね。
次に，1 次近似解を（4.58）の第 2 式に入れて
3
1
ϕ¨1 + p2 ϕ1 = −C1 a cos pt − cos3 pt = − C1 a + a3 cos pt − a3 cos 3pt
4
4
を得るが，右辺に cos pt が顔をだしているのでこれをゼロになるように C1 を決めてやるんだね。そうす
ると
3
3
C1 a + a = 0 −→ C1 = − a2
4
4
が得られる。関数 ϕ1 は次の微分方程式を解いて得られる。
1
ϕ¨1 + p2 ϕ1 = − a3 cos 3pt
4
16
(4.62)
これは（4.45）と同じ形をしているので，同じような計算をやり，
（4.60）の第 2 式の初期条件から
ϕ1 =
a3
(cos 3pt − cospt)
32p2
(4.63)
を得る。これで ϕ1 が求まったから（4.52）に入れると
x = ϕ0 + βϕ1 = a cos pt +
= a−
βa3
32p2
cos pt +
βa3
(cos 3pt − cospt)
32p2
βa3
cos 3pt
32p2
(4.64)
(4.65)
3
p2 = ω 2 + a2 β
4
(4.66)
が得られる。これは第 2 次の近似解だね。
次に第 3 次の近似解だけど，今まで得られた結果を整理して書くと
ϕ0 = a cos pt,
ϕ1 =
a3
(cos 3pt − cospt),
32p2
3
C1 = − a2
4
(4.67)
だね。これを（4.58）の第 3 式に入れると
ϕ¨2 + p2 ϕ2 = −C2 a cos pt − −
= −C2 a cos pt +
3a2
4
a3
a3
(cos 3pt − cospt) − 3a2 cos2 pt
(cos 3pt − cospt)
2
32p
32p2
12a5
12a5
(cos3 pt − cos pt) −
cos2 pt(cos3 pt − cos pt)
2
128p
32p2
= − C2 a +
12a5
128p2
cos pt +
48a5
60a5
3
cos
pt
−
cos5 pt
128p2
128p2
= − C2 a +
12a5
128p2
cos pt +
3a5
(20 cos3 pt − 16 cos5 pt)
128p2
= − C2 a +
12a5
128p2
cos pt +
3a5
(5 cos pt − cos 5pt)
128p2
= − C2 a −
3a5
128p2
cos pt −
3a5
cos 5pt
128p2
ここで cos pt の項をゼロとおくと
C2 =
(4.68)
3a4
128p2
(4.69)
を得る。尚，上の式の展開では次の三角関数の 5 倍角の公式9 を使った。
cos 5θ = 16 cos5 θ − 20 cos3 θ + 5 cos θ
ということで（4.68）は次の線形微分方程式となる。
ϕ¨2 + p2 ϕ2 = −
3a5
cos 5pt
128p2
(4.70)
特解として x = A cos 5pt とし，これを（4.70）に入れて A を求めると
ϕ¨2 + p2 ϕ2 = −25Ap2 cos 5pt + Ap2 cos 5pt = −
... A =
9
a5
1024p2
3a5
cos 5pt
128p2
(4.71)
3 倍角までメジャーなので 4 倍角の公式をついでに載せておきます。
cos 4θ = 4 sin θ cos θ − 8 sin3 θ cos θ, sin 4θ = 4 sin θ cos θ − 8 sin 3θ cos θ, sin θ の 5 倍角は sin 5θ = 16 sin5 θ − 20 sin3 θ + 5 sin θ
17
したがって一般解は
ϕ2 (t) = C1 cospt + C2 sin pt +
となるね。ここで初期条件 t = 0：ϕ2 (0) = 0,
C1 = −
となるので，
（4.72）は
ϕ2 (t) = −
a5
cos 5pt
1024p4
(4.72)
ϕ˙ 2 (0) = 0 を入れて定数係数 C1 , C2 を求めると
a5
,
1024p4
C2 = 0
a5
a5
cos pt +
cos 5pt
4
1024p
1024p4
(4.73)
したがって第 3 次近似の解は，
x = ϕ0 + βϕ1 + β 2 ϕ2
= a−
βa3
β 2 a5
−
2
32p
1024p4
cos pt +
βa3
β 2 a5
cos 3pt +
cos 5pt
2
32p
1024p4
3
3a4 2
p2 = ω 2 − C1 β − C2 β 2 = ω 2 + a2 β −
β
4
128p2
と得られる。4 次以降の近似解も同様なやり方で求めていける。もしその気があればフォローしてみて。
√
ところで，p2 を求める式は p2 の 2 次式となっているので，p2 は p2 = A ± B という形になるね。求める
p2 は多い値の方をとるよ。1 次，2 次，3 次の近似解を描くと Fig.5 のようになる。2 次と 3 次の近似解で
はもう殆ど差がなくなっているね。ついでに厳密解と近似解から得られる周期の一覧表を載せておくよ。
F ig.5
ω = 10, a = 20, β = 0.01
x
“
x= a−
“
x= a−
”
β 2 a5
βa3
− 1024p
cos pt + 32p
4
2 cos 3pt +
”
3
3
βa
βa
cos pt + 32p
2 cos 3pt
32p2
βa3
32p2
β 2 a5
cos5pt
1024p4
x = a cos ωt
t
0
周期 T
厳密解
1 次近似解
2 次近似解
3 次近似解
(m = 1)
a = 20, α = 1, β = 0.01
βa2 ) K(k) ,
4 1/(α +
√
2π/ α, (α ∼ ω 2 )
2.81080
6.28319
2π/ α + 34 a2 β
√
2π 2/ α + 34 a2 β + (α + 32 a2 αβ +
3.14159
15a4 2 1/2
32 β )
3 次近似解の周期が 2 次より大きがどこかおかしい (要注意!!)
え∼っと，最後にお約束していた変分原理の方法だね。
18
3.18001
（４）＋α：変分原理方式
• Ｋ氏：
（4.36）の非線形方程式
x
¨ + ω02 x + βx3 = 0
を与えるラグランジアン L は
L=
(4.74)
1 2
β
(x˙ − ω02 x2 ) − x4
2
4
(4.75)
∂L
∂x
d
∂l
− dx
∂ x˙ = 0 より（4.74）がすぐでてくるから，このラグランジ
アンでよいことがわかるよね。非線形パラメータ β がゼロならば x = A cos ω0 t という解がある。いま，β
が小さいとして角振動数は ω0 からどの程度ずれるかを変分原理を使って調べてみよう。処方によって作
用積分 I の変分をとる。
だね。オイラー・ラグランジュの式
t2
I=
t2
Ldt
−→
δI = δ
t1
Ldt = 0
(4.76)
t1
これを満たす関数 x は実際の運動を記述するということだね。お試し関数として x = A cos ωt とおく。変
分 δI の両端時刻 t1 , t2 は運動の始点と終点の時刻で，その両端では力学変数の変化はゼロとする。お試し
関数として周期解をとったから，作用積分を 1 周期にわたってとれば，時刻 t1 と t2 の途中で振幅 A がい
かように変化しようとも両端での x の変分はゼロとなるので，両端固定の条件が満たされる。そこで積分
を 1 周期にわたって実行すると
2π/ω
I=
Ldt =
0
1 2π 2
3 π 4
A (ω − ω02 ) −
A β
2 ω
16 ω
(4.77)
が得られる。
F ig.6
A
x = A cos ωt
t1
ωt
t2
変分原理より ω 2 = ω02 + (3/4)A2 β という関係式が得られる！
次ぎに δI = 0 となる，I が極値をとる A の値を求めると
1π
∂I
=
∂A2
2ω
3
ω 2 − ω02 − A2 β
4
=0
−→
3
ω 2 = ω02 + A2 β
4
が得られ，振幅 A と角振動数 ω の関係は既に得ている摂動計算の 2 次の近似解（4.50）と同一だね！
• コニー：なるほど，ラグランジアンをうまく見つければ変分原理のやり方も便利ね。ただ，先ほどの摂動
展開では，高次の近似解が手間さえ惜しまなければ次々と得られていくわけだけど，変分原理ではこれ以
上の近似解に迫るのは難しいようね。
• Ｋ氏：できるのかも知れないけど，僕は知らない。ケースによって使い分ければいいんだね。
さて，今日はここらでお開きにしようか。
• コニー：そうね，ちょうどキリがいいみたいだから，大変お疲れ様でした。まだまだ聞きたいことがある
のだけど，それは次回ということで楽しみにしておくわ。
• Ｋ氏：そうだね。次回は「非線形振動その２」という感じでもう少し話を進めようか。例えば速度の 2 乗
に比例する抵抗が作用する場合とか云々ということで。
19
• コニー：楽しみだわ。（ふと柱時計を見て）アッ，そろそろ彼と飲みにいく時間になったので，今日は本当
にありがとうございました，これで失礼することにします。また次回よろしくお願いしますね。
• Ｋ氏：了解。気がはやると思うけど，外は随分暗くなってきたから気をつけてね。それじゃまた次回とい
うことで。
(了)
——————————————————————————–
またお会いできる機会を楽しみに，
，
，
G OOD L U C K !
S E E Y OU A G A I N !
20
© Copyright 2024 ExpyDoc
