曲面上の運動 — 拘束とラグランジュ乗数

ビーズは輪の上だけを動く — 拘束が自由度をどう削り、ラグランジュ乗数が拘束力の正体をどう教えてくれるか。

はじめに

第1章でニュートン力学を整理したとき、私たちは粒子が3次元空間を自由に漂うと仮定していた。しかし実際の問題のほとんどはそうではない。ビーズは針金の輪の上にしか乗っていないし、振り子のおもりは長さ一定の糸の先にぶら下がっており、コインは滑らずに転がる。この章ではこうした 拘束 (constraint) を2種類に分け、その中でも最もよく現れる ホロノミック (holonomic) 拘束を、球面 S2S^2 上の運動という絵を通して触ってみる。そして第7章で本格的に扱う ラグランジュ乗数 (Lagrange multiplier) が単なる数学トリックではなく、拘束力そのもの を測っているという事実を先取りする。

本論 1 — 二種類の拘束

拘束は形式によって二つに分かれる。

ホロノミック拘束 は座標だけの(場合によっては時間も含む)方程式

f(q,t)=0f(q, t) = 0

で書ける。ここで q=(q1,…,qn)q = (q^1, \dots, q^n) は粒子の座標の組である。この1本の式が配位空間 (configuration space, あり得る位置の集合) の次元をひとつ削る。たとえば半径 RR の円形の輪の上のビーズは、平面で f(x,y)=x2+y2−R2=0f(x, y) = x^2 + y^2 - R^2 = 0 を満たすから、2次元が1次元の曲線に落ちる。同じく単位半径の球面 S2S^2 上の点は x2+y2+z2=1x^2 + y^2 + z^2 = 1 を満たし、3次元が2次元の曲面に落ちる。

非ホロノミック拘束 は速度に課せられる線型条件

∑iai(q) q˙i=0\sum_i a_i(q)\, \dot q^i = 0

で書かれるが、この式は座標だけの式に積分できない。つまり f(q)=0f(q) = 0 の形に書き直すことができない。古典的な例は滑らずに転がるコインの転がり条件 — コインは平面上どこへでも行ける(座標空間は減らない)が、各瞬間の速度の向きは転がっている方向に縛られる。本章ではホロノミック拘束のみを深く扱う。力学の教科書で出会う拘束のほとんどがこちら側で、非ホロノミックは第7章以降の独立した節で扱う方がきれいだからだ。

本論 2 — 輪の上のビーズ、最初から最後まで

質量 mm のビーズが半径 RR の 鉛直 円形の輪の上を滑る。重力加速度を gg とする。ビーズの位置を、輪の最下点から測った角度 θ\theta (シータ) ひとつで指定しよう。すると

x=Rsin⁡θ,y=R(1−cos⁡θ).x = R\sin\theta, \qquad y = R(1 - \cos\theta).

運動エネルギーと位置エネルギーはそのまま落ちる:

T=12mR2θ˙2,V=mgR(1−cos⁡θ).T = \tfrac{1}{2} m R^2 \dot\theta^2, \qquad V = m g R (1 - \cos\theta).

第7章でラグランジアン L=T−VL = T - V を正式に導入するが、ここではその結果を先取りして使う。定数項 mgRmgR は運動方程式に影響しないので落とすと

L=12mR2θ˙2+mgRcos⁡θ.L = \tfrac{1}{2} m R^2 \dot\theta^2 + m g R \cos\theta.

これを オイラー–ラグランジュ方程式

ddt ⁣(∂L∂θ˙)−∂L∂θ=0\frac{d}{dt}\!\left( \frac{\partial L}{\partial \dot\theta} \right) - \frac{\partial L}{\partial \theta} = 0

に入れると mR2θ¨+mgRsin⁡θ=0m R^2 \ddot\theta + m g R \sin\theta = 0、すなわち

θ¨+gRsin⁡θ=0.\ddot\theta + \frac{g}{R} \sin\theta = 0.

これがまさに 振り子の方程式 である。もともとは (x,y,z)(x, y, z) の3座標に輪の拘束式2本がかかった 3-DOF の問題だったが、適切な一般化座標 θ\theta を一本選ぶことで一気に 1-DOF に落としている。これがラグランジュ形式が学生を魅了する最初の場面 — 拘束を 代数的に解いて消してしまう のだ。

本論 3 — ラグランジュ乗数が指し示すもの

座標を一般化座標に取り替えると、拘束力は方程式から姿を消す。しかし時にはその拘束力そのものが知りたい — たとえば輪がビーズに及ぼす 垂直抗力 がどれほどの大きさで、輪がその力に耐えられるか、というのは工学的な問いとして実在する。

そのときは座標を減らさず、元の座標 (x,y)(x, y) をそのまま残したまま、拘束式 f(x,y)=x2+y2−R2=0f(x, y) = x^2 + y^2 - R^2 = 0 を 別個に 課す。運動方程式に1項を足して

mr¨=−∇V+λ ∇fm \ddot{\mathbf{r}} = -\nabla V + \lambda\, \nabla f

と書く。ここで λ\lambda (ラムダ) が ラグランジュ乗数 である。∇f\nabla f は拘束曲面に垂直なベクトルだから、追加された λ∇f\lambda \nabla f は曲面に 垂直な力、すなわち今探していた垂直抗力の方向そのものを向く。その大きさが — 符号も含めて — λ\lambda である。

言い換えれば、ラグランジュ乗数は単に拘束式を満たすために動員した数学的道具ではなく、拘束力の振幅そのもの を指す物理量である。正式な導出は第7章「変分原理と拘束」で行う。ここでは「ラグランジュ乗数 = 拘束力の大きさ」という一行を頭に刻んでおけば十分だ。

Pythonで確かめる

# 輪の上のビーズ: theta_ddot = -(g/R) sin(theta) を RK4 で積分.
import numpy as np
import matplotlib.pyplot as plt

g, R = 9.81, 0.2
dt, T_end = 1e-3, 4.0
N = int(T_end / dt)

def f(state):
    th, om = state
    return np.array([om, -(g / R) * np.sin(th)])

state = np.array([np.pi / 3, 0.0])   # theta0 = 60度, 静止から出発
ts = np.linspace(0, T_end, N + 1)
ths = np.empty(N + 1)
oms = np.empty(N + 1)
ths[0], oms[0] = state

for i in range(N):
    k1 = f(state)
    k2 = f(state + 0.5 * dt * k1)
    k3 = f(state + 0.5 * dt * k2)
    k4 = f(state + dt * k3)
    state = state + (dt / 6) * (k1 + 2*k2 + 2*k3 + k4)
    ths[i + 1], oms[i + 1] = state

fig, ax = plt.subplots(2, 1, sharex=True)
ax[0].plot(ts, ths); ax[0].set_ylabel(r'$\theta$ [rad]')
ax[1].plot(ts, oms); ax[1].set_ylabel(r'$\dot\theta$ [rad/s]')
ax[1].set_xlabel('t [s]')
plt.tight_layout(); plt.show()

軌跡はちょうど一周期を閉じて同じ振幅で戻ってくる。振幅が小さければほぼ正弦波だが、θ0=π/3\theta_0 = \pi/3 ほどの大きな振幅では山の近くで少し平らになる 非調和性 (anharmonicity) が目に付く。sin⁡θ\sin\theta の非線型性が作り出す効果だ。

次章へ

第3章: テンソルと共変微分では、本章で曲面と呼んだもの — 自由度の減ったなめらかな空間 — を 多様体 という正式な名前で再び迎え、その上でベクトルをどう微分すれば座標の取り方に依存しないかを見る。ラグランジアンが「座標を自由に取り替えても同じ方程式を与える」という事実が偶然ではない、ということもそのとき明らかになる。