力学編 第9章

自由な座標

拘束された運動の場合、「デカルト座標 𝒙 」ではなく「拘束に沿った自由な座標 𝜽 」に対する運動方程式を考えることもできる。具体的には、運動方程式は式(10)で与えられる。例えば、振り子の場合は式(12)である。

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

自由な座標 𝜃 での運動方程式(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.1-注1】自由な座標 𝜽 に対する運動方程式:式(10)

座標 𝑿 が、時間 𝑡 と自由な座標 𝜽 の関数 𝑿⁡(𝑡,𝜽) で表せる時、 𝜽 に関する運動方程式は以下のようになる: ¨𝜽=𝜇−1⁢𝜕𝑿𝜕𝜽T⁢[𝑭−𝑀⁡(𝑑𝑑⁢𝑡⁢𝜕𝑿𝜕𝜽)⁢˙𝜽⁢−𝑀⁡𝑑𝑑⁢𝑡⁢𝜕𝑿𝜕𝑡]𝜇≡𝜕𝑿𝜕𝜽T⁢𝑀⁡𝜕𝑿𝜕𝜽(10) ただし、 𝑭 は重力などの外力である。拘束条件が時間に依存しない場合、即ち、 𝑿⁡(𝜽) と書ける場合、緑字部分が消えて、以下のようになる: ¨𝜽=𝜇−1⁢𝑑⁢𝑿𝑑⁢𝜽T⁢[𝑭−𝑀⁡(𝑑𝑑⁢𝑡⁢𝑑⁢𝑿𝑑⁢𝜽)⁢˙𝜽]𝜇≡𝑑⁢𝑿𝑑⁢𝜽T⁢𝑀⁡𝑑⁢𝑿𝑑⁢𝜽(11)

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⁡𝜃+˙𝜙2⁢sin⁡𝜃⁢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]=⎡⎢ ⎢ ⎢ ⎢⎣−𝑙1⁢sin⁡𝜃1−𝑙1⁢cos⁡𝜃1−𝑙1⁢sin⁡𝜃1+𝑙2⁢sin⁡𝜃2−𝑙1⁢cos⁡𝜃1−𝑙2⁢cos⁡𝜃2⎤⎥ ⎥ ⎥ ⎥⎦𝑑⁢𝑿𝑑⁢𝜽=[𝜕𝑿𝜕𝜃1𝜕𝑿𝜕𝜃2]=⎡⎢ ⎢ ⎢ ⎢⎣𝑙1⁢cos⁡𝜃10𝑙1⁢sin⁡𝜃10𝑙1⁢cos⁡𝜃1𝑙2⁢cos⁡𝜃2𝑙1⁢sin⁡𝜃1𝑙2⁢sin⁡𝜃2⎤⎥ ⎥ ⎥ ⎥⎦𝑑𝑑⁢𝑡⁢𝑑⁢𝑿𝑑⁢𝜽=⎡⎢ ⎢ ⎢ ⎢ ⎢⎣−𝑙1⁢˙𝜃1⁢sin⁡𝜃10−𝑙1⁢˙𝜃1⁢cos⁡𝜃10−𝑙1⁢˙𝜃1⁢sin⁡𝜃1−𝑙2⁢˙𝜃2⁢sin⁡𝜃2−𝑙1⁢˙𝜃1⁢cos⁡𝜃1−𝑙2⁢˙𝜃2⁢cos⁡𝜃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) に代入すればよい。特にきれいになるわけではないので、結果は割愛する。

シミュレーション

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

【9.2-注1】3次のベジェ曲線

ベジェ曲線は、右図のように、2つの任意の点 𝒑1,𝒑2 の間を結ぶ曲線を表すために用いられる。具体的には、曲線の形状をコントロールするためのハンドルと呼ばれる2つの点 𝒉1,𝒉2 を用いて、以下のように定義される:( 𝒙⁡(0)=𝒑1,𝒙⁡(1)=𝒑2 ) 𝒙⁡(𝜃)=(1−𝜃)3⁢𝒑1+3⁢𝜃⁢(1−𝜃)2⁢𝒉1+3⁢𝜃2⁢(1−𝜃)⁢𝒉2+𝜃3⁢𝒑2(16)

補足

ベジェ曲線の連結

複雑な曲線を補間する際には、曲線を分割し、各区間をベジェ曲線で補間するという方法がとられる。右図のように、2つのベジェ曲線が点 𝒑2=𝒑′1 においてなめらかに接続する条件は、 𝒉2−𝒑2 と 𝒉′1−𝒑′1 が反平行になることである。