自由な座標

物体の運動を、拘束に沿った自由な座標 𝜽 を使って記述したい。

自由な座標 𝜃 での運動方程式(1)を求めたい

適当に曲げた針金を用意し、そこに物体を通す（右図）。針金に沿ったその物体の運動 𝒙𝑡 が知りたい。

拘束条件を書き下すのが難しい

これは、すでに第7章で扱った通常の拘束系であるが、1つ問題がある。これまでの章では、「拘束条件が 𝑮(𝒙)=𝟎 という形の方程式で表される」ことを前提としていたわけだが、今の場合、適当に曲げた針金の形状を表す 𝑮 を見つけることが困難なのである。

パラメータ表示のほうが自然

このような場合、右図のように、曲線上にパラメータ 𝜃 を取ることにより、物体の位置を 𝜃 で表すほうが簡単である。例えば2次元の振り子であれば、振れ角 𝜃 が使えるし、一般の形状であれば、ベジェ曲線（以下の9.2節の【9.2-注1】）が使える。この 𝜃 には拘束がかかっておらず自由な値をとれるので、 𝜃 を、自由な座標と呼ぶことにする。この用語は一般的ではない。一般化座標と呼ぶことが多いが、こちらは文字通りデカルト座標に限定しない一般の座標系を考える際の用語であり、拘束条件とは無関係の文脈でも用いられる。

自由な座標 𝜃 に対する運動方程式を求めたい

自由な座標 𝜃 の値を決めれば、物体の位置 𝒙 が求まる： 𝒙(𝜃) 。よって、物体の運動 𝒙𝑡 は、 𝜃 の時間変化 𝜃𝑡 から求めることができる： 𝒙𝑡=𝒙(𝜃𝑡) 。 𝜃𝑡 を直接計算する方法は難しくない。実際、前章までで扱ってきた運動方程式 𝑚¨𝒙=[⋯] を、 ¨𝜃 に対する運動方程式： ¨𝜃=[⋯](1) の形に変形するだけである（純粋に数学的な問題である）。あとは、この式(1)を任意の初期値 𝜃0,˙𝜃0 のもとで解けば、 𝜃𝑡 が求まり、 𝒙𝑡 も分かることになる。

この章の方針

そこでこの章では、「運動方程式(1)の導出」と、それを用いた「物体の運動の計算」の2つの節に分けて議論を行う：9.1自由な座標 𝜽 での運動方程式9.2物体の運動 𝑿𝑡 の計算この章からは、これまでのような「力を求めれば運動方程式を使って運動が計算できる」という素朴な考え方から少し離れ、応用的な話題となる。

9.1自由な座標 𝜽 での運動方程式

この節では、自由な座標 𝜽 を用いた場合の運動方程式(10)を導出する。これは、単なる座標変換である。

9.1.1求めたい式の確認

実際に計算に入る前に、やりたいことをまとめておく。

自由度・粒子が複数ある場合を考慮する

針金の場合、自由な座標 𝜃 の成分は1つだけだが、曲面拘束や多粒子の場合には、複数成分になる。そこで、以後、ベクトル表記 𝜽 を用いることにする。 𝜽 の成分の数を自由度という。例えば、次節で扱う平面上の二重振り子（右図）の場合、自由度は2である。物体の位置 𝒙 についても、これまで通り、複数の物体 𝒙1,𝒙2,… がある場合も考慮して、 𝑿 と書くことにする： 𝑿≡⎡⎢ ⎢⎣𝒙1𝒙2⋮⎤⎥ ⎥⎦

元となる運動方程式：式(3)

さて、前章で述べたように、拘束条件 𝑮(𝑡,𝑿)=𝟎(2) 時刻 𝑡 に依存しているのは、「拘束条件が時間とともに変化する場合」も考慮しているからである。が課されている時、運動方程式は、以下の形で書ける： 𝑀¨𝑿=𝑭+𝜕𝑮𝜕𝑿T𝝀⏟𝑭c(3) 右辺第2項は、拘束力 𝑭c が拘束面と「垂直」であることを示している。 𝑭c は求めることができ、具体的な形は前章の【8.2-注4】で与えた。

導出したい運動方程式：式(4)

この式(3)を変形して、 𝜽 に対する運動方程式 ¨𝜽=[⋯](4) の形にできれば、 𝜽𝑡 が計算できるようになる。これを求めたいわけである。

9.1.2自由な座標 𝜽 に対する運動方程式：式(10)

それでは実際に計算しよう。

基本的な方針

まず、拘束が時間とともに変化する場合(2)を考えているので、座標 𝑿 は、時刻 𝑡 と自由な座標 𝜽 の関数になる： 𝑿=𝑿(𝑡,𝜽)(5) 例えば、第7章の7.4節で扱った「動く中心点 𝒂(𝑡) を持つ振り子」の場合だと 𝒙(𝑡,𝜃)=[−𝑙sin⁡𝜃−𝑙cos⁡𝜃]−𝒂(𝑡) となる。この式(5)の両辺を 𝑡 で2階微分すれば、左辺は ¨𝑿 となり、右辺は ¨𝜽 含む。よって、それを運動方程式(3)に代入して ¨𝑿 を消去し、 ¨𝜽 について解けば、目的の式(4)の形になる。

¨𝜽 に対する運動方程式

まず、式(5)の両辺を時刻 𝑡 で2階微分すると、以下のようになる： ˙𝑿=𝑑𝑑𝑡𝑿(𝑡,𝜽)=𝜕𝑿𝜕𝑡+𝜕𝑿𝜕𝜽˙𝜽¨𝑿=𝑑𝑑𝑡(𝜕𝑿𝜕𝑡+𝜕𝑿𝜕𝜽˙𝜽)=𝑑𝑑𝑡𝜕𝑿𝜕𝑡+(𝑑𝑑𝑡𝜕𝑿𝜕𝜽)˙𝜽+𝜕𝑿𝜕𝜽¨𝜽(6) これが、 ¨𝑿,¨𝜽 の関係式である。この式を運動方程式(3)に代入すると 𝑀[𝑑𝑑𝑡𝜕𝑿𝜕𝑡+(𝑑𝑑𝑡𝜕𝑿𝜕𝜽)˙𝜽+𝜕𝑿𝜕𝜽¨𝜽]=𝑭+𝑭c となる。これを ¨𝜽 について解くには、両辺に 𝜕𝑿𝜕𝜽T を左乗し： 𝜕𝑿𝜕𝜽T𝑀𝜕𝑿𝜕𝜽⏟__⏟__⏟≡𝜇¨𝜽=𝜕𝑿𝜕𝜽T(𝑭+𝑭c−𝑀[𝑑𝑑𝑡𝜕𝑿𝜕𝑡+(𝑑𝑑𝑡𝜕𝑿𝜕𝜽)˙𝜽]) 𝜇−1 を左乗すればよい： ¨𝜽=𝜇−1𝜕𝑿𝜕𝜽T(𝑭+𝑭c−𝑀[𝑑𝑑𝑡𝜕𝑿𝜕𝑡+(𝑑𝑑𝑡𝜕𝑿𝜕𝜽)˙𝜽])(7) 正方行列 𝜇 が逆行列を持つことを示しておこう。「ある正方行列 𝐴 が逆行列を持つための条件」には様々な表現が知られているが、ここでは「 𝒂≠𝟎 ならば 𝐴𝒂≠𝟎 」を採用する。まず、 𝜕𝑿𝜕𝜽 は列フルランク（＝列が1次独立）なので、 𝒂≠𝟎 であれば 𝜕𝑿𝜕𝜽𝒂≡𝒃≠𝟎 である。よって、 𝒂T𝜇𝒂=𝒃T𝑀𝒃>0 、即ち、 𝜇𝒂≠𝟎 。 ◼

拘束条件を使うと、拘束力 𝑭c が消せるのでは？

ところで、拘束力は「拘束条件を満たすために働く力」である。自由な座標 𝜽 は、自動的に拘束条件を満たすのだから、 𝜽 に対する運動方程式には拘束力は現れないのではないだろうか。即ち式(7)の拘束力 𝑭c を含む項は、 𝟎 になることが期待される。どうすればこれを示せるだろうか。ここまでは、式(5)だけを使っており、「 𝜽 が自由な座標である」という条件(2)： 𝑮(𝑡,𝑿(𝑡,𝜽))=𝟎(8) は暗黙的に（「 𝜽 を使って任意の 𝑿 を表現できる」という形で）しか使っていない。そこで、この条件式(8)に着目してみよう。

拘束条件の微分形：式(9)

拘束条件(8)の微分を計算してみる。微小量 𝛿𝑡,𝛿𝜽 に対し、 𝑮 の変化 𝛿𝑮 は、1次近似の範囲で 𝟎 となるのだから： 𝟎=𝛿𝑮(𝑡,𝑿(𝑡,𝜽))≐𝜕𝑮𝜕𝑡𝛿𝑡+𝜕𝑮𝜕𝑿𝛿𝑿(𝑡,𝜽)≐𝜕𝑮𝜕𝑡𝛿𝑡+𝜕𝑮𝜕𝑿(𝜕𝑿𝜕𝑡𝛿𝑡+𝜕𝑿𝜕𝜽𝛿𝜽)∣𝛿𝑡,𝛿𝜽を括り出す=(𝜕𝑮𝜕𝑡+𝜕𝑮𝜕𝑿𝜕𝑿𝜕𝑡)𝛿𝑡+𝜕𝑮𝜕𝑿𝜕𝑿𝜕𝜽𝛿𝜽 これが任意の 𝛿𝑡,𝛿𝜽 で成り立つので、赤字部分と緑字部分はそれぞれゼロとなる： 𝜕𝑮𝜕𝑡+𝜕𝑮𝜕𝑿𝜕𝑿𝜕𝑡=0𝜕𝑮𝜕𝑿𝜕𝑿𝜕𝜽=𝟎⎫{ { {⎬{ { {⎭(9) これが、微分形での拘束条件である。なお、導出方法からも分かるように、第1式の 𝜕𝑮𝜕𝑡 は、 𝑮 を 𝑡,𝑿 の関数と見た時の 𝑡 での偏微分である（ 𝑡,𝜽 の関数と見た時ではない）。

運動方程式(7)から拘束力 𝑭c が消せる

拘束条件(9)において重要なのは、 𝜽 を含んでいる第2式である。これを使うと、運動方程式(7)から拘束力 𝑭c が消せる。実際、 𝑭c を含む項は 𝜕𝑿𝜕𝜽T𝑭c=𝜕𝑿𝜕𝜽T𝜕𝑮𝜕𝑿T𝝀∵ダランベールの原理=(𝜕𝑮𝜕𝑿𝜕𝑿𝜕𝜽)T𝝀∵𝐴T𝐵T=(𝐵𝐴)T=𝟎 となる。以上をまとめると、自由な座標 𝜽 は以下の【9.1-注1】の式(10)のようになる。

9.2物体の運動 𝑿𝑡 の計算

以上で、必要な議論がそろった。この節では、まず、自由な座標 𝜽 を用いた場合の運動の計算方法についてまとめる。その後、振り子と二重振り子の場合について、自由な座標 𝜽 を使った形で書き下す。これらについては、それぞれ第7章の7.3.2節と第8章の8.3.3節で、デカルト座標での運動方程式を示している。最後に、冒頭で触れた任意形状の針金上の運動の計算を行う。

9.2.1物体の運動 𝑿𝑡 の計算方法

物体の位置 𝑿 は、自由な座標 𝜽 から決まるのだから、物体の運動 𝑿𝑡 を知るには、 𝜽 の時間変化 𝜽𝑡 を求めれば十分である。実際に 𝜽𝑡 を求めるには、1 𝜽 を設定した上で、2初期値 𝜽0,˙𝜽0 を任意に与え、3運動方程式(10)を解けばよい。式(10)は、通常の運動方程式と同じ2階の微分方程式なので、これまでと同様にオイラー法（第3章の【3.1-注1】参照）などを用いて、数値的に解くことができる。

9.2.2例題1極座標での振り子の運動方程式：式(12)

2次元平面上の2次元振り子を考える。

自由な座標 𝜃 の設定

自由な座標 𝜃 として、右図のような極座標をとる（おもりが垂れ下がった位置が 𝜃=0 ）。座標 𝒙(𝜃) およびその微分は 𝒙(𝜃)=[𝑥𝑦]=[−𝑙sin⁡𝜃−𝑙cos⁡𝜃]𝑑𝒙𝑑𝜃=𝑙[cos⁡𝜃sin⁡𝜃]𝑑𝑑𝑡𝑑𝒙𝑑𝜃=𝑙˙𝜃[−sin⁡𝜃−cos⁡𝜃] となる。

運動方程式：式(12)

これらを運動方程式(11)：（ 𝑔=|𝒈| ） ¨𝜃=(𝑚∣𝑑𝒙𝑑𝜃∣2)−1𝑑𝒙𝑑𝜃T([0−𝑚𝑔]−𝑚𝑑𝑑𝑡𝑑𝒙𝑑𝜃˙𝜃) に代入すれば、極座標での運動方程式が得られる（緑字部分の項は消える）： ¨𝜃=−𝑔𝑙sin⁡𝜃(12) 数値計算については、既に第7章の7.3.2節で示した「デカルト座標で表したもの」と見かけ上変わらないので、割愛する。

参考3次元の場合

なお、3次元の場合、自由な座標は、方位角を 𝜙 として 𝒙(𝜃,𝜙)=⎡⎢ ⎢ ⎢⎣−𝑙sin⁡𝜃cos⁡𝜙−𝑙sin⁡𝜃sin⁡𝜙−𝑙cos⁡𝜃⎤⎥ ⎥ ⎥⎦ とすればよい。2次元の場合と同様の計算により、運動方程式は以下のようになる： [¨𝜃¨𝜙]=⎡⎢ ⎢ ⎢⎣−𝑔𝑙sin⁡𝜃+˙𝜙2sin⁡𝜃cos⁡𝜃−2˙𝜃˙𝜙cos⁡𝜃sin⁡𝜃⎤⎥ ⎥ ⎥⎦ 2次元の場合の運動方程式(12)に比べて複雑な式になり、 𝜃=0 のところに特異性がある（分母が 0 になる）。一方、デカルト座標で表せば、第7章の7.3.2節で見たように、2次元・3次元どちらでも同じ形になり、特異性もない。このように、自由な座標を使ったからといって必ずしも式が簡単になるわけではない。

9.2.3例題2極座標での二重振り子の運動方程式：式(14)

2次元平面上の二重振り子を考える。

自由な座標 𝜃 の設定

4次元2拘束なので、自由度は 4−2=2 である。よって、自由な座標 𝜽 は2成分である： 𝜽=[𝜃1𝜃2]T 。 𝜃1,𝜃2 を右図のようにとると、座標 𝑿(𝜽) およびその微分は 𝑿(𝜽)=[𝒙1𝒙2]=⎡⎢ ⎢ ⎢ ⎢⎣−𝑙1sin⁡𝜃1−𝑙1cos⁡𝜃1−𝑙1sin⁡𝜃1+𝑙2sin⁡𝜃2−𝑙1cos⁡𝜃1−𝑙2cos⁡𝜃2⎤⎥ ⎥ ⎥ ⎥⎦𝑑𝑿𝑑𝜽=[𝜕𝑿𝜕𝜃1𝜕𝑿𝜕𝜃2]=⎡⎢ ⎢ ⎢ ⎢⎣𝑙1cos⁡𝜃10𝑙1sin⁡𝜃10𝑙1cos⁡𝜃1𝑙2cos⁡𝜃2𝑙1sin⁡𝜃1𝑙2sin⁡𝜃2⎤⎥ ⎥ ⎥ ⎥⎦𝑑𝑑𝑡𝑑𝑿𝑑𝜽=⎡⎢ ⎢ ⎢ ⎢ ⎢⎣−𝑙1˙𝜃1sin⁡𝜃10−𝑙1˙𝜃1cos⁡𝜃10−𝑙1˙𝜃1sin⁡𝜃1−𝑙2˙𝜃2sin⁡𝜃2−𝑙1˙𝜃1cos⁡𝜃1−𝑙2˙𝜃2cos⁡𝜃2⎤⎥ ⎥ ⎥ ⎥ ⎥⎦ となる。

運動方程式：式(14)

これらを運動方程式(11)：（ 𝑔=|𝒈| ） ¨𝜽=(𝑑𝑿𝑑𝜽T𝑀𝑑𝑿𝑑𝜽)−1𝑑𝑿𝑑𝜽T⎛⎜ ⎜ ⎜ ⎜ ⎜ ⎜⎝⎡⎢ ⎢ ⎢ ⎢⎣0−𝑚1𝑔0−𝑚2𝑔⎤⎥ ⎥ ⎥ ⎥⎦−𝑀(𝑑𝑑𝑡𝑑𝑿𝑑𝜽)˙𝜽⎞⎟ ⎟ ⎟ ⎟ ⎟ ⎟⎠(13) に代入すれば、運動方程式が得られる： ¨𝜽=[1−cos⁡𝜃12−𝑙−1𝑟cos⁡𝜃12𝑙−1𝑟𝑚𝑟]𝑚𝑟−cos2⁡𝜃12[𝑚𝑟sin⁡𝜃1−𝑙𝑟˙𝜃22sin⁡𝜃2˙𝜃21][−𝑙−11𝑔sin⁡𝜃12]𝑙𝑟=𝑙1𝑙2,𝑚𝑟=𝑚1+𝑚2𝑚2,𝜃12=𝜃1−𝜃2(14) 数値計算をする場合には、行列計算はプログラム上で行えばよいので、式(14)の行列の積を手計算する必要はない（そもそも式(13)だけで十分）。数値計算結果については、既に第8章の8.3.3節で示した「デカルト座標で表したもの」と見かけ上変わらないので、割愛する。

9.2.4例題3ベジェ曲線上のおもりの運動方程式

任意の形状の曲線を拘束条件とする場合（例えば右図）、ベジェ曲線を用いることができる（以下の【9.2-注1】）。ここでは、外力として重力が働いているとする。

ベジェ曲線上の運動方程式：式(15)

運動方程式を求めるには、 𝒙(𝜃) およびその微分： 𝒙(𝜃)=式(16)𝑑𝒙𝑑𝜃=−3(1−𝜃)2𝒑1+3(1−3𝜃)(1−𝜃)𝒉1+3𝜃(2−3𝜃)𝒉2+3𝜃2𝒑2𝑑𝑑𝑡𝑑𝒙𝑑𝜃=6˙𝜃[(1−𝜃)𝒑1+(−2+3𝜃)𝒉1+(1−3𝜃)𝒉2+𝜃𝒑2] を、運動方程式(11)：（ 𝑔=|𝒈| ） ¨𝜃=(𝑚∣𝑑𝒙𝑑𝜃∣2)−1𝑑𝒙𝑑𝜃T([0−𝑚𝑔]−𝑚𝑑𝑑𝑡𝑑𝒙𝑑𝜃˙𝜃)(15) に代入すればよい。特にきれいになるわけではないので、結果は割愛する。

シミュレーション

数値計算を行うと、右図のようになる。