经典力学中的变分法与欧拉-拉格朗日方程
摘要：在本节课中，我们将通过变分法技术推导出分析力学中的欧拉-拉格朗日方程，并详细展示其在解决最速降线问题中的应用。
学习目标：
完成本节课后，学生将能够：
- 理解哈密顿的最小作用原理
- 推导欧拉-拉格朗日方程
- 解决使用欧拉-拉格朗日方程解决最速降线问题。
内容目录:
经典力学中为何使用变分法
变分问题的表述
欧拉-拉格朗日方程
最速降线问题
包含 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) 是一个在 x_1 和 x_2 处为零的 \mathcal{C}^1 类函数，因此包括该变化的函数 y(\alpha,x) 在积分路径的起点和终点与 y(x) 相同。
通过将包括变化 \eta(x) 的函数 y(\alpha,x) 替换进泛函 J 定义的积分中，而不是 y(x)，得到一个新的依赖于参数 \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} {} \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}
然后，如果我们观察第二个积分，将会发现可以通过分部积分来简化：
\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}），并假设运动在没有摩擦的情况下进行。在这种背景下，粒子被限制在不同的路径上，这些路径连接出发点和到达点，目的是找到哪条路径可以最小化旅行时间。
能量分析
为了解决这个问题，我们可以利用系统的能量守恒。系统的总能量保持不变，考虑到动能 E_{cin}=\frac{1}{2}mv^2 和重力势能 E_{pot,g}，其中 m 是粒子的质量，v 是其速度。对于势能，以原点为参考，因此 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)}{dx}\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{ydy}{\sqrt{2ay - y^2}} \\ \\ \vdash & x = \displaystyle \pm \int \dfrac{y}{\sqrt{2ay - y^2}}dy \end{array}
为了解决这个积分，可以考虑进行以下替换：
\begin{array} {} 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} {} & 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} {} & x(\theta) &=& a(\theta - \sin(\theta)) \\ & y(\theta) &=& a(1-\cos(\theta)) \end{array}}
选项 2:\boxed{\begin{array} {} & x(\theta) &=& - a(\theta - \sin(\theta)) \\ & y(\theta) &=& a(1-\cos(\theta)) \end{array}}
此问题的可行解由第二个选项给出，通过将常数 a 设为负值，我们得到一条符合成为解的必要条件的曲线。
解的最终调整
经过最后的调整，最速降线曲线具有以下参数形式：
\begin{array} {} 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[. 这最后一点很关键，因为它要求最速降线曲线表示为单摆线弧，如果粒子在返回零高度点时停止，解将不再有效。
为了将这些方程调整到问题中，我们需要找到满足以下系统的 \theta 和 b 的值：
\begin{array} {} 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 "有效值", "无效值"]
如果一切顺利，我们应该在输出中得到 "有效值"。这段代码将帮助你检查物理情况是否正确建模。
通过这些步骤，我们最终完全调整了我们的解曲线，它连接点 (x_1,y_1)=(0,0) 和 (x_2,y_2)=(1,-2). 最终曲线为：
\begin{array} {} 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"文件，其中详细说明了如何使用和理解算法，并提供了逐步解决最速降线问题的解释。希望你觉得它有用！