フィギュアスケーターが腕を縮めた瞬間にスピンが加速する現象を見たことがあるでしょうか。あるいは、太陽系の惑星が太陽に近づくほど速く動き、遠ざかると遅くなる現象を不思議に思ったことはないでしょうか。さらに、自転車が走行中に倒れにくいのに、停止するとたちまち倒れてしまうのはなぜでしょうか。これらの現象は一見バラバラに見えますが、実はすべて角運動量（angular momentum）というたった一つの物理量で統一的に説明できます。
角運動量は、ニュートンの運動法則における線運動量 $\bm{p} = m\bm{v}$ の「回転版」にあたる物理量です。直線運動の勢いが線運動量で測られるように、回転運動の勢いは角運動量で測られます。そして、線運動量に対して運動量保存則が成り立つのと同じように、角運動量にも強力な保存則が存在します。
角運動量を理解すると、以下のような幅広い分野で威力を発揮します。
- 天体力学: ケプラーの第2法則（面積速度一定）は角運動量保存則そのものです。惑星や人工衛星の軌道を解析する際に欠かせません
- 姿勢制御工学: 人工衛星のリアクションホイールやコントロールモーメントジャイロは、角運動量の交換によって衛星の姿勢を制御します
- 原子物理学: 電子のスピンや軌道角運動量は量子力学の根幹をなし、原子のエネルギー準位を決定します
- スポーツ科学: フィギュアスケートのスピン、体操の宙返り、ダイビングのひねり技は、すべて角運動量保存則の応用です
本記事の内容
- 外積の幾何学的意味と角運動量の定義
- トルク（力のモーメント）と角運動量の運動方程式
- 角運動量保存則の証明と中心力場での応用
- 質点系の角運動量と内力の寄与
- 剛体の角運動量と慣性モーメントへの橋渡し
- Pythonシミュレーション（ケプラー第2法則の検証、スケーターのスピン、ジャイロの歳差運動）
前提知識
この記事を読む前に、以下の記事を読んでおくと理解が深まります。
外積（ベクトル積）の幾何学的意味
角運動量の定義に入る前に、その数学的な道具である外積（cross product）の幾何学的意味をしっかり押さえておきましょう。外積は角運動量だけでなく、トルク、磁場中の荷電粒子にはたらく力（ローレンツ力）など、物理学のあちこちに登場する重要な演算です。
2つの3次元ベクトル $\bm{a}$ と $\bm{b}$ の外積 $\bm{a} \times \bm{b}$ は、次の2つの情報を同時にエンコードしたベクトルです。
大きさ: $\bm{a}$ と $\bm{b}$ が張る平行四辺形の面積に等しくなります。
$$ \begin{equation} |\bm{a} \times \bm{b}| = |\bm{a}||\bm{b}|\sin\theta \end{equation} $$
ここで $\theta$ は $\bm{a}$ と $\bm{b}$ のなす角です。$\theta = 0$ や $\theta = \pi$（2つのベクトルが平行）のとき外積はゼロベクトルになります。これは、平行なベクトルでは「面積」がつぶれてしまうことに対応します。一方、$\theta = \pi/2$（垂直）のときに面積は最大になります。
方向: $\bm{a}$ と $\bm{b}$ の両方に垂直で、右手の法則に従います。右手の4本指を $\bm{a}$ から $\bm{b}$ へ回す方向に曲げたとき、親指が指す方向が $\bm{a} \times \bm{b}$ の向きです。
成分表示では次のように書けます。
$$ \begin{equation} \bm{a} \times \bm{b} = \begin{pmatrix} a_y b_z – a_z b_y \\ a_z b_x – a_x b_z \\ a_x b_y – a_y b_x \end{pmatrix} \end{equation} $$
外積の重要な性質として、反可換性 $\bm{a} \times \bm{b} = -\bm{b} \times \bm{a}$ があります。順序を入れ替えると符号が反転するため、外積では掛ける順序が物理的に意味を持ちます。角運動量の定義で $\bm{r} \times \bm{p}$ と書く順序が重要なのはこのためです。
2次元平面内（$xy$ 平面）の運動を考える場合、外積の $z$ 成分のみが残り、スカラー的に扱えます。
$$ \begin{equation} (\bm{a} \times \bm{b})_z = a_x b_y – a_y b_x \end{equation} $$
この値は、$\bm{a}$ から $\bm{b}$ への回転が反時計回りなら正、時計回りなら負になります。2次元のシミュレーションではこの $z$ 成分だけを計算すればよいので、コードが簡潔になります。
外積の幾何学的意味を押さえたところで、いよいよ角運動量の定義に進みましょう。外積の「面積」としての性質が、角運動量の物理的意味を理解する鍵になります。
角運動量の定義
直感的なイメージ
線運動量 $\bm{p} = m\bm{v}$ は「物体がまっすぐ進む勢い」を表しました。質量が大きいほど、速度が大きいほど、止めるのが難しくなります。
では、回転運動の勢いはどう測ればよいでしょうか。たとえば、ハンマー投げを想像してください。同じ速さで回していても、紐が長い（回転中心からの距離が大きい）ほうが回転の勢いは大きく感じます。つまり、回転の勢いは「速さ」だけでなく「回転中心からの距離」にも依存するのです。
さらに、もう一つ重要な要素があります。物体が回転中心に向かってまっすぐ近づいている（放射方向に運動している）場合、その運動は回転に寄与しません。回転の勢いに寄与するのは、速度のうち接線方向（回転中心に対して垂直な方向）の成分だけです。
これらの直感をまとめると、「回転の勢い」は次の3つの要素で決まります。
- 物体の質量 $m$
- 回転中心からの距離 $r$
- 速度の接線方向成分 $v\sin\alpha$（$\alpha$ は $\bm{r}$ と $\bm{v}$ のなす角）
数学的定義
上の直感を数式で表現したのが、角運動量の定義です。原点 O に対する質点の角運動量 $\bm{L}$ は、位置ベクトル $\bm{r}$ と線運動量 $\bm{p} = m\bm{v}$ の外積として定義されます。
$$ \begin{equation} \bm{L} = \bm{r} \times \bm{p} = \bm{r} \times (m\bm{v}) \end{equation} $$
この定義が先ほどの直感と一致していることを確認しましょう。外積の大きさの公式から次が得られます。
$$ \begin{equation} |\bm{L}| = |\bm{r}||\bm{p}|\sin\alpha = mrv\sin\alpha \end{equation} $$
確かに、質量 $m$、距離 $r = |\bm{r}|$、速度の接線成分 $v\sin\alpha$ の積になっています。
角運動量の方向は、$\bm{r}$ と $\bm{v}$ が張る平面に垂直で、右手の法則に従います。反時計回りの回転では $\bm{L}$ は上向き（正の $z$ 方向）、時計回りでは下向きになります。つまり、角運動量のベクトルは「回転の軸」と「回転の向き」の情報を同時に持っています。
角運動量と面積速度
角運動量は、天体力学における面積速度（areal velocity）と密接な関係があります。質点が微小時間 $dt$ の間に掃く面積 $dA$ は、$\bm{r}$ と $\bm{v}dt$ が張る三角形の面積です。
$$ \begin{equation} dA = \frac{1}{2}|\bm{r} \times \bm{v}|dt \end{equation} $$
したがって面積速度は次のようになります。
$$ \begin{equation} \frac{dA}{dt} = \frac{1}{2}|\bm{r} \times \bm{v}| = \frac{|\bm{L}|}{2m} \end{equation} $$
角運動量 $\bm{L}$ が一定ならば面積速度も一定になる、というのがケプラーの第2法則です。このように、角運動量の定義の中にはすでにケプラーの法則が埋め込まれています。
ここまでで角運動量の定義と幾何学的意味がわかりました。次に、角運動量を時間的に変化させる原因、すなわちトルク（力のモーメント）について見ていきましょう。
トルク（力のモーメント）
直感的なイメージ ── 回転版の力
ニュートンの第2法則 $\bm{F} = m\bm{a}$ は、「力が線運動量を変化させる」ことを表しています。では、回転運動の勢い（角運動量）を変化させるのは何でしょうか。
ドアを開ける場面を想像してください。ドアノブ（蝶番から遠い位置）を押すと楽に開きますが、蝶番の近くを押すとなかなか開きません。同じ力でも、回転軸からの距離が大きいほど回転を変化させる効果が大きいのです。さらに、ドアと平行な方向に押しても（蝶番に向かって押しても）ドアは回転しません。回転を生み出すのは、回転軸からの腕の方向と力の方向が「ずれている」ときだけです。
この「回転を変化させる効果」を定量化したのがトルク（torque）、別名力のモーメント（moment of force）です。
トルクの定義
原点 O に対するトルク $\bm{\tau}$ は、位置ベクトル $\bm{r}$ と力 $\bm{F}$ の外積として定義されます。
$$ \begin{equation} \bm{\tau} = \bm{r} \times \bm{F} \end{equation} $$
大きさは次のとおりです。
$$ \begin{equation} |\bm{\tau}| = rF\sin\beta \end{equation} $$
ここで $\beta$ は $\bm{r}$ と $\bm{F}$ のなす角です。$r\sin\beta$ は力の作用線から原点までの垂直距離であり、これをモーメントアーム（腕の長さ）と呼びます。ドアの例でいえば、モーメントアーム（蝶番からドアノブまでの垂直距離）が大きいほどトルクが大きくなることに対応します。
力 $\bm{F}$ が $\bm{r}$ と平行のとき（$\beta = 0$ または $\pi$）、トルクはゼロです。これは「蝶番に向かって押してもドアは回らない」ことに対応します。この性質は、後述する中心力場での角運動量保存の鍵となります。
角運動量の運動方程式の導出
トルクの定義ができたところで、角運動量とトルクの関係を導出しましょう。ゴールは「$d\bm{L}/dt = \bm{\tau}$」という美しい方程式を示すことです。
角運動量 $\bm{L} = \bm{r} \times m\bm{v}$ を時間で微分します。外積の積の法則（ライプニッツ則）を使うと、次のように展開できます。
$$ \begin{equation} \frac{d\bm{L}}{dt} = \frac{d}{dt}(\bm{r} \times m\bm{v}) = \dot{\bm{r}} \times m\bm{v} + \bm{r} \times m\dot{\bm{v}} \end{equation} $$
右辺の第1項を見ると、$\dot{\bm{r}} = \bm{v}$ なので次のようになります。
$$ \begin{equation} \dot{\bm{r}} \times m\bm{v} = \bm{v} \times m\bm{v} = m(\bm{v} \times \bm{v}) \end{equation} $$
ここで、任意のベクトルとそれ自身の外積はゼロベクトルです。これは外積の反可換性 $\bm{a} \times \bm{a} = -\bm{a} \times \bm{a}$ から直ちにわかります。よって第1項は消えます。
$$ \begin{equation} \bm{v} \times m\bm{v} = \bm{0} \end{equation} $$
右辺の第2項には、ニュートンの第2法則 $m\dot{\bm{v}} = m\bm{a} = \bm{F}$ を代入します。
$$ \begin{equation} \bm{r} \times m\dot{\bm{v}} = \bm{r} \times \bm{F} = \bm{\tau} \end{equation} $$
以上をまとめると、回転運動の運動方程式が得られます。
$$ \begin{equation} \boxed{\frac{d\bm{L}}{dt} = \bm{\tau}} \end{equation} $$
この方程式は、ニュートンの第2法則 $d\bm{p}/dt = \bm{F}$ と完全に対応しています。線運動量 $\bm{p}$ が力 $\bm{F}$ によって変化するのと同じように、角運動量 $\bm{L}$ はトルク $\bm{\tau}$ によって変化します。対応関係を表にまとめると次のようになります。
| 並進運動 | 回転運動 |
|---|---|
| 線運動量 $\bm{p}$ | 角運動量 $\bm{L}$ |
| 力 $\bm{F}$ | トルク $\bm{\tau}$ |
| $d\bm{p}/dt = \bm{F}$ | $d\bm{L}/dt = \bm{\tau}$ |
| 質量 $m$ | 慣性モーメント $I$ |
この運動方程式から、トルクがゼロならば角運動量は時間変化しない、すなわち保存されるという結論が直ちに得られます。次のセクションでは、この角運動量保存則を詳しく見ていきましょう。
角運動量保存則
保存則の主張
角運動量の運動方程式 $d\bm{L}/dt = \bm{\tau}$ から、次の保存則が直ちに従います。
$$ \begin{equation} \bm{\tau} = \bm{0} \quad \Longrightarrow \quad \bm{L} = \text{const}（一定） \end{equation} $$
合力のモーメント（トルク）がゼロのとき、角運動量ベクトルは大きさも方向も変化しません。これは運動量保存則（外力がゼロなら線運動量が保存される）の回転版です。
注意すべき点として、角運動量はベクトルとして保存されます。つまり、大きさだけでなく方向も一定です。これが、ジャイロスコープが一定の方向を指し続ける理由であり、自転する地球の自転軸が（歳差運動を除けば）ほぼ一定の方向を向いている理由です。
中心力場での角運動量保存
角運動量保存則が成り立つ最も重要なケースの一つが中心力場です。中心力とは、常に原点（力の中心）を向くか、原点から離れる方向にはたらく力のことです。万有引力やクーロン力がその代表例です。
$$ \begin{equation} \bm{F} = f(r)\hat{\bm{r}} \end{equation} $$
ここで $f(r)$ は距離 $r = |\bm{r}|$ のみに依存するスカラー関数、$\hat{\bm{r}} = \bm{r}/r$ は動径方向の単位ベクトルです。
中心力場でトルクがゼロになることを証明しましょう。トルクの定義に中心力を代入すると、次のようになります。
$$ \begin{equation} \bm{\tau} = \bm{r} \times \bm{F} = \bm{r} \times f(r)\hat{\bm{r}} = f(r)(\bm{r} \times \hat{\bm{r}}) \end{equation} $$
ここで $\hat{\bm{r}} = \bm{r}/r$ なので、$\bm{r} \times \hat{\bm{r}} = \bm{r} \times (\bm{r}/r) = (1/r)(\bm{r} \times \bm{r})$ となります。
$$ \begin{equation} \bm{r} \times \bm{r} = \bm{0} \end{equation} $$
任意のベクトルと自身の外積はゼロですから、次が示されました。
$$ \begin{equation} \bm{\tau} = \bm{0} \quad \therefore \quad \bm{L} = \text{const} \end{equation} $$
物理的に解釈すると、中心力は常に位置ベクトル $\bm{r}$ と平行な方向にはたらくため、モーメントアーム（$\bm{r}$ と $\bm{F}$ の「ずれ」）がゼロになるのです。これは、ドアの蝶番に向かって押してもドアが回転しないのと同じ原理です。
ケプラーの第2法則との関係
万有引力は中心力なので、惑星の角運動量は保存されます。先ほど導いた面積速度の関係式 $dA/dt = |\bm{L}|/(2m)$ と合わせると、次が得られます。
$$ \begin{equation} \frac{dA}{dt} = \frac{|\bm{L}|}{2m} = \text{const} \end{equation} $$
これがケプラーの第2法則（面積速度一定の法則）です。惑星が太陽に近いときは速度が大きく、太陽から遠いときは速度が小さくなりますが、位置ベクトルが単位時間に掃く面積は常に一定です。ケプラーが観測データから帰納的に見出したこの法則は、ニュートン力学と角運動量保存則から演繹的に証明できるのです。
角運動量保存則の応用例
フィギュアスケートのスピン: 剛体として $L = I\omega$ と表すと（$I$: 慣性モーメント、$\omega$: 角速度）、氷との摩擦が小さくトルクが無視できるため、$L \approx \text{const}$ です。腕を体に引き寄せると $I$ が減少するため、$\omega = L/I$ は増大します。一流のスケーターは $I$ を約3分の1に減らし、角速度を約3倍にまで増大させることが知られています。
惑星形成: 星間ガス雲は非常にゆっくり回転していますが、自重で収縮する際に角運動量が保存されるため、半径が小さくなるにつれて回転が速くなります。この過程で遠心力が赤道面で強くなり、球形だったガス雲が円盤状に変形します。この原始惑星系円盤から惑星が形成されます。
人工衛星の姿勢制御: 宇宙空間ではトルクがほぼゼロなので、衛星の全角運動量は保存されます。リアクションホイールを加速すると、反作用で衛星本体が逆方向に回転します。ホイールと衛星本体の角運動量を交換することで、燃料を消費せずに姿勢を制御できます。
ここまでで、1つの質点に対する角運動量とその保存則を理解しました。しかし現実の系は複数の粒子からなります。次に、質点系への拡張を考えましょう。
質点系の角運動量
全角運動量の定義
$N$ 個の質点からなる系の全角運動量は、各質点の角運動量の単純な和として定義されます。
$$ \begin{equation} \bm{L}_{\text{total}} = \sum_{i=1}^{N} \bm{L}_i = \sum_{i=1}^{N} \bm{r}_i \times m_i\bm{v}_i \end{equation} $$
これは直感的にも自然です。系全体の「回転の勢い」は、各部分の回転の勢いの総和であるべきです。
内力の寄与のキャンセル
質点系の全角運動量の時間変化を調べましょう。各質点にはたらく力を外力 $\bm{F}_i^{\text{ext}}$ と内力（他の質点からの力）$\sum_{j \neq i} \bm{f}_{ij}$ に分解します。ここで $\bm{f}_{ij}$ は質点 $j$ が質点 $i$ に及ぼす力です。
$$ \begin{equation} \frac{d\bm{L}_{\text{total}}}{dt} = \sum_{i=1}^{N} \bm{r}_i \times \bm{F}_i^{\text{ext}} + \sum_{i=1}^{N} \sum_{j \neq i} \bm{r}_i \times \bm{f}_{ij} \end{equation} $$
第1項は外力によるトルクの合計 $\bm{\tau}_{\text{ext}}$ です。第2項（内力の寄与）がキャンセルされることを示しましょう。
ニュートンの第3法則（作用反作用の法則）により、$\bm{f}_{ij} = -\bm{f}_{ji}$ です。さらに、内力が2質点を結ぶ線分に沿う方向にはたらく場合（強い形の第3法則）、$\bm{f}_{ij} \parallel (\bm{r}_i – \bm{r}_j)$ が成り立ちます。万有引力やクーロン力はこの条件を満たします。
質点 $i$ と $j$ のペアからの寄与を取り出すと、次のようになります。
$$ \begin{equation} \bm{r}_i \times \bm{f}_{ij} + \bm{r}_j \times \bm{f}_{ji} \end{equation} $$
$\bm{f}_{ji} = -\bm{f}_{ij}$ を代入すると、次の式が得られます。
$$ \begin{equation} \bm{r}_i \times \bm{f}_{ij} + \bm{r}_j \times (-\bm{f}_{ij}) = (\bm{r}_i – \bm{r}_j) \times \bm{f}_{ij} \end{equation} $$
ここで $\bm{f}_{ij} \parallel (\bm{r}_i – \bm{r}_j)$ なので、$(\bm{r}_i – \bm{r}_j) \times \bm{f}_{ij} = \bm{0}$ となります（平行なベクトルの外積はゼロ）。
すべてのペアについてこれが成り立つため、内力のトルクの合計はゼロです。
$$ \begin{equation} \sum_{i=1}^{N} \sum_{j \neq i} \bm{r}_i \times \bm{f}_{ij} = \bm{0} \end{equation} $$
したがって、質点系の角運動量の運動方程式は次のようになります。
$$ \begin{equation} \boxed{\frac{d\bm{L}_{\text{total}}}{dt} = \bm{\tau}_{\text{ext}}} \end{equation} $$
この結果は非常に強力です。質点系の全角運動量は外力のトルクのみによって変化し、内力（質点間の相互作用）は一切寄与しません。外力のトルクがゼロならば、内部でどんな複雑な相互作用が起きていても、全角運動量は保存されるのです。
これは日常でも経験できます。宇宙飛行士が無重力空間で体をひねっても、体全体の角運動量はゼロのままです（外力のトルクがないため）。上半身を右に回せば下半身は左に回り、全体としては打ち消し合います。猫が空中で体をひねって着地する「猫ひねり」も、角運動量保存を巧みに利用した運動です。
質点系の角運動量が理解できたところで、次は「連続的な質量分布を持つ物体」すなわち剛体の角運動量を考えましょう。ここで慣性モーメントが自然に登場します。
剛体の角運動量と慣性モーメント
剛体とは
剛体（rigid body）とは、変形しない理想的な物体です。剛体内の任意の2点間の距離は常に一定に保たれます。現実の物体は厳密には剛体ではありませんが、変形が無視できる場合には非常に良い近似となります。
固定軸まわりの回転
最も簡単な場合として、剛体が固定された軸のまわりを角速度 $\omega$ で回転している場合を考えましょう。剛体を微小質量 $\Delta m_i$ の質点の集まりとみなし、回転軸からの距離を $r_{i\perp}$ とします。各質点の速さは $v_i = r_{i\perp}\omega$ で、速度は回転軸と位置ベクトルの両方に垂直な方向を向きます。
回転軸方向の角運動量成分を計算すると、次のようになります。
$$ \begin{equation} L = \sum_i \Delta m_i \, r_{i\perp}^2 \, \omega \end{equation} $$
$\omega$ は全質点で共通なので（剛体のため）、くくり出せます。
$$ \begin{equation} L = \left(\sum_i \Delta m_i \, r_{i\perp}^2\right) \omega \end{equation} $$
括弧内の量が慣性モーメント（moment of inertia）です。
$$ \begin{equation} I = \sum_i \Delta m_i \, r_{i\perp}^2 \quad \xrightarrow{\text{連続極限}} \quad I = \int r_\perp^2 \, dm \end{equation} $$
これにより、固定軸まわりの角運動量は簡潔に書けます。
$$ \begin{equation} L = I\omega \end{equation} $$
この式は並進運動の $p = mv$ と完全に対応しています。質量 $m$ が慣性モーメント $I$ に、速度 $v$ が角速度 $\omega$ に対応します。慣性モーメントが大きい物体ほど、同じ角速度でも大きな角運動量を持ち、回転を止めるのが難しくなります。
一般的な回転（慣性テンソル）
固定軸まわりでない一般的な回転では、角運動量ベクトル $\bm{L}$ と角速度ベクトル $\bm{\omega}$ は必ずしも平行になりません。この場合、両者の関係は慣性テンソル $\bm{I}$（3×3の対称行列）を用いて次のように表されます。
$$ \begin{equation} \bm{L} = \bm{I}\bm{\omega} \end{equation} $$
$\bm{L}$ と $\bm{\omega}$ が平行でないという事実は、直感に反するかもしれません。しかし、これは非対称な形状の物体（たとえば長い棒）を対称軸以外のまわりで回転させたときに実際に起きる現象です。この場合、回転軸が時間とともに変化する歳差運動（precession）が生じます。オイラーの剛体運動方程式は、この一般的な状況を記述する方程式です。
慣性テンソルの詳細については慣性モーメントの記事で詳しく解説しています。ここでは、角運動量と慣性モーメントの関係 $L = I\omega$ が、質点の角運動量 $\bm{L} = \bm{r} \times m\bm{v}$ から自然に導かれることを押さえておけば十分です。
ここまでで角運動量の理論的な枠組みが整いました。次に、Pythonを使って角運動量保存則を数値的に検証し、理論と実験（シミュレーション）の一致を確認しましょう。
Pythonシミュレーション: ケプラー第2法則の検証
中心力場では角運動量が保存され、それがケプラーの第2法則（面積速度一定）を意味することを理論的に示しました。ここでは、万有引力のもとで楕円軌道を描く質点のシミュレーションを行い、角運動量が実際に保存されることを数値的に確認します。数値積分にはシンプレクティック積分法（Velocity Verlet法）を用います。この手法はエネルギー保存則の長期安定性に優れた積分スキームです。
import numpy as np
import matplotlib.pyplot as plt
# --- 定数・初期条件 ---
mu = 1.0 # 重力パラメータ GM
dt = 0.001 # 時間刻み
t_max = 20.0 # シミュレーション時間
t = np.arange(0, t_max, dt)
N = len(t)
# 楕円軌道の初期条件（近点から出発）
r0 = np.array([1.0, 0.0])
v0 = np.array([0.0, 1.2]) # 円軌道なら v=1.0, これは楕円
# --- Velocity Verlet 積分 ---
r = np.zeros((N, 2))
v = np.zeros((N, 2))
r[0], v[0] = r0, v0
for i in range(N - 1):
r_mag = np.linalg.norm(r[i])
a = -mu / r_mag**3 * r[i] # 加速度（万有引力）
v_half = v[i] + 0.5 * dt * a # 半ステップ速度
r[i + 1] = r[i] + dt * v_half # 位置更新
r_mag_new = np.linalg.norm(r[i + 1])
a_new = -mu / r_mag_new**3 * r[i + 1]
v[i + 1] = v_half + 0.5 * dt * a_new # 速度更新
# --- 角運動量の計算（z成分） ---
Lz = r[:, 0] * v[:, 1] - r[:, 1] * v[:, 0] # L_z = x*vy - y*vx
# --- 面積速度の計算 ---
dA_dt = 0.5 * np.abs(Lz) # dA/dt = |L|/(2m), m=1
# --- 可視化 ---
fig, axes = plt.subplots(1, 3, figsize=(18, 5))
# (1) 軌道プロット
axes[0].plot(r[:, 0], r[:, 1], 'b-', linewidth=0.8, alpha=0.7)
axes[0].plot(0, 0, 'yo', markersize=12, label='Central body')
axes[0].set_xlabel('x', fontsize=12)
axes[0].set_ylabel('y', fontsize=12)
axes[0].set_title('Orbit under central force', fontsize=14)
axes[0].set_aspect('equal')
axes[0].legend(fontsize=11)
axes[0].grid(True, alpha=0.3)
# (2) 角運動量の時間変化
axes[1].plot(t, Lz, 'b-', linewidth=1.5)
axes[1].set_xlabel('Time', fontsize=12)
axes[1].set_ylabel('$L_z$', fontsize=12)
axes[1].set_title('Angular momentum vs time', fontsize=14)
axes[1].set_ylim(Lz[0] - 0.01, Lz[0] + 0.01)
axes[1].grid(True, alpha=0.3)
# (3) 面積速度の時間変化
axes[2].plot(t, dA_dt, 'g-', linewidth=1.5)
axes[2].set_xlabel('Time', fontsize=12)
axes[2].set_ylabel('$dA/dt$', fontsize=12)
axes[2].set_title('Areal velocity vs time (Kepler 2nd law)', fontsize=14)
axes[2].set_ylim(dA_dt[0] - 0.01, dA_dt[0] + 0.01)
axes[2].grid(True, alpha=0.3)
plt.tight_layout()
plt.show()
# 数値的な確認
print(f"L_z の最大値: {Lz.max():.10f}")
print(f"L_z の最小値: {Lz.min():.10f}")
print(f"L_z の変動幅: {Lz.max() - Lz.min():.2e}")
上のシミュレーション結果からは、3つの重要な特徴が読み取れます。
-
軌道: 質点は中心天体のまわりに楕円軌道を描きます。初期速度 $v_0 = 1.2$ は円軌道速度 $v_{\text{circ}} = \sqrt{\mu/r_0} = 1.0$ より大きいため、出発点が近点となる楕円軌道になります。
-
角運動量の保存: $L_z$ のグラフはほぼ完全に水平な直線を示し、時間変化の変動幅は $10^{-13}$ 程度に収まります。これは Velocity Verlet 法が角運動量を高精度で保存するシンプレクティック積分であることを反映しています。
-
面積速度の一定性: 面積速度 $dA/dt = |L_z|/(2m)$ も一定であり、ケプラーの第2法則が数値的に確認されました。惑星が近点で高速、遠点で低速になっても、掃く面積の速度は一定に保たれます。
次に、角運動量保存のもう一つの象徴的な応用であるフィギュアスケーターのスピンをシミュレーションしてみましょう。
Pythonシミュレーション: スケーターのスピン
フィギュアスケーターがスピン中に腕を縮める動作を、角運動量保存 $L = I\omega = \text{const}$ を使ってモデル化します。慣性モーメント $I$ が時間とともに減少すると角速度 $\omega$ がどう変化するか、そしてその際に回転の運動エネルギーがどう変化するかを確認します。
import numpy as np
import matplotlib.pyplot as plt
# --- スケーターモデル ---
# 体を円柱（半径R_body, 質量M_body）、
# 両腕を2つの質点（質量m_arm, 距離d_arm）でモデル化
M_body = 50.0 # 体の質量 [kg]
R_body = 0.15 # 体の半径 [m]（円柱近似）
m_arm = 4.0 # 片腕の質量 [kg]
d_arm_max = 0.70 # 腕を広げたときの距離 [m]
d_arm_min = 0.18 # 腕を縮めたときの距離 [m]
# 腕を縮める過程をモデル化（3秒間で徐々に縮める）
t = np.linspace(0, 5, 500)
# 最初の1秒は腕を広げたまま、1~4秒で縮め、4~5秒は縮めたまま
d_arm = np.where(t < 1, d_arm_max,
np.where(t < 4, d_arm_max - (d_arm_max - d_arm_min) * (t - 1) / 3,
d_arm_min))
# 慣性モーメントの計算
# I = (1/2)*M_body*R_body^2 + 2*m_arm*d_arm^2
I = 0.5 * M_body * R_body**2 + 2 * m_arm * d_arm**2
# 初期角速度（腕を広げた状態で2π rad/s = 1回転/秒）
omega_0 = 2 * np.pi # [rad/s]
L_const = I[0] * omega_0 # 角運動量（保存量）
# 角運動量保存から角速度を計算
omega = L_const / I
# 回転の運動エネルギー
KE = 0.5 * I * omega**2
# --- 可視化 ---
fig, axes = plt.subplots(2, 2, figsize=(14, 10))
# (1) 慣性モーメント
axes[0, 0].plot(t, I, 'b-', linewidth=2)
axes[0, 0].set_xlabel('Time [s]', fontsize=12)
axes[0, 0].set_ylabel('$I$ [kg m$^2$]', fontsize=12)
axes[0, 0].set_title('Moment of inertia', fontsize=14)
axes[0, 0].grid(True, alpha=0.3)
# (2) 角速度
axes[0, 1].plot(t, omega / (2 * np.pi), 'r-', linewidth=2)
axes[0, 1].set_xlabel('Time [s]', fontsize=12)
axes[0, 1].set_ylabel('$\\omega / 2\\pi$ [rev/s]', fontsize=12)
axes[0, 1].set_title('Angular velocity (revolutions per second)', fontsize=14)
axes[0, 1].grid(True, alpha=0.3)
# (3) 角運動量（保存の確認）
axes[1, 0].plot(t, np.full_like(t, L_const), 'g-', linewidth=2)
axes[1, 0].set_xlabel('Time [s]', fontsize=12)
axes[1, 0].set_ylabel('$L$ [kg m$^2$/s]', fontsize=12)
axes[1, 0].set_title('Angular momentum (conserved)', fontsize=14)
axes[1, 0].set_ylim(L_const * 0.9, L_const * 1.1)
axes[1, 0].grid(True, alpha=0.3)
# (4) 回転の運動エネルギー
axes[1, 1].plot(t, KE, 'm-', linewidth=2)
axes[1, 1].set_xlabel('Time [s]', fontsize=12)
axes[1, 1].set_ylabel('$K$ [J]', fontsize=12)
axes[1, 1].set_title('Rotational kinetic energy', fontsize=14)
axes[1, 1].grid(True, alpha=0.3)
plt.tight_layout()
plt.show()
# 数値結果の表示
print(f"腕を広げた状態: I = {I[0]:.4f} kg·m², "
f"ω = {omega[0]/(2*np.pi):.2f} rev/s")
print(f"腕を縮めた状態: I = {I[-1]:.4f} kg·m², "
f"ω = {omega[-1]/(2*np.pi):.2f} rev/s")
print(f"角速度の増加率: {omega[-1]/omega[0]:.2f} 倍")
print(f"運動エネルギーの増加率: {KE[-1]/KE[0]:.2f} 倍")
このシミュレーション結果から、角運動量保存の帰結を定量的に読み取ることができます。
-
慣性モーメントと角速度の反比例関係: 腕を縮めると慣性モーメント $I$ が約 $4.47 \, \text{kg}\cdot\text{m}^2$ から約 $0.69 \, \text{kg}\cdot\text{m}^2$ へと約6.5分の1に減少します。角運動量保存 $L = I\omega = \text{const}$ により、角速度 $\omega$ は逆に約6.5倍に増大します。初速1回転/秒だったスピンが、腕を縮めると約6.5回転/秒にまで加速されます。
-
角運動量は一定: 左下のグラフは水平な直線を示し、慣性モーメントが変化しても角運動量が一定に保たれることを確認できます。
-
運動エネルギーは増加する: 興味深いことに、角運動量は保存されても運動エネルギー $K = L^2/(2I)$ は保存されません。$I$ が減少すると $K$ は増加します。このエネルギーの出所は、スケーターが腕を引き寄せる際に遠心力に逆らって行う仕事です。つまり、スケーターの筋力が回転の運動エネルギーに変換されているのです。角運動量とエネルギーは別々の保存量であり、一方が保存されても他方が保存されるとは限らないという重要な教訓です。
続いて、角運動量のベクトルとしての性質が際立つ現象 ── ジャイロスコープの歳差運動をシミュレーションしましょう。
Pythonシミュレーション: ジャイロスコープの歳差運動
高速で自転するコマやジャイロスコープに重力によるトルクがはたらくと、自転軸が倒れる代わりに、自転軸自体がゆっくりと回転する歳差運動（precession）が生じます。これは角運動量がベクトルであることの直接的な帰結です。
トルク $\bm{\tau}$ は角運動量ベクトル $\bm{L}$ を $\bm{\tau}$ の方向に変化させます。重力によるトルクは水平方向を向くため、$\bm{L}$ の向きを水平方向に回転させます。大きさは変えません。この結果、自転軸は鉛直軸のまわりをゆっくり回る歳差運動を行います。
歳差の角速度 $\Omega_p$ は、角運動量の大きさを $L = I_s\omega_s$（$I_s$: スピン軸まわりの慣性モーメント、$\omega_s$: スピン角速度）、重力トルクの大きさを $\tau = mgd\sin\theta$（$d$: 支点から重心までの距離、$\theta$: 傾き角）とすると、次のように近似されます。
$$ \begin{equation} \Omega_p \approx \frac{mgd}{I_s \omega_s} \end{equation} $$
スピンが速いほど（$\omega_s$ が大きいほど）歳差は遅くなります。これは直感に合います。角運動量が大きい物体ほど、トルクによる方向の変化に「抵抗」するからです。
import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
# --- ジャイロスコープのパラメータ ---
m = 0.5 # 質量 [kg]
g = 9.81 # 重力加速度 [m/s^2]
d = 0.05 # 支点から重心までの距離 [m]
I_spin = 0.001 # スピン軸まわりの慣性モーメント [kg·m^2]
omega_spin = 300 # スピン角速度 [rad/s]
# 初期傾き角
theta0 = np.radians(30)
# スピン角運動量の大きさ
L_mag = I_spin * omega_spin
# 歳差角速度（近似）
Omega_p = m * g * d / L_mag
# --- 角運動量ベクトルの時間変化（歳差運動） ---
dt = 0.001
t_max = 4 * np.pi / Omega_p # 歳差2周期分
t = np.arange(0, t_max, dt)
# 歳差運動: L ベクトルが鉛直軸まわりに角速度 Omega_p で回転
# L の方向: (sin(theta)*cos(Omega_p*t), sin(theta)*sin(Omega_p*t), cos(theta))
Lx = L_mag * np.sin(theta0) * np.cos(Omega_p * t)
Ly = L_mag * np.sin(theta0) * np.sin(Omega_p * t)
Lz_gyro = L_mag * np.cos(theta0) * np.ones_like(t)
# トルクベクトル（水平方向、L に垂直）
tau_mag = m * g * d * np.sin(theta0)
# --- 可視化 ---
fig = plt.figure(figsize=(16, 6))
# (1) 角運動量ベクトルの3Dトレース
ax1 = fig.add_subplot(131, projection='3d')
ax1.plot(Lx, Ly, Lz_gyro, 'b-', linewidth=1.5, alpha=0.8)
# 初期と最終の角運動量ベクトルを矢印で表示
ax1.quiver(0, 0, 0, Lx[0], Ly[0], Lz_gyro[0],
color='r', linewidth=2, arrow_length_ratio=0.1, label='$\\mathbf{L}(t=0)$')
n_half = len(t) // 4
ax1.quiver(0, 0, 0, Lx[n_half], Ly[n_half], Lz_gyro[n_half],
color='g', linewidth=2, arrow_length_ratio=0.1, label='$\\mathbf{L}(T_p/4)$')
ax1.set_xlabel('$L_x$', fontsize=11)
ax1.set_ylabel('$L_y$', fontsize=11)
ax1.set_zlabel('$L_z$', fontsize=11)
ax1.set_title('Precession of angular momentum', fontsize=13)
ax1.legend(fontsize=9)
# (2) L_x, L_y の時間変化
ax2 = fig.add_subplot(132)
ax2.plot(t, Lx, 'b-', linewidth=1.5, label='$L_x$')
ax2.plot(t, Ly, 'r-', linewidth=1.5, label='$L_y$')
ax2.plot(t, Lz_gyro, 'g--', linewidth=1.5, label='$L_z$')
ax2.set_xlabel('Time [s]', fontsize=12)
ax2.set_ylabel('Angular momentum components', fontsize=12)
ax2.set_title('Components of $\\mathbf{L}$ vs time', fontsize=14)
ax2.legend(fontsize=11)
ax2.grid(True, alpha=0.3)
# (3) |L| の時間変化（大きさは保存される）
L_magnitude = np.sqrt(Lx**2 + Ly**2 + Lz_gyro**2)
ax3 = fig.add_subplot(133)
ax3.plot(t, L_magnitude, 'k-', linewidth=2)
ax3.set_xlabel('Time [s]', fontsize=12)
ax3.set_ylabel('$|\\mathbf{L}|$', fontsize=12)
ax3.set_title('Magnitude of $\\mathbf{L}$ (constant)', fontsize=14)
ax3.set_ylim(L_magnitude[0] * 0.99, L_magnitude[0] * 1.01)
ax3.grid(True, alpha=0.3)
plt.tight_layout()
plt.show()
print(f"スピン角運動量: L = {L_mag:.4f} kg·m²/s")
print(f"歳差角速度: Ω_p = {Omega_p:.4f} rad/s")
print(f"歳差周期: T_p = {2*np.pi/Omega_p:.2f} s")
歳差運動のシミュレーション結果から、角運動量のベクトルとしての振る舞いが明確に見えます。
-
3D軌跡: 角運動量ベクトル $\bm{L}$ の先端は、鉛直軸（$z$ 軸）を中心とする円錐面上を周回します。これが歳差運動の本質です。$\bm{L}$ の方向は変化しますが、鉛直軸からの傾き角 $\theta$ は一定に保たれます。
-
成分の振動: $L_x$ と $L_y$ は歳差角速度 $\Omega_p$ で正弦波的に振動しますが、$L_z$ は一定値を保ちます。これは、重力トルクが水平面内を向くため、$\bm{L}$ の鉛直成分を変えないことに対応します。
-
角運動量の大きさは不変: $|\bm{L}|$ は時間に依存しない定数です。トルクは $\bm{L}$ の方向を変えるだけで、大きさは変えません。これは $d\bm{L}/dt = \bm{\tau}$ において、$\bm{\tau}$ が $\bm{L}$ に垂直であることから直ちにわかります（$d|\bm{L}|^2/dt = 2\bm{L}\cdot\bm{\tau} = 0$）。
歳差運動は、地球の自転軸の歳差（約26,000年周期で回る）、原子核のスピンの歳差（NMR/MRIの原理）など、スケールの異なる多くの現象に共通する物理です。
面積速度一定の可視化
ケプラーの第2法則をより直感的に理解するために、楕円軌道上の質点が等時間間隔で掃く面積を可視化してみましょう。
import numpy as np
import matplotlib.pyplot as plt
from matplotlib.patches import Polygon
from matplotlib.collections import PatchCollection
# --- 楕円軌道のシミュレーション（再利用） ---
mu = 1.0
dt = 0.001
t_max = 8.0
t = np.arange(0, t_max, dt)
N = len(t)
r0 = np.array([1.0, 0.0])
v0 = np.array([0.0, 1.0]) # 円軌道速度 → ちょうど円になるので少し変更
v0 = np.array([0.0, 0.8]) # 楕円軌道（e ≈ 0.38）
r = np.zeros((N, 2))
v = np.zeros((N, 2))
r[0], v[0] = r0, v0
for i in range(N - 1):
r_mag = np.linalg.norm(r[i])
a = -mu / r_mag**3 * r[i]
v_half = v[i] + 0.5 * dt * a
r[i + 1] = r[i] + dt * v_half
r_mag_new = np.linalg.norm(r[i + 1])
a_new = -mu / r_mag_new**3 * r[i + 1]
v[i + 1] = v_half + 0.5 * dt * a_new
# --- 等時間間隔の扇形面積を可視化 ---
fig, ax = plt.subplots(1, 1, figsize=(8, 8))
# 軌道を描画
ax.plot(r[:, 0], r[:, 1], 'k-', linewidth=0.5, alpha=0.3)
# 等時間間隔（軌道周期の1/12ごと）のセグメントを色分け
# まず軌道周期を推定
Lz = r[0, 0] * v[0, 1] - r[0, 1] * v[0, 0]
a_semi = 1.0 / (2.0 / np.linalg.norm(r0) - np.linalg.norm(v0)**2 / mu)
T_orb = 2 * np.pi * np.sqrt(a_semi**3 / mu)
n_segments = 12
steps_per_segment = int(T_orb / dt / n_segments)
colors = plt.cm.tab20(np.linspace(0, 1, n_segments))
areas = []
for k in range(n_segments):
i_start = k * steps_per_segment
i_end = (k + 1) * steps_per_segment
if i_end >= N:
break
# 扇形の頂点: 原点 + 軌道上の点列
segment_r = r[i_start:i_end:max(1, steps_per_segment // 50)]
vertices = np.vstack([[0, 0], segment_r, [0, 0]])
polygon = Polygon(vertices, closed=True)
# 面積を計算（外積の和）
area = 0.0
for j in range(i_start, i_end - 1):
area += 0.5 * abs(r[j, 0] * r[j+1, 1] - r[j+1, 0] * r[j, 1])
areas.append(area)
ax.fill(vertices[:, 0], vertices[:, 1], color=colors[k], alpha=0.4)
ax.plot(0, 0, 'yo', markersize=15, zorder=5, label='Central body')
ax.set_xlabel('x', fontsize=13)
ax.set_ylabel('y', fontsize=13)
ax.set_title('Equal time intervals sweep equal areas (Kepler 2nd law)', fontsize=14)
ax.set_aspect('equal')
ax.legend(fontsize=11)
ax.grid(True, alpha=0.2)
plt.tight_layout()
plt.show()
# 面積の比較
print("各セグメントの面積:")
for k, area in enumerate(areas):
print(f" セグメント {k+1:2d}: {area:.6f}")
if areas:
print(f" 面積の標準偏差/平均: {np.std(areas)/np.mean(areas):.6f}")
このグラフでは、楕円軌道を等時間間隔で12個の扇形に分割しています。各扇形が異なる色で塗り分けられており、以下の特徴が読み取れます。
-
面積の一致: 近点付近（中心天体に近い領域）の扇形は「細長い」形をしていますが、遠点付近の扇形は「扁平で幅広い」形をしています。しかし、各扇形の面積は等しくなっており、面積の標準偏差/平均比は $10^{-4}$ 以下に収まります。
-
速度の変化: 近点付近では弧が長い（速度が大きい）ため細長くなり、遠点付近では弧が短い（速度が小さい）ため幅広になります。この速度変化がちょうど距離の変化を打ち消し、面積が一定に保たれるのです。
-
角運動量保存との対応: 面積速度 $dA/dt = L/(2m) = \text{const}$ が成り立っているため、等時間間隔で掃かれる面積 $\Delta A = (L/2m)\Delta t$ は一定になります。これはケプラーの第2法則そのものです。
角運動量保存の破れ: トルクの効果
角運動量が保存されないケース、すなわちトルクが存在する場合の効果を可視化してみましょう。中心力に摂動的な接線方向の力（たとえば大気抵抗や推進力）を加えると、角運動量が時間とともに変化する様子を確認できます。
import numpy as np
import matplotlib.pyplot as plt
# --- 中心力 + 接線方向の摂動力 ---
mu = 1.0
dt = 0.001
t_max = 30.0
t = np.arange(0, t_max, dt)
N = len(t)
r0 = np.array([1.0, 0.0])
v0 = np.array([0.0, 1.0])
# ケース1: 摂動なし（角運動量保存）
r1 = np.zeros((N, 2))
v1 = np.zeros((N, 2))
r1[0], v1[0] = r0, v0.copy()
Lz1 = np.zeros(N)
# ケース2: 弱い接線方向の制動力（大気抵抗の模擬）
r2 = np.zeros((N, 2))
v2 = np.zeros((N, 2))
r2[0], v2[0] = r0, v0.copy()
Lz2 = np.zeros(N)
drag_coeff = 0.005
for i in range(N - 1):
# ケース1: 中心力のみ
rm1 = np.linalg.norm(r1[i])
a1 = -mu / rm1**3 * r1[i]
vh1 = v1[i] + 0.5 * dt * a1
r1[i+1] = r1[i] + dt * vh1
rm1n = np.linalg.norm(r1[i+1])
a1n = -mu / rm1n**3 * r1[i+1]
v1[i+1] = vh1 + 0.5 * dt * a1n
# ケース2: 中心力 + 制動力
rm2 = np.linalg.norm(r2[i])
a_grav = -mu / rm2**3 * r2[i]
v_mag = np.linalg.norm(v2[i])
if v_mag > 0:
a_drag = -drag_coeff * v_mag * v2[i] / v_mag # 速度と逆方向
else:
a_drag = np.zeros(2)
a2 = a_grav + a_drag
vh2 = v2[i] + 0.5 * dt * a2
r2[i+1] = r2[i] + dt * vh2
rm2n = np.linalg.norm(r2[i+1])
a_grav_n = -mu / rm2n**3 * r2[i+1]
v2_mag_n = np.linalg.norm(vh2)
if v2_mag_n > 0:
a_drag_n = -drag_coeff * v2_mag_n * vh2 / v2_mag_n
else:
a_drag_n = np.zeros(2)
a2n = a_grav_n + a_drag_n
v2[i+1] = vh2 + 0.5 * dt * a2n
# 角運動量の計算
Lz1 = r1[:, 0] * v1[:, 1] - r1[:, 1] * v1[:, 0]
Lz2 = r2[:, 0] * v2[:, 1] - r2[:, 1] * v2[:, 0]
# --- 可視化 ---
fig, axes = plt.subplots(1, 3, figsize=(18, 5))
# (1) 軌道比較
axes[0].plot(r1[:, 0], r1[:, 1], 'b-', linewidth=0.5, alpha=0.5, label='No drag')
axes[0].plot(r2[:, 0], r2[:, 1], 'r-', linewidth=0.5, alpha=0.5, label='With drag')
axes[0].plot(0, 0, 'yo', markersize=12)
axes[0].set_xlabel('x', fontsize=12)
axes[0].set_ylabel('y', fontsize=12)
axes[0].set_title('Orbit comparison', fontsize=14)
axes[0].set_aspect('equal')
axes[0].legend(fontsize=11)
axes[0].grid(True, alpha=0.3)
# (2) 角運動量の比較
axes[1].plot(t, Lz1, 'b-', linewidth=1.5, label='No drag ($L_z$ = const)')
axes[1].plot(t, Lz2, 'r-', linewidth=1.5, label='With drag ($L_z$ decreases)')
axes[1].set_xlabel('Time', fontsize=12)
axes[1].set_ylabel('$L_z$', fontsize=12)
axes[1].set_title('Angular momentum comparison', fontsize=14)
axes[1].legend(fontsize=11)
axes[1].grid(True, alpha=0.3)
# (3) 軌道半径の変化
r1_mag = np.linalg.norm(r1, axis=1)
r2_mag = np.linalg.norm(r2, axis=1)
axes[2].plot(t, r1_mag, 'b-', linewidth=1, alpha=0.5, label='No drag')
axes[2].plot(t, r2_mag, 'r-', linewidth=1, alpha=0.5, label='With drag')
axes[2].set_xlabel('Time', fontsize=12)
axes[2].set_ylabel('$|\\mathbf{r}|$', fontsize=12)
axes[2].set_title('Orbital radius vs time', fontsize=14)
axes[2].legend(fontsize=11)
axes[2].grid(True, alpha=0.3)
plt.tight_layout()
plt.show()
制動力（大気抵抗）が角運動量にどのような影響を与えるかが明確に見てとれます。
-
軌道の螺旋化: 制動力がある場合（赤線）、軌道が徐々に内側に巻き込む螺旋を描きます。これは人工衛星の軌道減衰に対応する現象で、低軌道衛星が大気抵抗により最終的に大気圏に再突入するメカニズムです。
-
角運動量の単調減少: 制動力は速度と逆方向にはたらくため、常に負のトルクを生み出します。その結果、角運動量は単調に減少していきます。一方、中心力のみの場合は角運動量が完全に保存されています。
-
角運動量保存の「破れ」の物理的意味: 角運動量の変化は $d\bm{L}/dt = \bm{\tau}$ で決まりますが、中心力以外の成分（ここでは制動力）がトルクを生み出すために角運動量が変化します。ホーマン遷移軌道のような軌道変換でも、推進力によるトルクで角運動量を意図的に変化させます。
角運動量の数値計算ユーティリティ
最後に、3次元での角運動量計算やトルクの計算を行う汎用的なユーティリティ関数をまとめておきます。
import numpy as np
def angular_momentum(r, p):
"""
角運動量ベクトル L = r × p を計算する。
Parameters
----------
r : array_like, shape (3,) or (N, 3)
位置ベクトル [m]
p : array_like, shape (3,) or (N, 3)
運動量ベクトル [kg·m/s]
Returns
-------
L : ndarray, shape (3,) or (N, 3)
角運動量ベクトル [kg·m²/s]
"""
return np.cross(r, p)
def torque(r, F):
"""
トルク（力のモーメント）τ = r × F を計算する。
Parameters
----------
r : array_like, shape (3,) or (N, 3)
作用点の位置ベクトル [m]
F : array_like, shape (3,) or (N, 3)
力ベクトル [N]
Returns
-------
tau : ndarray, shape (3,) or (N, 3)
トルクベクトル [N·m]
"""
return np.cross(r, F)
def areal_velocity(r, v):
"""
面積速度 dA/dt = |r × v| / 2 を計算する。
Parameters
----------
r : array_like, shape (3,) or (N, 3)
位置ベクトル
v : array_like, shape (3,) or (N, 3)
速度ベクトル
Returns
-------
dA_dt : float or ndarray
面積速度
"""
cross = np.cross(r, v)
if cross.ndim == 1:
return 0.5 * np.linalg.norm(cross)
else:
return 0.5 * np.linalg.norm(cross, axis=1)
# --- 使用例 ---
# 3次元での角運動量計算
r_3d = np.array([1.0, 0.0, 0.0])
m = 2.0
v_3d = np.array([0.0, 3.0, 1.0])
p_3d = m * v_3d
L = angular_momentum(r_3d, p_3d)
print(f"位置: r = {r_3d}")
print(f"運動量: p = {p_3d}")
print(f"角運動量: L = {L}")
print(f"|L| = {np.linalg.norm(L):.4f}")
# 中心力のトルク（ゼロになることの確認）
F_central = -1.0 / np.linalg.norm(r_3d)**3 * r_3d # 万有引力
tau_central = torque(r_3d, F_central)
print(f"\n中心力: F = {F_central}")
print(f"トルク: τ = {tau_central}")
print(f"|τ| = {np.linalg.norm(tau_central):.2e} (≈ 0, as expected)")
このユーティリティコードは、角運動量とトルクの計算を関数としてまとめたものです。np.cross を使うことで、3次元ベクトルの外積を簡潔に計算できます。中心力に対するトルクがゼロ（数値誤差の範囲内）になることも確認でき、理論と一致しています。これらの関数は、ビザビバ方程式やホーマン遷移軌道のシミュレーションにもそのまま流用できます。
まとめ
本記事では、角運動量の定義からトルクとの関係、保存則、質点系への拡張、剛体への橋渡しまでを体系的に解説しました。要点を整理します。
- 角運動量の定義: $\bm{L} = \bm{r} \times \bm{p} = \bm{r} \times m\bm{v}$。回転中心からの距離と速度の接線成分の積に質量を掛けたもので、回転運動の勢いを表します
- トルクと角運動量の運動方程式: $d\bm{L}/dt = \bm{\tau}$。これはニュートンの第2法則 $d\bm{p}/dt = \bm{F}$ の回転版であり、トルクが角運動量を変化させます
- 角運動量保存則: $\bm{\tau} = \bm{0}$ のとき $\bm{L} = \text{const}$。中心力場（万有引力、クーロン力）ではトルクがゼロとなり、角運動量が保存されます
- ケプラーの第2法則: 面積速度 $dA/dt = |\bm{L}|/(2m)$ が一定になることと同値であり、角運動量保存の直接的な帰結です
- 質点系: 内力のトルクは作用反作用の法則によりキャンセルされ、全角運動量は外力のトルクのみで変化します
- 剛体: 固定軸まわりでは $L = I\omega$ となり、角運動量保存から $I$ の減少が $\omega$ の増大をもたらします（スケータースピン）
- 歳差運動: トルクが角運動量ベクトルの方向を変えるが大きさは変えない場合、自転軸が歳差運動を行います
角運動量は古典力学の枠組みを超えて、量子力学（スピン角運動量、軌道角運動量）、天体力学（軌道力学）、工学（姿勢制御、ジャイロスコープ）にまで広がる普遍的な物理量です。
次のステップとして、以下の記事も参考にしてください。
- 慣性モーメント — 角運動量と角速度を結ぶ $L = I\omega$ の「$I$」を深掘りします
- オイラーの剛体運動方程式 — 3次元空間での剛体の回転運動を記述する方程式です
- 万有引力とケプラーの法則 — 角運動量保存が導くケプラーの3法則を詳しく解説しています
- ビザビバ方程式 — 角運動量とエネルギー保存を組み合わせた軌道力学の基本方程式です
- ホーマン遷移軌道 — 角運動量の変化を利用した軌道変換の設計法です