古典力学における変分法とオイラー＝ラグランジュ方程式
Resumen:
この授業では、変分計算の手法を用いて解析力学のオイラー＝ラグランジュ方程式の導出を再確認し、その結果を基にブラキストクローナ問題の解法への応用を詳述する。
Objetivos de Aprendizaje:
Al concluir esta clase el estudiante será capaz de:
- 理解する ハミルトンの最小作用の原理
- 証明する オイラー＝ラグランジュ方程式
- 解く オイラー＝ラグランジュ方程式を用いてブラキストクローナ問題を
目次：
古典力学における変分法の意義
変分問題の定式化
オイラー＝ラグランジュ方程式
ブラキストクローナ問題
Wolfram アルゴリズムを含む GitHub リポジトリ
古典力学における変分法の意義
ニュートン力学では、変分法を用いることでより効果的に取り組むことができる多くの問題が存在する。このアプローチはラグランジュの方程式およびハミルトンの最小作用の原理にとって基本である。本質的に、この方法はある量を最大化または最小化する軌道を見つけることにある。例えば、二点間の軌道のうち移動距離や移動時間を最小にするものを求めることができる。このアプローチの例がフェルマーの原理であり、これは光が常に通過時間を最小にする経路に沿って進むことを示しており、これから光の屈折に関する スネルの法則 が導かれる。
変分法は古典力学において多くの利点を有する。例えば、対称性を持つ系に対しては解析的な厳密解が得られ、より複雑な系に対しては変分摂動理論を通じて近似解が得られる。また、力を微分方程式の形で表すことが困難な状況では、最小作用の原理が古典力学の問題をより効率的に解く方法を提供する。要するに、変分法はニュートンの法則の代替的な定式化、物理法則の統一、問題解決の効率向上、実験結果予測の精度向上をもたらす基本的なツールである。
変分問題の定式化
変分法では、汎関数の値を極値とする関数 y(x) を求めることに重点を置く。
J(x,y(x))=\displaystyle \int_{x_1}^{x_2} f\left(x,y(x),\frac{dy(x)}{dx}\right)dx,
その最大値または最小値を求めるためである。この式において、汎関数 J は関数 y(x) とその導関数 dy(x)/dx に依存し、積分範囲の上下限は固定されている。積分を極値とするために、関数 y(x) に変分を加え、汎関数の値を極とする関数を求める。例えば、積分が最小値となる場合には、y(x) にどれほど近い関数であっても、その近傍の任意の関数は汎関数の値を増大させる。
「近接関数」の概念を定式化するために、すべての可能な関数 y に対してパラメトリックな表示 y=(\alpha,x) を割り当て、\alpha=0 のとき y(0,x)=y(x) が J を極値とする関数であるとする。このことは次のように表現できる。
y(\alpha, x) = y(x) + \alpha \eta(x),
ここで \eta(x) は \mathcal{C}^1 級の関数で、x_1 と x_2 で消えるものであり、この変分を含む関数 y(\alpha,x) は積分の始点と終点において y(x) と一致する。
この変分 \eta(x) を含む関数 y(\alpha,x) を y(x) の代わりに汎関数 J を定義する積分に代入すると、パラメータ \alpha に依存する新たな汎関数が得られる。
J(x,y(\alpha, x)) = \displaystyle \int_{x_1}^{x_2} f\left(x,y(\alpha,x), \dfrac{d}{dx}y(\alpha,x)\right)dx
局所的な極が存在するためには、次の条件が満たされる必要がある。
\left.\dfrac{\partial J(x,y(\alpha,x))}{\partial \alpha}\right|_{\alpha=0} = 0
任意の関数 \eta(x) に対して。
オイラー＝ラグランジュ方程式
導関数 \partial J(x,y(\alpha,x))/\partial \alpha を解析すると、次のようになる。
\begin{array}{rll} {}\dfrac{\partial J(x,y(\alpha,x))}{\partial \alpha} &=&\dfrac{\partial}{\partial \alpha} \displaystyle \int_{x_1}^{x_2} f\left(x,y(\alpha,x),\dfrac{dy(\alpha, x)}{dx}\right)dx \\ \ &=&\displaystyle \int_{x_1}^{x_2} \left(\dfrac{\partial f}{\partial x}\dfrac{\partial x}{\partial \alpha} + \dfrac{\partial f}{\partial y(\alpha, x)}\dfrac{\partial y(\alpha, x)}{\partial \alpha} + \dfrac{\partial f }{ \partial \frac{dy(\alpha,x)}{dx}} \dfrac{\partial \frac{dy(\alpha,x)}{dx}}{\partial \alpha}\right)dx \end{array}
ここから、次のことに注意することが重要である。
\begin{array}{rll} \dfrac{\partial x}{\partial \alpha} &=& 0 \\ \\ \dfrac{\partial y(\alpha,x)}{\partial \alpha} &=& \dfrac{\partial}{\partial \alpha} \left(y(x) + \alpha \eta(x) \right) = \eta(x) \\ \\ \dfrac{\partial}{\partial \alpha}\left( \dfrac{dy(\alpha, x)}{dx} \right)&=& \dfrac{\partial}{\partial \alpha} \left(\dfrac{dy(x)}{dx} + \alpha\dfrac{d\eta(x)}{dx} \right) = \dfrac{d\eta}{dx} \end{array}
したがって、この式は次のように簡略化される。
\begin{array}{rll} {} \dfrac{\partial J(x,y(\alpha,x))}{\partial \alpha} &=& \displaystyle \int_{x_1}^{x_2} \left(\dfrac{\partial f}{\partial y(\alpha,x)}\eta(x) + \dfrac{\partial f}{\partial \frac{dy(\alpha,x)}{dx}} \dfrac{d\eta(x)}{dx} \right)dx \\ \\ &=&\displaystyle \int_{x_1}^{x_2} \dfrac{\partial f}{\partial y(\alpha,x)}\eta(x) dx + \int_{x_1}^{x_2} \dfrac{\partial f}{\partial \frac{dy(\alpha,x)}{dx}} \dfrac{d\eta(x)}{dx} dx \end{array}
次に、第2の積分を見ると、部分積分を用いて簡略化できることがわかる。
\begin{array}{rll} \displaystyle \int_{x_1}^{x_2} \dfrac{\partial f}{\partial \frac{dy(\alpha,x)}{dx}} \dfrac{d\eta}{dx} dx &=& \left. \dfrac{\partial f}{\partial \frac{dy(\alpha,x)}{dx}} \eta(x)\right|_{x_1}^{x_2} - \displaystyle \int_{x_1}^{x_2}\eta(x) \dfrac{d}{dx}\left( \dfrac{\partial f}{\partial \frac{dy(\alpha, x)}{dx}} \right) dx\\ \\ &=& - \displaystyle \int_{x_1}^{x_2}\eta(x) \dfrac{d}{dx}\left( \dfrac{\partial f}{\partial \frac{dy(\alpha, x)}{dx}} \right)dx \end{array}
よって
\begin{array}{rll} {} \dfrac{\partial J(x,y(\alpha,x))}{\partial \alpha} &=& \displaystyle \int_{x_1}^{x_2} \left[ \eta(x) \dfrac{\partial f}{\partial y(\alpha, x)} - \eta(x) \dfrac{d}{dx}\left( \dfrac{\partial f}{\partial \frac{dy(\alpha,x)}{dx}} \right) \right]dx \\ \\ &=& \displaystyle \int_{x_1}^{x_2} \left[ \dfrac{\partial f}{\partial y(\alpha, x)} - \dfrac{d}{dx}\left( \dfrac{\partial f}{\partial \frac{dy(\alpha,x)}{dx}} \right) \right] \eta(x) dx \end{array}
したがって \left.\dfrac{\partial J (x,y(\alpha, x))}{\partial \alpha}\right|_{\alpha=0} = 0 という条件と、\eta(x) が x_1 と x_2 で消えるという唯一の条件に従う任意の関数であることから、次が得られる。
\dfrac{\partial f}{\partial y(0, x)} - \dfrac{d}{dx}\left( \dfrac{\partial f}{\partial \frac{dy(0,x)}{dx}}\right) = \dfrac{\partial f}{\partial y(x)} - \dfrac{d}{dx}\left( \dfrac{\partial f}{\partial \frac{dy(x)}{dx}}\right) = 0.
最後に、この表現において表記を簡略化すると、いわゆるオイラー＝ラグランジュ方程式に至る。
\boxed{\dfrac{\partial f}{\partial y}= \dfrac{d}{dx}\left( \dfrac{\partial f}{\partial y^\prime} \right)},
これは、汎関数 J が極値を取るための必要条件をはるかに簡潔に表したものである。
ブラキストクローナ問題
問題の定式化
ブラキストクローナ問題は変分計算によって解かれる機械力学の古典的な問題である。状況は次のとおりである。一定の力場の作用下で運動する物体が、初期点 (x_1,y_1) から終点 (x_2,y_2) へ移動するとし、初期点の高さが終点よりも高い場合、粒子が可能な限り短時間で終点に到達するためにはどのような軌跡を辿るべきか、という問いである。
解の定式化
ブラキストクローナ問題を解くには、状況を単純化して考えると役に立つ。そこで、出発点 (x_1, y_1) を座標の原点に固定し、到達点 (x_2,y_2) は原点の右側で \hat{x} 軸の下にあるものとする。
この場合、重力によって生成される -\hat{y} 方向に働く下向きの力場を考え、運動が摩擦なしで行われると仮定する。この文脈では、粒子が出発点と到達点を結ぶさまざまな軌道に従うものとして、それらの中で移動時間を最小にする軌道を求めることを目的とする。
エネルギーの検討
この問題を解くために、重力系のエネルギー保存を利用できる。粒子の質量 m と速度 v に対する運動エネルギー E_{cin}=\frac{1}{2}mv^2 と重力ポテンシャルエネルギー E_{pot,g} の双方を考慮すると、系の全エネルギーは一定である。ポテンシャルエネルギーの基準は原点に置き、E_{pot,g}(y=0)=0 とし、任意の高さ y では E_{pot,g}(y)=mgy となる。
粒子は原点から速度ゼロで出発するので、その全エネルギーは零に等しい。したがって次が成り立つ。
E_{cin} + E_{pot,g}=0
粒子が基準点より下に落下するので、ポテンシャルエネルギーは負であり、運動エネルギーは正である。このようにしてエネルギー保存の式から速度 v を解くと次の結果を得る。
\begin{array}{rl} {} &\dfrac{1}{2}mv^2 + (-mgy) = 0 \\ \\ \vdash &\dfrac{1}{2}mv^2 = mgy \\ \\ \vdash &v^2 = 2gy \\ \\ \vdash &v = \sqrt{2gy} \end{array}
このようにして、任意の高さ y での粒子の速度をその軌道上の位置に応じて計算することができる。
移動時間の検討
速度を求めた後、移動要素 ds=\sqrt{dx^2 + dy^2} を利用して次のように移動時間の微小要素を構築することができる。
\begin{array}{rl} {} dt &= \dfrac{ds}{v} = \dfrac{\sqrt{dx^2 + dy^2}}{\sqrt{2gy}}\\ \\ &= \sqrt{\dfrac{dx^2 + dy^2}{2gy} } \end{array}
したがって、(x_1,y_1) から (x_2,y_2) までの移動時間は次の積分によって求められる。
\begin{array}{rl} {} t &= \displaystyle \int_{(x_1,y_1)}^{(x_2,y_2)} dt \\ \\ &= \displaystyle \int_{(x_1,y_1)}^{(x_2,y_2)} \sqrt{\dfrac{dx^2 + dy^2}{2gy}} \\ \\ &= \displaystyle \dfrac{1}{\sqrt{2g}}\int_{y_1}^{y_2} \sqrt{\dfrac{1+ \left(\dfrac{dx}{dy}\right)^2 }{y}}dy \\ \\ \end{array}
変分問題の定式化
この最後の式により、時間を次の形の汎関数として表すことができた。
{}t = J(y,x(y)) = \displaystyle \int_{y_1}^{y_2} f\left(y,x(y),\dfrac{dx(y)}{dy} \right) dy
ここで
f\left(y,x(y), \dfrac{dx(y)}{dy}\right) = \sqrt{\dfrac{1+ \left(\dfrac{dx(y)}{dy} \right)^2}{y}}
ここで \sqrt{2g} という因子は無視できる。というのも、J を最適化することは \sqrt{2g}J を最適化するのと全く同じだからである。
以上を踏まえ、先ほどと同じ手順に従ってオイラー＝ラグランジュ方程式を構成すると、最終的に次の式に到達する。
\dfrac{\partial f}{\partial x} = \dfrac{d}{dy} \dfrac{\partial f}{\partial x^\prime}
しかし、ここでは \dfrac{\partial f}{\partial x} = 0 であることがわかるので、
\dfrac{d}{dy}\dfrac{\partial f}{\partial x^\prime} = 0,
言い換えると
\dfrac{\partial f}{\partial x^\prime} = \dfrac{1}{\sqrt{2a}},
ここで a は任意の定数であり、後の展開にとって「便利」であることからこのように記した。
変分問題の解法
この最後の式に f を代入すると、次のようになる。
\begin{array}{rl} {} &\dfrac{\partial }{\partial x^\prime} \sqrt{\dfrac{1+ x^{\prime 2}}{y}} = \dfrac{1}{\sqrt{2a}} \\ \\ \vdash & \dfrac{1}{2}\left( \dfrac{1 + x^{\prime 2} }{y} \right)^{-1/2} \left(\dfrac{2x^\prime}{y} \right) = \dfrac{1}{\sqrt{2a}} \\ \\ \vdash & \dfrac{1}{2}\sqrt{\dfrac{y}{1 + x^{\prime 2}}} \left(\dfrac{2x^\prime}{y} \right) = \dfrac{1}{\sqrt{2a}} \\ \\ \vdash & \sqrt{\dfrac{4x^{\prime 2} y}{4y^2 (1 + x^{\prime 2})} } = \sqrt{\dfrac{1}{2a}} \\ \\ \vdash & \dfrac{y x^{\prime 2} }{y^2 (1 + x^{\prime 2})} = \dfrac{1}{ 2a} \\ \\ \vdash & 2ayx^{\prime 2} = y^2 + y^2 x^{\prime 2} \\ \\ \vdash & x^{\prime 2} (2ay - y^2) = y^2 \\ \\ \vdash & \left(\dfrac{dx}{dy}\right)^2 = \dfrac{y^2}{2ay - y^2} \\ \\ \vdash & \dfrac{dx}{dy} = \pm \sqrt{\dfrac{y^2}{2ay - y^2}} \\ \\ \vdash & dx = \pm \dfrac{y\,dy}{\sqrt{2ay - y^2}} \\ \\ \vdash & x = \displaystyle \pm \int \dfrac{y}{\sqrt{2ay - y^2}}\,dy \end{array}
この積分を解くには、次の置換を行う方法を検討できる。
\begin{array}{rl} {} y &= a[1-\cos(\theta)] \\ dy &= a\sin(\theta) d\theta \end{array}
これにより次のようになる。
\begin{array}{rl} {} x= & \pm \displaystyle \int \dfrac{y}{\sqrt{2ay - y^2}}\,dy = \displaystyle \int \dfrac{a[1-\cos(\theta)]a\sin(\theta)}{\sqrt{2a^2[1-\cos(\theta)] - a^2[1-\cos(\theta)]^2 }}\,d\theta \\ \\ & = \pm \displaystyle \int \dfrac{a^2[1-\cos(\theta)]\sin(\theta)}{\sqrt{a^2[1-\cos(\theta)]\left\{ 2 - [1-\cos(\theta)] \right\} }}\,d\theta \\ \\ & = \pm \displaystyle \int \dfrac{a[1-\cos(\theta)]\sin(\theta)}{\sqrt{[1-\cos(\theta)] [1 + \cos(\theta)] }}\,d\theta \\ \\ & = \pm \displaystyle \int \dfrac{a[1-\cos(\theta)]\sin(\theta)}{\sqrt{ 1-\cos^2(\theta)}}\,d\theta \\ \\ & = \pm \displaystyle \int \dfrac{a[1-\cos(\theta)]\sin(\theta)}{\sin(\theta)}\,d\theta \\ \\ & = \pm \displaystyle \int a[1-\cos(\theta)]\,d\theta \\ \\ & = \pm a(\theta - \sin(\theta)) + C \end{array}
ブラキストクローナ曲線は極座標におけるパラメトリック曲線として表すことができ、原点を出発点とするサイクロイドと一致することがわかる。
\begin{array}{rl} {} x(\theta) &= \pm a(\theta - \sin(\theta)) \\ y(\theta) &= a(1-\cos(\theta)) \end{array}
積分定数 C は、軌道が原点から始まるという初期条件を満たすために消去された。また、定数 a を調整して曲線が最後に (x_2,y_2) を通るようにすることができる二つの方程式が存在することがわかる。これらの方程式は次のとおりである。
選択肢 1： \boxed{\begin{array}{rl} {} x(\theta) &= a(\theta - \sin(\theta)) \\ y(\theta) &= a(1-\cos(\theta)) \end{array}}
選択肢 2： \boxed{\begin{array}{rl} {} x(\theta) &= - a(\theta - \sin(\theta)) \\ y(\theta) &= a(1-\cos(\theta)) \end{array}}
この問題に対する実行可能な解は第二の選択肢によって与えられ、定数 a を負の値に設定することにより、解として必要な条件を満たす曲線が得られる。
解の最終調整
最後の調整を施した後、ブラキストクローナ曲線は次のパラメトリック形式をとる。
\begin{array}{rl} x(\theta) &= b(\theta - \sin(\theta)) \\ y(\theta) &= -b(1-\cos(\theta)) \end{array}
ここでは a=-b と置き換え、0 \lt b とした。この曲線の周期は 2b\pi であり、x_2 \in ]0,2b\pi[[ および y_2 \in ]-2b,0[[ という条件を満たさなければならない。これは極めて重要であり、ブラキストクローナ曲線がサイクロイドの1本のアークとして表されなければならないことを意味する。というのも、粒子が高さゼロの点に戻って安定すると解が有効でなくなるからである。
この方程式を問題に適用するには、次の連立方程式を満たす \theta と b の値を求める必要がある。
\begin{array}{rl} {} x_2 &= b(\theta - \sin(\theta))\\ y_2 &= - b(1-\cos(\theta)) \end{array}
この非線形システムは解析的な解を持たないように見えるので、Wolfram Mathematica による数値的手法を用いることにする。以下に問題を解くための手順を示す。
ステップ 1: システムの設定
解くべきシステムを構成する方程式を設定する
eq1 = x2 == b*(theta - Sin[theta])
eq2 = y2 == -b*(1 - Cos[theta])
ステップ 2: 到達点を定義する
粒子が軌道の終端で到達する点を設定する。ここでは (x_2,y_2)=(1,-2) とする。これらの値を変更して同様の別の設定を試すことができる。
x2val = 1; y2val = -2;
ステップ 3: 求める値を数値的に計算する
問題の解を数値的に計算するために 「FindRoot」 関数を使用する
sol = FindRoot[{eq1, eq2} /. {x2 -> x2val, y2 -> y2val}, {{b,1}, {theta, 1}}]
ここでは数値解の初期近似として b=1 および \theta=1 を使用した。その結果、b\approx 2.4056 および \theta \approx 1.40138 という解が得られる。
ステップ 4: 結果の検証
これらの解が物理的に意味を持つためには、x_2 \in ]0,2b\pi[[ および y_2 \in ]-2b, 0[[ である必要があることを思い出そう。次の手順によりこれを迅速に確認できる。
まず、求めた解から b と \theta の値を抽出する。
bval = sol[[1, 2]]; thetaval = sol[[2, 2]];
次に、確認処理を実行するよう指示する。
If[0 < x2val < 2*Pi*bval && -2*bval < y2val < 0, "Valores válidos", "Valores inválidos"]
うまくいけば出力は「Valores válidos」となるはずである。このコード片は、物理状況が正しくモデル化されているかどうかを確認するのに役立つ。
これらの手順により、出発点 (x_1,y_1)=(0,0) と到達点 (x_2,y_2)=(1,-2) を結ぶ解曲線を完全に調整することができた。得られた曲線は次の通りである。
\begin{array}{rl} {} x(\theta) &\approx 2.4056(\theta - \sin(\theta)) \\ y(\theta) &\approx -2.4056(1-\cos(\theta)) \end{array}\;\;;\theta\in [0, 1.40138]
そのグラフは次のように見える。
Wolfram アルゴリズムを含む GitHub リポジトリ
ブラキストクローナ問題の解決に用いた完全なコード（Wolfram Mathematica で開発されたアルゴリズムを含む）は、私の GitHub リポジトリからダウンロードおよび閲覧が可能である。このリポジトリにはインタラクティブなノートブック形式のコードが格納された .nb ファイルと、コードを直接見たい人向けのプレーンテキスト版である .m ファイルが含まれている。
GitHub からリポジトリを こちら からダウンロードできる。
コードに加えて、このリポジトリにはアルゴリズムの使用方法や理解方法に関する詳細な説明を含む「README」ファイルと、ブラキストクローナ問題の解法を段階的に解説した資料が含まれている。役に立てば幸いである。