1.4 计算作用量

为了阐明上述思想，并介绍它们作为计算机程序的表述形式，我们考虑最简单的力学系统——一个在三维空间中运动的自由粒子。欧拉和拉格朗日发现，对于自由粒子，动能沿粒子实际路径的时间积分，小于在同一点之间沿任何替代路径的相同积分：自由粒子按照平稳作用量原理运动，前提是我们将拉格朗日量取为动能。质量为 m、速度为 的粒子的动能为 (1/2) m v2，其中 v 是 的大小。在这种情况下，我们可以选择广义坐标为通常的直角坐标。

遵循欧拉和拉格朗日，自由粒子的拉格朗日量为26

其中形式参数 x 命名了关于给定直角坐标系的位置分量的元组，形式参数 v 命名了速度分量的元组。27

我们可以将这个公式表达为一个过程：

(define ((L-free-particle mass) local)
  (let ((v (velocity local)))
    (* 1/2 mass (dot-product v v))))

该定义表明，L-free-particle 是一个以质量作为参数并返回一个以局部元组 local 为参数的过程，它通过过程 velocity 提取广义速度，并利用该速度计算拉格朗日量的值。

假设我们令 q 表示一个将时间映射到位置分量的坐标路径函数：29

我们可以进行如下定义30

(define q
  (up (literal-function 'x)
      (literal-function 'y)
      (literal-function 'z)))

其中 literal-function 创建了一个过程，该过程表示一个单参数函数，除了给定的符号名称外没有已知的其他性质。符号 q 现在命名了一个单实参（时间）的过程，该过程产生一个三分量元组，表示该时刻的坐标。例如，我们可以对符号时间 t 求值该过程如下：

(print-expression (q 't))
(up (x t) (y t) (z t))

过程 print-expression 生成表达式的可打印形式，并在打印之前对表达式进行化简。

坐标路径 Dq 的导数是那个将时间映射到速度分量的函数：

我们可以创建并使用函数的导数。31 例如，我们可以写出：

(print-expression ((D q) 't))
(up ((D x) t) ((D y) t) ((D z) t))

函数 接受一个坐标路径，并返回一个时间的函数，该函数给出局部元组（ t, q(t), Dq(t), ... ）。我们通过过程 Gamma 来实现这个 。以下是 Gamma 的功能：

(print-expression ((Gamma q) 't))
(up t
    (up (x t) (y t) (z t))
    (up ((D x) t) ((D y) t) ((D z) t)))

因此，复合函数 L o 是一个时间的函数，它返回该路径上该点的拉格朗日量的值：32

(print-expression 
  ((compose (L-free-particle 'm) (Gamma q)) 't))
(+ (* 1/2 m (expt ((D x) t) 2))
   (* 1/2 m (expt ((D y) t) 2))
   (* 1/2 m (expt ((D z) t) 2)))

过程 show-expression 与 print-expression 类似，不同之处在于它将化简后的表达式转换为传统的中缀形式并显示结果。33 我们在本书中使用这种显示方法来生成方框表达式。过程 show-expression 也会生成 print-expression 返回的前缀形式，但我们通常不展示它。34

(show-expression
 ((compose (L-free-particle 'm) (Gamma q)) 't))

根据方程 (1.11)，我们可以计算从时间 t1 到时间 t2 的拉格朗日作用量：

(define (Lagrangian-action L q t1 t2)
  (definite-integral (compose L (Gamma q)) t1 t2))

Lagrangian-action 接受一个计算拉格朗日量的过程 L、一个计算坐标路径的过程 q，以及起始和结束时间 t1 和 t2。这里使用的 definiteintegral 接受一个函数和两个边界 t1 和 t2，并在从 t1 到 t2 的区间上计算该函数的定积分。35 注意，Lagrangian-action 的定义不依赖于任何特定的坐标系，甚至不依赖于构型空间的维数。从拉格朗日量的坐标表示和坐标路径计算作用量的方法不依赖于坐标系。

我们现在可以计算自由粒子沿一条路径的作用量。例如，考虑一个粒子沿直线 t ( 4t + 7, 3t + 5, 2t + 1 ) 匀速运动。36 我们将路径表示为一个过程

(define (test-path t)
  (up (+ (* 4 t) 7)
      (+ (* 3 t) 5)
      (+ (* 2 t) 1)))

对于质量为 3 的粒子，我们得到从 t = 0 到 t = 10 的作用量为37

(Lagrangian-action (L-free-particle 3.0) test-path 0.0 10.0)
435.

练习 1.4 拉格朗日作用量 对于自由粒子，一个合适的拉格朗日量为38

假设 x 是自由粒子的匀速直线路径，使得 xa = x(ta) 且 xb = x(tb)。证明解路径上的作用量为

最小作用量路径

我们已经知道自由粒子的实际路径是匀速直线运动。根据欧拉和拉格朗日，沿直线测试路径的作用量小于沿邻近路径的作用量。设 q 是一条直线测试路径，其作用量为 S[q](t1, t2)。设 q + 是一条邻近路径，由 q 加上一个路径变分 乘以实参数 得到。39 变分路径上的作用量为 S[q + ](t1, t2)。欧拉和拉格朗日发现，对于任何在端点处为零的 和任何非零小量 ，有 S[q + ](t1, t2) > S[q](t1, t2)。

让我们通过数值计算来验证这一点：对测试路径进行变分，添加一定量的在端点 t = t1 和 t = t2 处为零的测试函数。为了构造一个在端点处为零的函数 ，给定一个足够良态的函数 ，我们可以使用 (t) = (t - t1)(t - t2)(t)。这可以实现为：

(define ((make-eta nu t1 t2) t)
  (* (- t t1) (- t t2) (nu t)))

我们可以使用此函数来计算自由粒子在从给定路径变分而来的路径上的作用量，作为 的函数：40

(define ((varied-free-particle-action mass q nu t1 t2) epsilon)
  (let ((eta (make-eta nu t1 t2)))
    (Lagrangian-action (L-free-particle mass)
                       (+ q (* epsilon eta))
                       t1
                       t2)))

变分路径上的作用量，其中 (t) = ( sin t, cos t, t2 ) 且 = 0.001，正如预期的那样，大于测试路径上的作用量：

((varied-free-particle-action 3.0 test-path 
                              (up sin cos square)
                              0.0 10.0)
 0.001)
436.29121428571153

我们可以数值计算作用量取最小值时 的值。我们在例如 -2 到 1 之间搜索：41

(minimize
 (varied-free-particle-action 3.0 test-path 
                              (up sin cos square) 
                              0.0 10.0)
 -2.0 1.0)
(-1.5987211554602254e-14 435.0000000000237 5)

我们得到的结果完全符合预期—— 的最佳值为零，42 并且作用量的最小值就是沿直线路径的作用量。

寻找使作用量最小化的轨迹

我们已经使用变分原理来确定给定的轨迹是否可实现。我们也可以使用变分原理来寻找轨迹。给定一组由有限个参数指定的轨迹，我们可以搜索参数空间，通过寻找使作用量最小化的轨迹，来找到集合中最接近实际轨迹的轨迹。通过选择一组好的逼近函数，我们可以任意接近真实轨迹。43

构造一条具有固定端点的参数化路径的一种方法是使用经过端点以及若干中间点的多项式。中间点位置的变化会改变路径；变分路径的参数就是中间位置的坐标。过程 make-path 使用拉格朗日插值多项式构造这样的路径。44 调用过程 make-path 时需要五个参数：(make-path t0 q0 t1 q1 qs)，其中 q0 和 q1 是端点，t0 和 t1 是对应的时间，qs 是中间点的列表。

指定了一条参数化路径后，我们可以构造一个参数化作用量，它就是沿该参数化路径计算得到的作用量：

(define ((parametric-path-action Lagrangian t0 q0 t1 q1) qs)
  (let ((path (make-path t0 q0 t1 q1 qs)))
    (Lagrangian-action Lagrangian path t0 t1)))

我们可以通过寻找使作用量最小化的参数来找到近似的解路径。我们使用一个现成的多维最小化过程来执行这一最小化：45

(define (find-path Lagrangian t0 q0 t1 q1 n)
  (let ((initial-qs (linear-interpolants q0 q1 n)))
    (let ((minimizing-qs
           (multidimensional-minimize
            (parametric-path-action Lagrangian t0 q0 t1 q1)
            initial-qs)))
      (make-path t0 q0 t1 q1 minimizing-qs))))

过程 multidimensional-minimize 接受一个计算待最小化函数（此处为作用量）的过程（在此是调用 parametricpathaction 的返回值）以及一个参数的初始猜测。这里我们选择初始猜测为两个端点之间直线上的等间距点，该值由 linear-interpolants 计算得到。

为了说明这一策略的应用，我们将寻找谐振子的轨迹，其拉格朗日量为46

其中质量为 m，弹簧常数为 k。该拉格朗日量由以下过程实现：

(define ((L-harmonic m k) local)
  (let ((q (coordinate local)) 
        (v (velocity local)))
    (- (* 1/2 m (square v)) (* 1/2 k (square q)))))

我们可以找到谐振子在 m = 1、k = 1 条件下，从 q(0) = 1 到 q(/2) = 0 的近似路径，如下所示：47

(define q (find-path (L-harmonic 1.0 1.0) 0. 1. :pi/2 0. 3))

我们知道该谐振子（m = 1, k = 1）的轨迹为

其中振幅 A 和相位 由初始条件决定。对于所选的端点条件，解为 q(t) = cos(t)。近似路径应在 0 到 /2 范围内逼近余弦函数。图 1.1 显示了该过程中多项式逼近的误差。使用三个中间点时的最大误差小于 1.7 x 10^-4。正如预期，随着中间点数量的增加，逼近误差减小。对于四个中间点，误差大约改善了 15 倍。

练习 1.5 求解过程 我们可以通过修改过程 parametric-path-action，使其在每次计算作用量时绘制路径，来观察最小化过程的进展。请尝试如下操作：

(define win2 (frame 0. :pi/2 0. 1.2))

(define ((parametric-path-action Lagrangian t0 q0 t1 q1) 
         intermediate-qs)
    (let ((path (make-path t0 q0 t1 q1 intermediate-qs)))
      ;; display path
      (graphics-clear win2)
      (plot-function win2 path t0 t1 (/ (- t1 t0) 100))
      ;; compute action
      (Lagrangian-action Lagrangian path t0 t1)))

(find-path (L-harmonic 1. 1.) 0. 1. :pi/2 0. 2)

练习 1.6 最小化作用量 假设我们试图通过为一个不可能的问题最小化作用量来获得一条路径。例如，假设我们有一个自由粒子，我们对其施加与粒子自由运动不一致的端点速度条件和位置条件。形式体系是否能抵御这种不恰当的攻击？你可能会发现通过编写程序并观察结果会很有启发。

26 这里我们是在进行函数定义。一个定义为任意选定的形式参数指定了函数的值。可以更改形式参数的名称，只要新名称不与定义中的任何其他符号冲突。例如，以下定义指定了完全相同的自由粒子拉格朗日量：

27 拉格朗日量形式上是局部元组的函数，但任何特定的拉格朗日量只依赖于局部元组的一个有限初始段。我们通过显式声明函数所依赖的局部元组初始段中的元素名称来定义局部元组的函数。

28 我们将局部元组表示为一个复合数据结构，其分量包括时间、广义坐标、广义速度以及可能的高阶导数。我们不想被这些结构中分量的打包和解包细节所困扰，因此我们提供了相关的实用工具。构造函数 ->local 接受时间、坐标和速度，并返回表示局部元组的数据结构。选择器 time、coordinate 和 velocity 从局部结构中提取相应的部分。过程 time 与过程 (component 0) 相同，类似地，coordinate 相当于 (component 1)，velocity 相当于 (component 2)。

29 注意。q 定义中的 x 与上面自由粒子拉格朗日量定义中用作形式参数的 x 不同。字母表中只有那么多字母，因此我们不得不重复使用它们。我们将注意在何处赋予符号新的含义。

30 坐标或速度分量的元组通过过程 up 构造。元组 q 的分量 i 为 (ref q i)。所有索引从零开始。单词 up 是为了提醒我们，在数学记法中这些分量由上标索引。还有通过下标索引的 down 分量元组。参见关于记法的附录。

31 函数的导数产生函数。例如，((D cube) 2) => 12 和 ((D cube) 'a) => (* 3 (expt a 2))。

32 在我们的系统中，算术运算符对符号表达式和数值都是通用的；算术过程可以统一地处理数字或表达式。例如，给定过程 (define (cube x) (* x x x))，我们可以得到它对数值 (cube 2) => 8 或文字符号 (cube 'a) => (* a a a) 的值。

33 该显示结果由 TEX 生成。

34 对于非常复杂的表达式，Scheme 的前缀记法通常更好，但化简几乎总是有用的。我们可以将化简和中缀显示的功能分开。我们将在后面看到这样的例子。

35 Scmutils 包含多种数值积分过程。本节中的示例是通过有理函数外推欧拉-麦克劳林公式计算的，相对误差容限为 10^-10。

36 对于实际的物理情况，我们需要为这些量指定单位，但在此示例中我们不做指定。

37 这里我们使用带小数点的数字来指定参数。这迫使表示形式为浮点数，这对于数值计算是高效的。如果要进行符号代数计算，数字必须为精确整数或有理分数，以便表达式能够可靠地化简为最简形式。这种数字的指定不带小数点。

38 速度的平方大小是 · ，即速度与自身的矢量点积。分量结构的平方定义为各分量平方之和，因此我们直接写为 v2 = v · v。

39 注意，我们是在对函数进行算术运算。我们扩展了算术运算，使得两个相同类型（相同定义域和值域）的函数的组合成为相同定义域上的函数，该函数将自变量函数的值域中的值组合起来。例如，如果 f 和 g 是 t 的函数，那么 fg 就是函数 t f(t)g(t)。函数的常数倍是这样一个函数，其值为常数乘以函数对每个自变量的值：cf 是函数 t c f(t)。

40 注意，我们正在对过程进行加法运算。与我们将算术运算扩展到函数类似，算术运算也被扩展到兼容的过程。

41 minimize 的参数包括：一个实现所讨论的单变量函数的过程，以及待搜索区域的下界和上界。Scmutils 包含多种数值最小化方法；这里使用的是 Brent 算法，误差容限为 10^-5。minimize 返回一个包含三个数字的列表：第一个是取得最小值时的自变量，第二个是得到的最小值，第三个是获得该最小值所需的最小化算法迭代次数。

42 是的，-1.5987211554602254e-14 在最小化器要求的容限下即为零。而 435.0000000000237 也可以认为就是之前得到的 435。

43 有很多好的方法可以构造这样一组参数化的逼近轨迹。可以使用样条或高阶插值多项式；可以使用切比雪夫多项式；也可以使用傅里叶分量。选择取决于想要逼近的轨迹类型。

44 以下是实现 make-path 的一种方法：

(define (make-path t0 q0 t1 q1 qs)
  (let ((n (length qs)))
    (let ((ts (linear-interpolants t0 t1 n)))
      (Lagrange-interpolation-function
       (append (list q0) qs (list q1))
       (append (list t0) ts (list t1))))))

过程 linear-interpolants 生成一个元素列表，该列表对前两个参数进行线性插值。我们在这里使用此过程来指定 ts，即 n 个在 t0 和 t1 之间均匀分布的中间时间点，路径在这些时间点处被指定。正在调整的参数 qs 是这些中间时刻的位置。过程 Lagrange-interpolation-function 接受一个值列表和一个时间列表，并生成一个计算通过这些点的拉格朗日插值多项式的过程。

45 这里使用的最小化器是 Nelder-Mead 下山单纯形法。与数值过程通常的情况一样，nelder-mead 过程的接口很复杂，包含大量可选参数，使用户能够有效地控制误差。在本文中，我们通过将 nelder-mead 封装到更友好的 multidimensional-minimize 中来实现特化。不幸的是，你终有一天必须学会处理复杂的数值过程。

46 别担心。我们知道你还不知道为什么这是正确的拉格朗日量。我们将在 1.6 节中讨论这个问题。

47 按照约定，命名常量以冒号开头。名为 :pi 和 :-pi 的常量的含义与其名称所预期的一致。