非ホロノミック拘束下のラグランジュ方程式を導く二つの流儀

新しい解析力学の教科書が出版された。
　近藤慶一著 『解析力学講義ー古典力学を超えてー』
非ホロノミック拘束の章があるので購入した。
私は非ホロノミック拘束系の応用に興味があって勉強中なのだがこれについて書かれているテキストは多くなく私にとってはMUST_BUYである。早速、該当の章を読んでみると、手持の邦書と違う運動方程式を導く手法が書かれているが、二つの手法は等価ではないことが知られている。
　実際自分の手で計算して異同を確かめたことはないのでやってみた。私の理解のレベルで簡略に書いているので詳細は以下の参考文献を参照していただきたい。
 この一文を書いている終盤でネット上にある演習問題の解答を見ると手持ちの邦書と同じ解き方で解説されている。

参考書籍
　近藤[2022] 近藤慶一著 『解析力学講義ー古典力学を超えてー』
　前野[2013] 前野昌弘著『よくわかる解析力学』
　山本[1998] 山本義隆・中村孔一著『解析力学Ⅰ』
　A.M.Bloch[2003] ”Nonholonomic Mechanics and Control"
参照論文
Flannery M R [2005] "The enigma of nonholonomic constraints"
(ズバリの内容である。珍しく微分幾何を使って書かれていないので助かる。)

その他関連論文
Leon A [2012] ”A historical review on nonholomic mechanics”
　この問題の議論の歴史的展開についてくわしく書いてある。
Fernandez and Bloch[2008] Equivalence of the dynamics of nonholonomic and variational nonholonomic systems for certain initial data
　表題からこれもズバリだが微分幾何全開で書かれていて今の私には読みこなせない。

１．本稿の主旨

非ホロノミック拘束を含む力学系について解析力学的に運動方程式を導くには二つの流儀がある。拘束条件由来の項をラグランジアンの段階で入れるか、ラグランジュ方程式にしてから入れるか。この二者について英語論文では長い議論が続いている。私が読み取ったところでは後者は物理的に正しい運動方程式を与えるが、前者はそうとは限らない。近藤[2022]は本文では前者を説き演習問題の解答では後者を説いている。前野[2013]山本[1998]は後者を説く。A.M.Bloch[2003]は両者を説く。　演習として、簡単な例で両者を使って運動方程式を導き比較してみる。

２．ホロノミック拘束下のラグランジアン

これについてはどの文献も同じことが書いてある。
ホロノミック拘束条件が要請される場合にラグランジアンから運動方程式（ラグランジュ方程式）を導く手法。（ホロノミック拘束条件：荒っぽく言えば速度項を含まず変位項だけからなる拘束条件）
一般座標を$${q_i}$$、その微分を$${\dot{q_i}}$$ とすればオイラー・ラグランジュ方程式はラグランジアン$${L}$$に対して

$$
{\frac{d}{d t}\left(\frac{\partial L}{\partial \dot q^i}\right)-\frac{\partial L}{\partial q^i}=0\hspace{20pt}-(1)}
$$

となる。拘束関係式を$${b^j(q)=0}$$ とし、対応するラグランジュ乗数を$${\lambda_j}$$と置けば、オイラー・ラグランジュ方程式は

$$
{\frac{d}{d t}\left(\frac{\partial L}{\partial \dot q^i}\right)-\frac{\partial L}{\partial q^i}=\sum_{j}{}\lambda_j\frac{\partial b^j}{\partial q^i}\hspace{20pt}-(2)}
$$

この式をよく見ると 拡張ラグランジアンを$${\tilde L}$$として

$$
　{\tilde L=L-\sum_{j}{}\lambda_j\ b^j\hspace{20pt}-(3)}
$$

と置くと

$$
{\frac{d}{d t}\left(\frac{\partial \tilde L}{\partial \dot q^i}\right)-\frac{\partial \tilde L}{\partial q^i}=0\hspace{20pt}-(4)}
$$

となって、拘束がないオイラー・ラグランジュ方程式と同じ形に書くことができる。

３．非ホロノミック拘束の取り扱い

非ホロノミック拘束は、ホロノミック拘束と同じ手法で取り扱えるのであろうか？ 　ここでは非ホロノミック拘束の中でも速度に線形な　

$$
{\sum_{k=1}^{n}{a^j_k(q)}{\dot q^k}=0\hspace{20pt}-(5)}
$$

と表される$${j=1 ... m}$$の$${m}$$個の拘束を考える。$${n}$$は一般座標qiの個数である。以下変位のみの関数$${a^j_k(q)}$$の引数を省略して$${a^j_k}$$と書く。
ここで各テキストの内容は分かれる。
3.1 ダランベールの原理に基づく方法
拘束条件由来の項をラグランジュ方程式に付加する方法。
前野[2013],山本[1998],近藤[2022]の演習問題の解答にはこの方法が書かれている。

$$
{\frac{d}{d t}\left(\frac{\partial L}{\partial \dot q^i}\right)-\frac{\partial L}{\partial q^i}=\sum_{j=1}^{m}\lambda_ｊ{a^j_i}\hspace{20pt}-(6)}
$$

として、拘束に由来する上式右辺の項を付加することで運動方程式が得られるとする。
3.2 拘束条件をとりこんだ拡張ラグランジアンとする方法
近藤[2022]の本文はこの手法で書かれている。近藤[2022]では拡張ラグランジアンを

$$
{\tilde L=L-\sum_{j}{}\lambda_ja^j_k\hspace{20pt}-(7)}
$$

と置き(表記はテキストとは変えている）、これから導かれるオイラー・ラグランジュ方程式と拘束条件を連立させることにより運動方程式が得られるとしている。しかし、この式のままでは期待する結果にならないので、他の文献と見比べ、

$$
{\tilde L=L+\sum_{j}{}\sum_{k}^{}\lambda_ja^j_k\dot q^k \hspace{20pt}-(8)}
$$

と修正して考える。ここからオイラー・ラグランジュ方程式を求めるとラグランジュ乗数$${\lambda}$$は時間の関数なので

$$
{\frac{d}{d t}\left(\frac{\partial L}{\partial \dot q^i}\right)-\frac{\partial L}{\partial q^i}+\sum_{j=1}^{m}\lambda_ｊ\left(\frac{d}{dt}a^j_i-\sum_{k=1}^{n}\frac{\partial a^j_k}{\partial q^i}\dot q^k\right)+\sum_{j=1}^{m}\dot \lambda_ｊa^j_i＝０\hspace{10pt}-(9)}
$$

となる。

3.3  ２つの手法を比較している論文の紹介
Flannery M R [2005]は「拡張ラグランジアンとする方法」と「ダランベールの原理に基づく方法」の両方を取り上げて比較している。二つの手法から導かれる運動方程式は一致しないことがあると書かれている。2021年にもこの問題についての論文は発表されていて、英語圏ではこの問題は今だにホットであり続けているようである。拘束条件つきの最適化問題として見れば、拡張ラグランジアンで書ける場合とそうではない場合を分けるものは何なのかというのは興味深い。
最新の知見を解説するのは私の手に余る、参考論文を見てほしい。

４．計算例

水平面上をすべらずにころがる直立円盤の運動を実際に解いてみる。
座標系は右手系として水平面に$${x,y}$$、鉛直上方に$${z}$$をとり、円盤のころがり方向がx軸の正方向からz軸の右ネジ回転する角度$${\psi}$$,円盤のころがり角を$${\phi}$$とする。円盤の質量を$${m}$$、半径を$${R}$$、転がり軸まわりの慣性モーメントを$${I}$$、垂直軸まわりの慣性モーメントを$${J}$$とする。
　運動方程式の導出については 補に書く。
 数値計算結果

　系の諸元値
　$${m=1,R=1,I=1,J=0.5}$$
 　状態の初期値
　　$${x=1,y=1,\phi=\pi/2,\theta=0,}$$
　　$${\dot x=0,\dot y=1,\dot\phi=1,\dot \theta=1}$$
数値積分法はルンゲクッタ法の４次
時間刻み 0.02 で700点（円軌道を２周回余）計算した結果をプロットした。
数式の操作、数値計算にはmaximaを用い、プロットはgnuplot。
4.1 ダランベールの原理に基づく方法
後述の6.補.1の式(19)~(24)を用いて計算する。x,yの軌道のみあげる。
円盤の重心(x,y)の軌道
　半径１で円を描く。

4.2 拡張ラグランジアンとする手法
後述の6.補.2の式(36)~(41)を使って数値計算する
－１）ダランベールの原理に基づく方法と一致する条件を選んだ場合
　6.補.2に後述するが、数値計算の初期値の選び方で4.1と一致する結果が得られる。しかし、この初期値は系の物理的初期条件としては決まらない。
　ラグランジュ乗数の初期値を$${\mu_x=0,\mu_y=-1}$$とする。
計算誤差の出方が違うが4.3.1と同じ結果である。

－２）別の初期値を選んだ場合
　ラグランジュ乗数の初期値を$${\mu_x=0,\mu_y=0}$$とする
　前二者とは全く違う軌道になっている

５．最後に

　完全に古典力学の範囲内の問題であるが、100年たっても決着がつかずに議論が続いているというのはおもしろい。奥が深い。
　微分幾何に基づく解析力学を勉強中であるが、これが理解できると見え方が変わりそうな気がして楽しみである。

６.補.運動方程式の導出

6.1 ダランベールの原理に基づく方法で解く。
拘束を考慮しないラグランジアンは

$$
{L＝\frac{1}{2}m(\dot x^2+\dot y^2)+\frac{1}{2}I\dot\theta^2+\frac{1}{2}J\dot\varphi^2\hspace{20pt}-(10)}
$$

拘束条件は

$$
{\dot x-R(cos\varphi)\dot \theta=0\hspace{20pt}-(11)} \\{\dot y-R(sin\varphi)\dot \theta=0\hspace{20pt}-(12)}
$$

の二つとなる。拘束条件$${(11)}$$に対するラグランジュ乗数を$${\lambda_x}$$,拘束条件$${12}$$に対し$${\lambda_x}$$とする。
各変数についてのラグランジュ方程式は

$$
{x:\ \ m\ddot x=\lambda_x\hspace{20pt}-(13)}\\
{y:\ \ m\ddot y=\lambda_y\hspace{20pt}-(14)}\\
{\phi:\ \ J\ddot \varphi=0\hspace{20pt}-(15)}\\
{\theta:\ I\ddot \theta=\lambda_xRcos(\varphi)-\lambda_yRsin(\varphi)\hspace{20pt}-(16)}
$$

$${(11,12)}$$を１回微分して

$$
{\ddot x-R(cos\theta)\ddot \theta+R(sin\varphi)\dot\varphi\dot\theta=0\hspace{20pt}-(17)} \\
{\ddot y-R(sin\varphi)\ddot \theta-R(sin\varphi)\dot\varphi\dot\theta=0\hspace{20pt}-(18)}
$$

$${(13)}$$～$${(18)}$$を連立させて解き数値計算用に正規形というか右辺に二次微分とラグランジュ乗数が現れない形にする。

$$
{\ddot x=-R\dot\varphi\dot\theta sin\varphi\hspace{20pt}-(19)}\\
{\ddot y=R\dot\varphi\dot\theta cos\varphi\hspace{20pt}-(20)}\\
{\ddot \varphi=0\hspace{20pt}-(21)}\\
{\ddot \theta=0\hspace{20pt}-(22)}\\
{\lambda_x=-mR\dot\varphi\dot\theta sin\varphi\hspace{20pt}-(23)}\\
{\lambda_y=mR\dot\varphi\dot\theta cos\varphi\hspace{20pt}-(24)}
$$

を得る。重心は等速で円を描いていて円盤の軸が常に円運動の中心を向いているような重心回りの回転をしている。
6.2拡張ラグランジアンとする手法
拘束条件をとりこんだ拡張ラグランジアンをつくる手法で解く。拡張ラグランジアンを$${\hat L}$$、拘束条件に対応するラグランジュ乗数は$${ \mu_x,\mu_y}$$として

$$
{\hat L＝\frac{1}{2}m(\dot x^2+\dot y^2)+\frac{1}{2}I\dot\theta^2+\frac{1}{2}J\dot\varphi^2}\\
{+\mu_x(\dot x-R\dot\theta cos\varphi)+\mu_y(\dot y-R\theta sin\varphi)\hspace{20pt}-(25)}
$$

各変数についてのラグランジュ方程式は

$$
{x:\ \ m\ddot x＋\dot\mu_x=0\hspace{20pt}-(26)}\\
{y:\ \ m\ddot y+\dot\mu_y=0\hspace{20pt}-(27)}\\
{\phi:\ \ J\ddot \varphi+R\dot\theta(-\mu_x sin\varphi+\mu_y cos\varphi)=0\hspace{20pt}-(28)}\\
{\theta:\ \ I\ddot \theta-R\frac{d}{dt}\left(\mu_xcos\varphi+\mu_ysin\varphi\right)=0\hspace{20pt}-(29)}\\
{\mu_x:\ \ \ -\dot x+R(cos\varphi)\dot \theta=0\hspace{20pt}-(30)}\\
{\mu_y:\ \ \ -\dot y+R(sin\varphi)\dot \theta=0\hspace{20pt}-(31)}
$$

$${(30),(31)}$$は拘束条件そのものである。$${(26)(27)}$$を積分すると

$$
{x:\ \ \mu_x=-m\dot x＋C_x\hspace{20pt}-(32)}\\
{y:\ \ \mu_y=-m\dot y+C_y\hspace{20pt}-(33)}\\
$$

$${C_x,C_y}$$は積分定数。これに$${(30)(31)}$$を代入すると

$$
{x:\ \ \mu_x=-mR\dot\theta cos\varphi＋C_x\hspace{20pt}-(32)}\\
{y:\ \ \mu_y=-mR\dot\theta sin\varphi+C_y\hspace{20pt}-(33)}\\
$$

更に$${(32)(33)}$$を$${(30)(31)}$$に代入して整理すると

$$
{J\ddot \varphi=R\dot\theta(C_xsin\varphiーC_ycos\varphi)\hspace{20pt}-(34)}\\
{(I+MR^2)\ddot \theta=R\dot\varphi(-C_xcos\varphi+C_ysin\varphi)\hspace{20pt}-(35)}
$$

$${(21)(22)}$$と$${(34)(35)}$$を比較すると$${C_x=C_y=0}$$であれば一致する。しかし、$${C_x,C_y}$$は初期値からも拘束条件からも決まらないのでどう設定すべきなのか応用観点からは使いづらい解である。数値計算用に正規形を求めておく。

$$
{ \ddot x}=-{ R\dot \varphi}{ \dot \theta}\sin \varphi
-{{R^2\dot \varphi}\cos\varphi\left(-\mu_x\sin
\varphi+\mu_y\cos\varphi\right)
\over{mR^2+I}}{\hspace{10pt}-(36)}\\
{ \ddot y}={ R\dot\varphi} { \dot \theta \cos \varphi} -
{{R^2 \dot \varphi}\sin\varphi\left( -\mu_x \sin \varphi+
\mu_y\cos \varphi \right)
\over{m R^2+I}}{\hspace{10pt}-(37)}\\
{ \ddot \varphi}={{R\dot \theta\left(-\mu_x \sin \varphi+\mu_y \cos \varphi\right) }\over{J}}{\hspace{10pt}-(38)}\\
{ \ddot \theta}={{R\dot\varphi\left(
-\mu_x \sin \varphi+ \mu_y \cos \varphi
\right) }\over{m R^2+I}}{\hspace{10pt}-(39)}\\
{ \ddot \mu_x}={ mR\dot \varphi}{ \dot \theta}\sin \varphi
+{mR^2{ \dot \varphi}\cos\varphi\left(-\mu_x\sin
\varphi+\mu_y\cos\varphi\right)
\over{mR^2+I}}{=\left(m\ddot x \right)\hspace{10pt}-(40)}\\
{ \ddot \mu_y}={ -m R\dot\varphi} { \dot \theta} \cos \varphi +
{{mR^2{ \dot \varphi}\sin\varphi\left( -\mu_x \sin \varphi+\mu_y\cos \varphi 
\right)
}\over{m R^2+I}}{=\left(m\ddot y \right)\hspace{10pt}-(41)}
$$