一般化座標系に於けるLagrangeの運動方程式の理解

arg maxarg min\providecommandrecterf\providecommand\providecommand\providecommandPr

経緯

 物理以外の参考書等(例えば制御工学)で拘束条件付きの問題に対して最初から一般化座標系を導入して拘束条件なしのLagrangeの運動方程式を立てている様子を見て、それが直角座標系に於けるNewtonの運動方程式と等価であることを納得したくて暫く考えていた。その結果をここに書き残しておく。

はじめに

 まず直角座標系に於いて拘束条件付きのLagrangeの運動方程式がNewtonの運動方程式と同値であることを理解する。そして、可逆な座標変換、所謂「一般化座標系」への変換を行ってもLagrangeの運動方程式が表面的には形を変えない(Lagrange関数Lの数式は変わるが、それを相変わらずLと書く限り形を変えない)ことを示す。そして拘束条件を表す制約式に現れる一般化座標変数の個数が最小となるような変換を選ぶと、一般化座標変数の一部が拘束条件から容易に求まることと、残りの一般化座標変数について拘束条件が消える可能性について述べる。最後に、そのことを確認できる例題を一つ紹介し、数値計算結果を示す。

問題設定

 Nを質点の個数とし、mi(i=1,,N)i番目の質点(以下「質点i」と呼ぶ)の質量とする。各質点に働く力は保存力、および質点の速度と直角に働く(つまり質点に対して仕事をしない)拘束力のみであるとする。ここでの拘束力とは、質点がある曲面や曲線に沿ってしか運動できず、かつそれら曲面,曲線を構成する物体(拘束物)と質点の間に摩擦力が生じない状況で、質点が拘束物から受ける力を想定する(例えば滑らかなワイヤに通されたビーズの運動や、半径が異なる、中心を共有する2つの滑らかな球面の間に挟まった鉄球の運動などを想像するとよい)。

 直角座標系では質点1個あたり3つの座標変数があるので、系全体では3N個の座標変数がある。これらをr=[r1,r2,,r3N]R3Nとする。質点iの座標は[r3i2,r3i1,r3i]となる。個々の質点の座標を分けて考えるのではなくrにまとめて表現することで、質点系を3N次元空間を運動する1個の高次元の質点として扱える。このようにしても数学的には元の運動方程式と変わらない。

保存力

 直角座標系に於いて、質点系に保存力をもたらすポテンシャルをVC:rR3NVC(r)Rとする(CはCartesianのC)(空間に固く固定された電荷などポテンシャルの原因が複数あるときはそれらの和をVCとすればよい)。質点系が受ける保存力はrVC(r)となる。

Lagrange関数

 直角座標系に於けるLagrange関数をLC:(t,r,r˙)R×R3N×R3NLC(t,r,r˙)Rとすると

LC(t,r,r˙)=i=1N12middt[r3i2,r3i1,r3i]22VC(r)=12r˙MCr˙VC(r)

ここにMCは「質量行列」と呼ばれる対角行列であり、M3i2,3i2=M3i1,3i1=M3i,3i=miである(直角座標系では対角であるが、他の座標系ではそうとも限らない)。

拘束条件

 質点系の拘束条件は複数あり得る。NCをその個数とする(CはConstraintのC)。これらの制約が、直角座標系においてC1級の関数gC,i:rR3NgC,i(r)R(i=1,,NC)を用いてgC,i(r)=0と表される場合を考える(gの添え字のCはCartesianのC)。但し各制約式gC,i(r)=0は、それが満足される領域でgC,i0となるように選ぶものとする(後で判るが、勾配が消えると計算上都合が悪い)。例えば質点1に対する制約条件として直線への拘束r12+r22=0を考えると、制約条件を満たす領域r1=r2=0で左辺の勾配[2r1,2r2,0]0になる。しかし同じ拘束条件を、平面を表す2つの制約式r1=0,r2=0の連立として表現すれば勾配はそれぞれ[1,0,0],[0,1,0]0となる。同様に曲線も曲面の交線として表現する。今回は不等式制約(物理的には壁に進入しない等に相当する)は扱わない。冒頭で仮定したように、質点が曲面から受ける拘束力は曲面に垂直、すなわちgC,iに平行である。

直角座標系に於けるLagrangeの運動方程式

 以上に述べたことから、直角座標系に於ける拘束条件付きのNewtonの運動方程式は次のようになり、これは拘束条件付きのLagrangeの運動方程式(Lagrange’s equations of the first kind (第一種Lagrange方程式)とも呼ばれる)と同値である。

(1)rLC(t,r,r˙)ddtr˙LC(t,r,r˙)+i=1NCλi(t)rgC,i(r)=0(2)gC,i(r)=0(i=1,,NC)

ここにλiは質点系の位置と加速度, 質点系に働く保存力に依存する、究極的にはtの関数である。質点が拘束条件のレールからはみ出ないようにgC,iに平行な力が働いて、質点系に働く慣性抗力(加速度と質量の積の-1倍)と質点系に働く保存力をキャンセルしていることは確かだが、いつどんな大きさで働くのかは自明ではない。そのためλirと同時に解かねばならない。

可逆な座標変換をしてもLagrangeの運動方程式の形が変わらないこと

 産業応用上解く必要に迫られる力学の問題は解析的に解けないものが殆どであり、数値解析に頼ることになるのだが、その場合でも前出の式(1),(2)をできるだけ解きやすくなるように座標変換することは必要である。拘束条件を利用して解くべき座標変数の個数を減らせれば問題の次元が下がって扱いやすくなる。

 rR3NqR3Nの間の可逆なC2級の座標変換を考える。つまりrqC2級関数r:qR3Nr(q)R3Nとして表せて、かつqC2級関数q:rR3Nq(r)R3Nとして表せる状況を考える。

 まずLCt,q,q˙の関数になる。なぜなら
ri˙=dri(q)dt=j=1Nriqjqj˙
となり、riqjqの関数なのでr˙q,q˙の関数となるからである。全く同様にしてq˙r,r˙の関数となる。後々の計算の為に座標系q上でのLagrange関数Lを次式で定義しておく。
L(t,q,q˙):=LC(t,r(q),r˙(q,q˙))
LLCは関数値は等しいが、写像としての定義が異なる。物理学ではそういうものを区別せずに、異なる座標系で同じ記号Lを使いまわすのが慣例だが、ここでは式変形から手品を排除して全てをはっきりと見る為に敢えて執拗に数学の形式に拘ることにする。(「いやしかし、rは既にtの関数としてr(t)と表しているのに、qの関数としてr(q)として表すのはいい加減ではないか?」という指摘もあろう。その通りで、有言実行するには新たに記号を定義して区別するしかないが、それでは却って式の可読性が下がる。結局、厳密さによる混乱の回避の恩恵と式の可読性のトレードオフで戦略を決めている。)LLCで表したのと同じ理屈で次式も成り立つ。LC(t,r,r˙)=L(t,q(r),q˙(r,r˙))

 座標変換によってrLCがどうなるか考える。

riLC(t,r,r˙)=riL(t,q(r),q˙(r,r˙))=j=1N(qjL(t,q,q˙))qjri+j=1N(qj˙L(t,q,q˙))qj˙ri

ここでN次正方行列AAi,j=qirjで定義すると、これは関数rのJacobianであり、関数rは可逆だからAは正則である。同様にN次正方行列B(正則とは限らない)をBi,j=qi˙rjで定義すると次式を得る。

(2)rLC=AqL+Bq˙L

 次に座標変換によってr˙LCがどうなるか考える。
ri˙LC(t,r,r˙)=ri˙L(t,q(r),q˙(r,r˙))=j=1N(qj˙L(t,q,q˙))qj˙ri˙
ここで
qi˙=dqi(r)dt=j=1Nqirjrj˙
より
qi˙rj˙=qirj=Ai,j
であるから次式を得る。
r˙LC=Aq˙L
これをさらにtで微分すると

(3)ddtr˙LC=(ddtA)q˙L+Addtq˙L

 次に座標変換によってrgC,iがどうなるか考える。Lを定義したときと同じように、gi(q):=gC,i(r(q))とするとgC,i(r)=gi(q(r))である。
rigC,i(r)=rig(q(r))=j=1Nqjg(q)qjri
よって次式が成り立つ。

(4)rgC,i=Aqgi

以上の式(2),(3),(4)(1)に適用すると次式を得る。
A(qLddtq˙L+i=1NCλi(t)qgi)+(B(ddtA))q˙L=0
BddtA=Oとなることを示す。まず
Bi,j=qi˙rj=rjdqi(r)dt=rjk=1Nqirkrk˙=k=1N2qirjrkrk˙
次に
dAi,jdt=ddtqirj=k=1N2qirkrjrk˙=k=1N2qirjrkrk˙=Bi,j
である。左から2つ目の等号は、qirjrの関数である(qirの関数なのでそれをrjで偏微分してもrの関数である)ことから従う。右から2つ目の等号はqrC2級関数であと仮定しておいたことから従う。以上2つの式からBddtA=Oであるので結局次式を得る。
A(qLddtq˙L+qgi)=0
先述の通りAは正則なので両辺に左からA1を掛けて次式を得る
qLddtq˙L+i=1NCλi(t)qgi=0
これと制約条件を合わせて、座標変換後の拘束条件付きのLagrangeの運動方程式は次式である。

(3)qLddtq˙L+i=1NCλi(t)qgi=0(4)gi(q)=0(i=1,,NC)

Lagrange関数と拘束条件式の記号が変わった以外は直角座標系のときと同じ形の式である。

座標変換を工夫すると制約なしのLagrangeの運動方程式が得られる可能性について

 式(3)3N行の連立微分方程式であるが、qの取り方を工夫することでqgkの幾つかの行を0にできる。q1,,q3Nのうちgk(q)(k=1,,NC)に現れるものの個数が最小になるようなqを見つけると都合が良い。もしqigk(q)に現れなければqgki行目は0である。そしてgk(q)=0(k=1,,NC)に現れるqの成分が簡単に求まる場合、残った成分について拘束なしのLagrangeの運動方程式(Euler–Lagrange方程式とかLagrange’s equations of the second kind (第二種Lagrange方程式)と呼ばれる)を解けば良いので問題が格段に扱いやすくなる。

 物理学を道具として用いる分野の技術書等(例えば制御工学)では普通、最初から拘束条件が消えるような一般化座標系への変換を考えてLagrange関数を構築して拘束なしのLagrangeの運動方程式を鮮やかに立てているが、これは「拘束条件が消えるような一般化座標系への変換を考え」る過程で自動的に「gk(q)=0(k=1,,NC)に現れるqの成分が簡単に求ま」って、そして「qigk(q)に現れなければqgki行目は0」となって、最初から拘束なしのLagrangeの運動方程式に辿り着ているのである。

問題例

 N=1,NC=1とし、質点の質量をmとする。拘束条件として潰れた楕円形の閉曲面0=gC(r):=r12+14r22+19r321を考える。r3軸の負の向きに重力加速度gが働くとする。

拘束条件の可視化

Lagrange関数は次式である。

LC(t,r,r˙)=12mr˙2gr3

直角座標系に於ける拘束条件付きのLagrangeの運動方程式を具体的に書き下すと次のようになる。

{[00g]m[r1¨r2¨r3¨]+λ(t)[2r112r229r3]=[000]r12+14r22+19r32=1

人間に解けるような問題ではないから計算機で数値計算するとして、拘束条件が厄介である。そこで次のような座標変換を考える。

{r1=q3sinq1cosq2r2=2q3sinq1sinq2r3=3q3cosq1q30

そうすると拘束条件は

0=g(q)=gC(r(q))=q321

となりq3=1が決まる。そしてqgは次式になる。

qg(q)=[002q3]

そしてLagrange関数をt,q,q˙の関数で表してq,q˙で偏微分してLagrangeの運動方程式を立てるのだが、実はここから先の計算は恐ろしく大変である。筆者の手作業では1日掛かっても終わりそうにないのでMathematicaにやらせたが、Lagrangeの運動方程式に至ってはここに書けないほど長い数式になった。簡単に表せる部分のみ具体的に表示し、他は記号的に書くと次のようになる。

(5)[q1L(t,q,q˙)q2L(t,q,q˙)q3L(t,q,q˙)]ddt[q1˙L(t,q,q˙)q2˙L(t,q,q˙)q3˙L(t,q,q˙)]+λ(t)[002q3]=[000](6)q3=1

(5)の3行目はL(t,q,q˙)q3=1を代入する前に計算しなければならないことに注意する(1,2行目はq3=1を先に代入してから計算しても構わない)。そもそも式(5)q3が変数であるとして導出されたのだから、(5)の計算前にq3を特定の値で固定してはいけない。例えばf:xx2という写像があって、何らかの条件からx(t)1であることが解っているからといってx(t)f(x(t))=0とするのはおかしい。x(t)f(x(t))=2x(t)なのでx(t)1であってもx(t)f(x(t))=20である。「関数の定義」と「関数値」は別物であって、「関数の定義」は「関数値」よりも多くの情報を保持している。偏微分は「関数の値」ではなく「関数の定義」に対して定められている。

(5)の上2行はqgの要素が0になっており、λ(t)の影響を受けない。既にq3(t)=1は解っていて未知の座標変数はq1(t),q2(t)の2つだけであるから式(5)の上2行からこれらが唯一に決定されるように思えるし、そうであって欲しい。そして実際にそうでなくてはならない。仮に式(5)の上2行だけでq1(t),q2(t)が唯一に決定されないとする。これ以外に利用できる情報は3行目の式しか残されていなが、それを考慮に入れたところで、λ(t)が同時に解くべき未知変数に追加されるため、結局解けないのである。古典力学的に正しく定義された問題に対しては解が唯一に決まることは疑いようのない事実であり、「q1(t),q2(t)が唯一に決まらない」状況は物理的に許されず、式(5)の上2行だけでq1(t),q2(t)が唯一に決定されるより他にあり得ないのである。そして3行目からλ(t)が決定される。

先にも述べた、解析力学の理解を前提とした技術書では座標変換で一般化座標変数の一部を制約条件から事前に求めてしまい(制約条件と同じ個数の座標変数が代数的に決まることを想定する)、それらを最初から定数とした上で、Lagrange関数を残りの一般化座標の関数と見做して拘束なしの運動方程式を解いている。こうしてλ(t)を表に出さずに解軌道に辿り着く。

この一般化座標系に於けるLagrangeの運動方程式をMathematicaで解いた結果を元の座標系に変換して可視化したものが次の図である。

数値計算で求めた軌道

初期条件はq1(0)=π/3,q2(0)=π/4,q1˙(0)=q2˙(0)=0すなわちr(0)=[6/4,6/2,3/2],r˙(0)=0である。この位置に固定されていた質点が時刻0に解放されて楕円曲面に拘束されながら重力によって落下し、底に辿り着くと勢いで壁を登っていき、速度が落ちて再び落下していくのを繰り返す様子が見える。

以下は、このシミュレーションを行ったMathematicaノートブックである。

投稿者: motchy

DSP and FPGA engineer working on measuring instrument.

コメントを残す