1.7 动力学状态的演化

拉格朗日方程是路径必须满足的常微分方程。它们可用于检验一条假想路径是否为系统的可实现路径。然而，我们也可以利用它们从初始条件出发来逐步构建一条路径。

系统的“状态”被定义为必须加以规定、以便后续演化得以确定的信息。回想一下我们的杂耍者：他必须以某种方式抛出球棒，才能使其完成期望的运动。杂耍者可以控制球棒的初始位置和取向，以及初始速度和旋转。杂耍以及类似系统的经验表明，初始构型和构型的变化率足以确定后续运动。其他系统可能需要构型的更高阶导数。

对于用一组广义坐标和速度表示的拉格朗日量，我们已经证明拉格朗日方程是二阶常微分方程。如果这些微分方程可以解出最高阶导数，并且满足适当的条件，[73] 那么初值问题存在唯一解：即在给定时刻给出解及其低阶导数的值，存在唯一的解函数。给定无冗余坐标，拉格朗日方程满足这些条件。[74] 因此，轨迹由任意时刻的广义坐标和广义速度所确定。这就是指定动力学状态所需的信息。

路径的完整局部描述由该路径在某一时刻的所有导数组成。给定从局部元组的初始片段出发、根据低阶导数计算高阶导数的规则，就可以重构路径的完整局部描述。系统的状态由局部元组中能够推断出完整局部描述的那一初始片段来指定。完整局部描述给出了路径在该时刻附近的行为。实际上，我们所需要的只是一个计算下一个更高阶导数的规则；由此我们可以得到所有其余导数。假设系统的状态由元组 ( t, q, v ) 给出。如果给定一个计算加速度的规则 a = A(t, q, v)，那么

因此我们有

以此类推。因此局部元组的高阶导数分量由函数 Dt A、Dt2 A、... 给出。这些函数中的每一个都依赖于局部元组的低阶导数分量。从状态推导路径所需要的全部就是一个函数，它从状态给出局部描述的下一更高阶导数分量。我们利用拉格朗日方程来寻找这个函数。

首先，我们展开拉格朗日方程

使得二阶导数显式出现：

求解此方程组的 D2q，可得到沿解路径 q 的广义加速度：

其中 [ 2 2 L o ] 是一个可用对称方阵表示的结构，因此我们可以计算其逆矩阵。给出加速度的函数是

其中 = I2 是速度分量的选择器。

局部元组中用于指定状态的那一初始片段称为局部状态元组，或更简单地称为状态元组。

我们可以将给出加速度（作为状态元组的函数）的函数表示为以下过程。该过程接受一个计算拉格朗日量的过程，并返回一个以状态元组为参数并返回加速度的过程。[75]

(define (Lagrangian->acceleration L)
  (let ((P ((partial 2) L))
        (F ((partial 1) L)))
    (/ (- F
          (+ ((partial 0) P) 
             (* ((partial 1) P) velocity)))
       ((partial 2) P))))

一旦我们有了从坐标和速度计算加速度的方法，就可以给出从状态计算状态导数的规则。对于时刻 t 的状态 ( t, q(t), Dq(t) )，状态的导数为 ( 1, Dq(t), D2q(t) ) = ( 1, Dq(t), A(t, q(t), Dq(t)) )。过程 Lagrangian->state-derivative 接受一个拉格朗日量，返回一个接受状态并返回状态导数的过程：

(define (Lagrangian->state-derivative L)
  (let ((acceleration (Lagrangian->acceleration L)))
    (lambda (state)
      (up 1
          (velocity state)
          (acceleration state)))))

我们用状态元组的 up-元组来表示状态，该元组由确定状态的局部元组初始片段的各分量组成。

例如，谐振子的参数化状态导数为

(define (harmonic-state-derivative m k)
  (Lagrangian->state-derivative (L-harmonic m k)))

(print-expression
 ((harmonic-state-derivative 'm 'k)
  (up 't (up 'x 'y) (up 'v_x 'v_y))))
(up 1 (up v_x v_y) (up (/ (* -1 k x) m) (/ (* -1 k y) m)))      

拉格朗日方程是约束可实现路径 q 的二阶微分方程组。我们可以利用状态导数将拉格朗日方程表示为约束可实现坐标路径 q 和速度路径 v 的一阶微分方程组：

(define ((Lagrange-equations-first-order L) q v)
  (let ((state-path (qv->state-path q v)))
    (- (D state-path)
       (compose (Lagrangian->state-derivative L)
                state-path))))

(define ((qv->state-path q v) t)
  (up t (q t) (v t)))

例如，我们可以求出二维谐振子运动方程的一阶形式：

(show-expression
 (((Lagrange-equations-first-order (L-harmonic 'm 'k))
   (up (literal-function 'x)
       (literal-function 'y))
   (up (literal-function 'v_x)
       (literal-function 'v_y)))
  't))

拉格朗日方程残差结构第一个元素中的零只是时间均匀流逝的同义反复：时间函数就是恒等函数，因此其导数为 1，残差为零。第二个元素中的方程约束速度路径为坐标路径的导数。第三个元素中的方程根据施加的力给出了速度的变化率。

数值积分

给出状态关于状态导数的函数的一阶常微分方程组，可以通过积分从给定的初态出发找到状态路径。数值积分器通过图 1.6 所示的过程来寻找这些微分方程的近似解。Lagrangian->state-derivative 所产生的状态导数可被用于数值积分一阶常微分方程组的软件包。

过程 state-advancer 可用于在指定时刻找到系统的状态，给定包括初始时间在内的初态，以及一个参数化状态导数过程。[76] 例如，为了推进二维谐振子的状态，我们编写[77]

(print-expression
  ((state-advancer harmonic-state-derivative 2 1)
   (up 0 (up 1 2) (up 3 4))
   10
   1.e-12))
(up 10. 
    (up 3.712791664584467 5.420620823651575)
    (up 1.6148030925459906 1.8189103724750977))

state-advancer 的参数有：参数化状态导数 harmonic-state-derivative 和状态导数参数（质量 2 和弹簧常数 1）。返回的过程接受初态 (up 0 (up 1 2) (up 3 4))、目标时间 10 和相对误差容限 1.e-12。输出是指定最终时刻状态的近似值。

考虑上述具有周期性驱动的受驱摆。我们选择 ys(t) = a cos t。

(define ((periodic-drive amplitude frequency phase) t)
  (* amplitude (cos (+ (* frequency t) phase))))
(define (L-periodically-driven-pendulum m l g a omega)
  (let ((ys (periodic-drive a omega 0)))
    (L-pend m l g ys)))

该系统的拉格朗日方程为

(show-expression
 (((Lagrange-equations
    (L-periodically-driven-pendulum 'm 'l 'g 'a 'omega))       
   (literal-function 'theta))
  't))

周期性受驱摆的参数化状态导数为

(define (pend-state-derivative m l g a omega)
  (Lagrangian->state-derivative
    (L-periodically-driven-pendulum m l g a omega)))

(show-expression
  ((pend-state-derivative 'm 'l 'g 'a 'omega)
   (up 't 'theta 'thetadot)))

为了研究受驱摆的演化，我们需要一种机制，能够在系统演化的某个时间区间内监控系统的各个方面。过程 evolve 提供了这一功能，它重复使用 state-advancer 将状态推进到所需的时刻。过程 evolve 接受一个参数化状态导数及其参数，返回一个过程，该过程从指定的初态开始将系统演化到多个其他时刻，并在这些时刻监控状态的某个方面。为了生成角度随时间变化的图形，我们制作一个监控过程，在演化过程中生成图形：[78]

(define ((monitor-theta win) state)
  (let ((theta ((principal-value :pi) (coordinate state))))
    (plot-point win (time state) theta)))

(define plot-win (frame 0. 100. :-pi :pi))
((evolve pend-state-derivative 
         1.0                   ;m=1kg
         1.0                   ;l=1m
         9.8                   ;g=9.8m/s²
         0.1                   ;a=1/10 m
         (* 2.0 (sqrt 9.8)) )  ;omega
 (up 0.0                       ;t₀=0
     1.                        ;theta₀=1 radian
     0.)                       ;thetadot₀=0 radians/s
 (monitor-theta plot-win)
 0.01                          ;step between plotted points
 100.0                         ;final time
 1.0e-13)                      ;local error tolerance

图 1.7 显示了受驱摆在两条轨道上的角度 随时间的变化。两次运行的初始条件完全相同，只是其中一次给摆锤一个微小的初始速度，大小为 10-10 m/s，大约为每秒钟一个原子宽度。两条轨道的初始段难以区分。大约 75 秒后，两条轨道开始偏离并变得完全不同。这种对初始条件微小变化的极度敏感性是所谓混沌行为的特征。稍后，我们将使用其他工具（如李雅普诺夫指数、相空间和庞加莱截面）进一步研究这个例子。