基于<bold>Lagrange</bold>插值多项式拟合的力学系统的变分积分子

孔新雷，杨雪; Kong Xinlei; Yang Xue

网刊加载中。。。

使用Chrome浏览器效果最佳，继续浏览，你可能不会看到最佳的展示效果，

确定继续浏览么?

复制成功，请在其他浏览器进行阅读

摘要

变分积分子是通过直接离散变分原理得到的一类特殊的动力学系统的数值差分格式，较之传统差分格式呈现出明显的计算优越性.由离散Euler-Lagrange方程的形式可知，变分积分子的构造过程最终归结为计算离散Lagrange函数的偏导数，其中离散Lagrange函数是Lagrange函数在单个时间步长的积分，通常由经典求积公式近似得到.根据离散Lagrange函数的积分表达式，解析计算其偏导数会随之衍生一个新的且与连续Euler-Lagrange方程密切关联的积分，因此，构造变分积分子就可以不再以通过经典求积公式得到的具体形式的离散Lagrange函数为前提，而是可以直接基于一组离散结点近似新衍生的积分.在这些离散结点处，如果进一步让系统的拟合轨迹严格满足Euler-Lagrange方程，即运动方程，那么新的积分自动为零，相应地，计算离散Lagrange函数的偏导数就简化为计算连续Lagrange函数关于速度变量的偏导数.这种新的构造方式同时结合了连续和离散的Euler-Lagrange方程，不仅让最终得到的差分格式仍然继承了变分积分子特有的优越计算性能，而且在同阶精度的情况下具有更小的局部误差.

关键词

Euler-Lagrange方程; 变分积分子; Lagrange插值多项式

引言

Lagrange系统是由Euler-Lagrange方程

\frac{\partial L}{\partial x_{i}} - \frac{d}{d t} \frac{\partial L}{\partial {\dot{x}}_{i}} = 0, i = 1,2, \dots, n

（1）

所描述的一类力学系统，其中 $x = (x_{1}, x_{1}, \dots, x_{n})$ 代表位形流形 $M$ 的局部坐标， $L (x, \dot{x})$ 是系统的Lagrange函数，一般为系统的动能和势能之差.如果进一步考虑非保守力，那么Euler-Lagrange方程的形式会发生些许变化^［

1］，但这种情况在本文中暂不研究.当系统（1）是一组非线性微分方程时，对其解析求解是非常困难甚至是不可能的，因此需要借助于数值手段，通常是采用经典差分格式，用差分方程代替微分方程实现数值求解.

然而，人为地采用经典差分格式离散系统的运动方程，不可避免地会破坏连续系统具有的几何结构，引入不必要的数值耗散.因此，在“数值算法应尽可能地保持原问题的本质特征”^［

2］的指导原则下，针对系统（1），人们提出并设计了保结构算法^{［参考文献 2

百度学术}2，3］.较早的保结构算法主要是保辛差分格式.借助于Legendre变换，Euler-Lagrange方程可以转化为Hamilton正则方程，而基于Hamilton系统的辛结构，利用生成函数法可以构造系统的保辛差分格式^{［参考文献 2

百度学术}2］.变分积分子是另外一种构造思路的保结构算法.考虑到Euler-Lagrange方程可以由Hamilton原理诱导，仿照这一连续情形，在对系统（1）构造数值算法时可以先离散Hamilton原理，由离散后的Hamilton原理可以导出离散Euler-Lagrange方程.作为离散变分原理的产物，离散Euler-Lagrange方程在作为数值差分格式时，真实地继承了连续系统具有的几何特性，不仅是自然保辛的，而且满足Noether定理进而保持系统的动量^{［参考文献 4

百度学术}4，5］.自然简便的构造方式和优越显著的计算性能，使得变分积分子迅速发展，不仅局限于系统（1），对于非保守系统^{［参考文献 5

百度学术}5，6］、非光滑系统^{［参考文献 7

百度学术}7］、完整和非完整约束系统^{［参考文献 8-10}8-10］、Birkhoff系统^{［参考文献 11-13}11-13］、连续介质力学^{［参考文献 14

百度学术}14，15］、随机系统^{［参考文献 16

百度学术}16，17］等都可以构造系统的变分积分子，进行高效的数值模拟仿真.

本文结合局部路径拟合^［

18，19］方法提出一种新的变分积分子的设计构造方式.考虑到变分积分子的经典构造过程最终归结为计算离散Lagrange函数的偏导数，本文先对相关涉及到的偏导数进行解析计算，发现一个新的依赖于Euler-Lagrange方程的积分随之产生.在对该积分实现科学有效的近似之后，离散Lagrange函数的偏导数就转化为连续Lagrange函数关于速度变量的偏导数，这样就省略了利用经典求积公式计算离散Lagrange函数的过程.因此，这种新的构造方式不仅让最终得到的差分格式仍然继承了变分积分子特有的优越计算性能，而且显著地简化了其构造过程，特别是对于高阶精度的变分积分子.

1 变分积分子

如引言中所述，Lagrange力学系统是由Euler-Lagrange方程

\frac{\partial L}{\partial x_{i}} - \frac{d}{d t} \frac{\partial L}{\partial {\dot{x}}_{i}} = 0, i = 1,2, \dots, n

所描述的一类力学系统.根据Hamilton原理，Euler-Lagrange方程的解 $x = x (t)$ 恰好对应作用泛函 $\int_{0}^{T} L (x, \dot{x}) d t$ 的驻点，即对于任意满足 $δ x (0) = δ x (T) = 0$ 的变分 $δ x (t)$ ， $δ \int_{0}^{T} L (x, \dot{x}) d t = 0$ 等价于 $\frac{\partial L}{\partial x} - \frac{d}{d t} \frac{\partial L}{\partial \dot{x}} = 0$ .这一事实是构造Lagrange系统变分积分子的主要依据.

变分积分子是一类特殊的数值差分格式.与经典差分格式不同，构造力学系统的变分积分子，出发点不是直接离散运动方程，而是离散能够诱导运动方程的变分原理.

具体到Lagrange力学系统，首先考虑离散Lagrange函数

L_{d} (x^{k}, x^{k + 1}) = \int_{t^{k}}^{t^{k + 1}} L (x, \dot{x}) d t

（2）

其中， $x^{k}$ 是 $x (t^{k})$ 的近似， $τ = t^{k + 1} - t^{k}$ 是选定的时间步长.这样

\int_{0}^{T} L (x, \dot{x}) d t = L_{d} (x^{k}, x^{k + 1})

其中 $N = T / τ$ 反映了离散划分结点的疏密程度.仿照连续情形的Hamilton原理，对离散作用泛函 $\sum_{k = 0}^{N - 1} L_{d} (x^{k}, x^{k + 1})$ 做变分并令其等于零，即

δ \sum_{k = 0}^{N - 1} L_{d} (x^{k}, x^{k + 1}) = 0

进一步注意到 $δ x^{0} = δ x (0) = 0$ 以及 $δ x^{N} = δ x (T) = 0$ ，可得离散Euler-Lagrange方程

\partial_{1} L_{d} (x^{k}, x^{k + 1}) + \partial_{2} L_{d} (x^{k - 1}, x^{k}) = 0

（3）

写成分量形式为

\frac{\partial L_{d} (x^{k}, x^{k + 1})}{\partial x_{i}^{k}} + \frac{\partial L_{d} (x^{k - 1}, x^{k})}{\partial x_{i}^{k}} = 0, i = 1,2, \dots, n

当 $\partial_{2} L_{d}$ 满足非退化条件时，离散Euler-Lagrange方程（3）就确定了差分格式

Φ : (x^{k - 1}, x^{k}) \to (x^{k}, x^{k + 1})

这种由离散变分原理所得到的数值算法 $Φ$ 就被称作变分积分子.独特的构造方式避免了人为的离散运动方程所带来的数值耗散，使得变分积分子能够真实地继承连续系统的解所具有的几何性质，例如保辛性，即算法 $Φ$ 满足 $Φ^{*} Ω = Ω$ ，其中微分形式

Ω = \sum_{i, j = 1}^{N - 1} \frac{\partial^{2} L_{d} (x^{k}, x^{k + 1})}{\partial x_{i}^{k} \partial x_{j}^{k + 1}} d x_{i}^{k} \land d x_{j}^{k + 1}

另外，与连续情形类似，当离散Lagrange函数 $L_{d}$ 进一步满足李群作用不变性时，算法 $Φ$ 还保持动量（映射）不变.根据后误差分析理论，变分积分子所表现出的优越计算性能，如高精度、长时间稳定、保持系统守恒量等，都与其保结构特性有着直接关系.当力学系统受到若干完整约束，即位形坐标 $x = (x_{1}, x_{1}, \dots, x_{n})$ 满足代数方程组

f (x) = {(f_{1} (x), f_{2} (x), \dots, f_{p} (x))}^{T} = 0

时，借鉴上述先离散后变分的思想，并结合Lagrange乘子法，同样可以构造完整约束系统

\begin{array}{l} \frac{\partial L}{\partial x} - \frac{d}{d t} \frac{\partial L}{\partial \dot{x}} = \sum_{j = 1}^{p} λ_{j} \frac{\partial f_{j}}{\partial x}, \\ f (x) = 0 \end{array}

（4）

的变分积分子

\begin{array}{l} \partial_{1} L_{d} (x^{k}, x^{k + 1}) + \partial_{2} L_{d} (x^{k - 1}, x^{k}) = <λ^{k}, \nabla f (x^{k})>, \\ f (x^{k}) = 0, \end{array}

（5）

其中，

<λ^{k}, \nabla f (x^{k})> = \sum_{i = 1}^{p} λ_{j}^{k} \frac{\partial f_{j} (x^{k})}{\partial x^{k}} τ

值得注意的是，如果引入扩展的离散Lagrange函数

\begin{array}{l} {\bar{L}}_{d} ({\bar{x}}^{k}, {\bar{x}}^{k + 1}) = {\bar{L}}_{d} (x^{k}, λ^{k}, x^{k + 1}, λ^{k + 1}) \\ = L_{d} (x^{k}, x^{k + 1}) - <λ^{k}, f (x^{k})> \end{array}

或

\begin{array}{l} {\bar{L}}_{d} ({\bar{x}}^{k}, {\bar{x}}^{k + 1}) = {\bar{L}}_{d} (x^{k}, λ^{k}, x^{k + 1}, λ^{k + 1}) \\ = L_{d} (x^{k}, x^{k + 1}) - <λ^{k + 1}, f (x^{k + 1})> \end{array}

那么，方程组（5）与以 ${\bar{L}}_{d} ({\bar{x}}^{k}, {\bar{x}}^{k + 1})$ 为离散Lagrange函数的离散Euler-Lagrange方程

\partial_{1} {\bar{L}}_{d} ({\bar{x}}^{k}, {\bar{x}}^{k + 1}) + \partial_{2} {\bar{L}}_{d} ({\bar{x}}^{k - 1}, {\bar{x}}^{k}) = 0

等价^［

5］.

2 基于Lagrange插值多项式拟合的力学系统的变分积分子

根据离散Euler-Lagrange方程（3）的形式可知，构造力学系统的变分积分子最后归结为计算离散Lagrange函数的偏导数—— $\partial_{1} L_{d}$ 和 $\partial_{2} L_{d}$ .为了计算这些偏导数，首先引入自由变量 $μ$ ，在实际计算偏导数时 $μ$ 可以代表 $x_{i}^{k}$ 或 $x_{i}^{k + 1}$ .由于离散Lagrange函数

L_{d} (x^{k}, x^{k + 1}) = \int_{t^{k}}^{t^{k + 1}} L (x, \dot{x}) d t

那么，

\frac{\partial L_{d} (x^{k}, x^{k + 1})}{\partial μ} = \int_{t^{k}}^{t^{k + 1}} \sum_{i = 1}^{n} [\frac{\partial L (x, \dot{x})}{\partial x_{i}} \frac{\partial x_{i}}{\partial μ} + \frac{\partial L (x, \dot{x})}{\partial {\dot{x}}_{i}} \frac{\partial {\dot{x}}_{i}}{\partial μ}] d t

又因为 $\frac{\partial {\dot{x}}_{i}}{\partial μ} = \frac{d}{d t} \frac{\partial x_{i}}{\partial μ}$ ，所以

\frac{\partial L_{d} (x^{k}, x^{k + 1})}{\partial μ} = \int_{t^{k}}^{t^{k + 1}} \sum_{i = 1}^{n} [\frac{\partial L (x, \dot{x})}{\partial x_{i}} \frac{\partial x_{i}}{\partial μ} + \frac{\partial L (x, \dot{x})}{\partial {\dot{x}}_{i}} \frac{d}{d t} \frac{\partial x_{i}}{\partial μ}] d t

进一步地，利用分部积分公式可得

\frac{\partial L_{d} (x^{k}, x^{k + 1})}{\partial μ} = \int_{t^{k}}^{t^{k + 1}} \sum_{i = 1}^{n} (\frac{\partial L (x, \dot{x})}{\partial x_{i}} - \frac{d}{d t} \frac{\partial L (x, \dot{x})}{\partial {\dot{x}}_{i}}) \frac{\partial x_{i}}{\partial μ} d t + \sum_{i = 1}^{n} \frac{\partial L (x, \dot{x})}{\partial {\dot{x}}_{i}} {\frac{\partial x_{i}}{\partial μ}|}_{t^{k}}^{t^{k + 1}}

在上式中，记

S = \int_{t^{k}}^{t^{k + 1}} \sum_{i = 1}^{n} (\frac{\partial L (x, \dot{x})}{\partial x_{i}} - \frac{d}{d t} \frac{\partial L (x, \dot{x})}{\partial {\dot{x}}_{i}}) \frac{\partial x_{i}}{\partial μ} d t

（6）

则显然变分积分子构造过程中利用经典求积公式近似积分（2）的工作就转化为近似计算积分（6）.而如果进一步要求在离散时间结点 $\{t^{k} + c_{α} τ| c_{α} = α / m, α = 1, \dots, m - 1\}$ 处Euler-Lagrange方程（1）成立，那么近似计算积分（6）可得 $S = 0$ ，相应地

\frac{\partial L_{d} (x^{k}, x^{k + 1})}{\partial μ} = \sum_{i = 1}^{n} \frac{\partial L (x, \dot{x})}{\partial {\dot{x}}_{i}} {\frac{\partial x_{i}}{\partial μ}|}_{t^{k}}^{t^{k + 1}}

（7）

同理，当考虑时间区间 $[t^{k - 1}, t^{k}]$ 时，

\frac{\partial L_{d} (x^{k - 1}, x^{k})}{\partial μ} = \sum_{i = 1}^{n} \frac{\partial L (x, \dot{x})}{\partial {\dot{x}}_{i}} {\frac{\partial x_{i}}{\partial μ}|}_{t^{k - 1}}^{t^{k}}

假设在时间区间 $[t^{k}, t^{k + 1}]$ 内，

x (t) = l_{0} (t) x^{k} + \sum_{α = 1}^{m - 1} l_{α} (t) η^{α} + l_{m} (t) x^{k + 1}

（8）

其中，

\begin{array}{l} l_{α} (t) = \prod_{\begin{matrix} β = 0 \\ β \neq α \end{matrix}}^{m} \frac{t - (t^{k} + c_{β} τ)}{(t^{k} + c_{α} τ) - (t^{k} + c_{β} τ)} \\ α = 0,1, \dots, m \end{array}

而 $η^{α} = {(η_{1}^{α}, η_{2}^{α}, \dots η_{n}^{α})}^{T}$ 是拟合轨迹 $x (t)$ 在离散时间结点 $\{t^{k} + c_{α} τ| c_{α} = α / m, α = 1, \dots, m - 1\}$ 处的值，并且满足Euler-Lagrange方程，那么根据 $x (t)$ 的表达式可以计算 $\dot{x} (t)$ 、 $\ddot{x} (t)$ 以及

\begin{array}{l} {\frac{\partial x_{i} (t)}{\partial x_{i}^{k}}|}_{t = t^{k + 1}} = l_{0} (t^{k + 1}) = 0, \\ {\frac{\partial x_{i} (t)}{\partial x_{i}^{k}}|}_{t = t^{k}} = l_{0} (t^{k}) = 1 \end{array}

这样，由（7）式可得

\begin{array}{l} \frac{\partial L_{d} (x^{k}, x^{k + 1})}{\partial x_{i}^{k}} = {\frac{\partial L (x, \dot{x})}{\partial {\dot{x}}_{i}} \frac{\partial x_{i}}{\partial x_{i}^{k}}|}_{t^{k}}^{t^{k + 1}} \\ = - \frac{\partial L (x (t^{k}), \dot{x} (t^{k}))}{\partial {\dot{x}}_{i}} \end{array}

（9）

同理，当考虑时间区间 $[t^{k - 1}, t^{k}]$ 时，

\begin{array}{l} \frac{\partial L_{d} (x^{k - 1}, x^{k})}{\partial x_{i}^{k}} = {\frac{\partial L (x, \dot{x})}{\partial {\dot{x}}_{i}} \frac{\partial x_{i}}{\partial x_{i}^{k}}|}_{t^{k - 1}}^{t^{k}} \\ = \frac{\partial L (x (t^{k}), \dot{x} (t^{k}))}{\partial {\dot{x}}_{i}} \end{array}

（10）

利用（9）式和（10）式，联立方程组

\frac{\partial L_{d} (x^{k}, x^{k + 1})}{\partial x_{i}^{k}} + \frac{\partial L_{d} (x^{k - 1}, x^{k})}{\partial x_{i}^{k}} = - \frac{\partial L (x (t^{k}), \dot{x} (t^{k}))}{\partial {\dot{x}}_{i}} + \frac{\partial L (x (t^{k}), \dot{x} (t^{k}))}{\partial {\dot{x}}_{i}} = 0

（11）

{{\frac{\partial L (x, \dot{x})}{\partial x_{i}}|}_{t = t^{k} + c_{α} τ} - \frac{d}{d t} \frac{\partial L (x, \dot{x})}{\partial {\dot{x}}_{i}}|}_{t = t^{k} + c_{α} τ} = 0

（12）

其中， $i = 1,2, \dots, n$ ， $α = 1,2, \dots, m - 1$ ，则首先由方程（12）可求得 $η^{α} = η^{α} (x^{k}, x^{k + 1})$ ，平行地，当时间区间为 $[t^{k - 1}, t^{k}]$ 时， $η^{α} = η^{α} (x^{k - 1}, x^{k})$ ；然后将求得的 $η^{α}$ 代入 $x (t^{k})$ 和 $\dot{x} (t^{k})$ ，并进一步带入（9）式和（10）式，那么方程（11）就转化关于 $x^{k - 1}$ 、 $x^{k}$ 和 $x^{k + 1}$ 的代数方程，即变分积分子 $Φ : (x^{k - 1}, x^{k}) \to (x^{k}, x^{k + 1})$ .

注1：单个步长区间 $[t^{k}, t^{k + 1}]$ 内，在对局部轨迹 $x (t)$ 进行插值拟合时，除了Lagrange插值多项式（8）之外，还可以采用其他的插值方式，例如二项插值

x (t) = b_{0, m} (ε) x^{k} + \sum_{α = 1}^{m - 1} b_{α, m} (ε) η^{α} + b_{m, m} (ε) x^{k + 1}

（13）

其中， $b_{α, m} (ε) = C_{m}^{α} ε^{α} {(1 - ε)}^{m - α}$ ， $ε = \frac{t - t^{k}}{τ}$ .由（13）式可得

\begin{array}{l} {\frac{\partial x_{i} (t)}{\partial x_{i}^{k}}|}_{t = t^{k + 1}} = {b_{0, m} (ε)|}_{ε = 1} = 0, \\ {\frac{\partial x_{i} (t)}{\partial x_{i}^{k}}|}_{t = t^{k}} = {b_{0, m} (ε)|}_{ε = 0} = 1 \end{array}

代入偏导数 $\partial_{1} L_{d}$ 和 $\partial_{2} L_{d}$ ，同样可以得到

\begin{array}{l} \frac{\partial L_{d} (x^{k}, x^{k + 1})}{\partial x_{i}^{k}} = - \frac{\partial L (x (t^{k}), \dot{x} (t^{k}))}{\partial {\dot{x}}_{i}}, \\ \frac{\partial L_{d} (x^{k - 1}, x^{k})}{\partial x_{i}^{k}} = \frac{\partial L (x (t^{k}), \dot{x} (t^{k}))}{\partial {\dot{x}}_{i}} \end{array}

进一步联立连续和离散的Euler-Lagrange方程（11-12）也可得到变分积分子.前文曾提及，对于完整约束系统（4），其变分积分子（5）等价于以 ${\bar{L}}_{d} ({\bar{x}}^{k}, {\bar{x}}^{k + 1})$ 为离散Lagrange函数的离散Euler-Lagrange方程

\partial_{1} {\bar{L}}_{d} ({\bar{x}}^{k}, {\bar{x}}^{k + 1}) + \partial_{2} {\bar{L}}_{d} ({\bar{x}}^{k - 1}, {\bar{x}}^{k}) = 0

（14）

而（14）式同样只包含了偏导数 $\partial_{1} {\bar{L}}_{d}$ 和 $\partial_{2} {\bar{L}}_{d}$ ，因此可以仿照自由系统利用局部路径拟合方法构造完整约束系统的变分积分子.此时，

{\bar{L}}_{d} ({\bar{x}}^{k}, {\bar{x}}^{k + 1}) = \int_{t^{k}}^{t^{k + 1}} L (x, \dot{x}) d t - \int_{t^{k}}^{t^{k + 1}} \sum_{j = 1}^{p} λ_{j} (t) f_{j} (x) d t

对应地

\begin{array}{l} \frac{\partial {\bar{L}}_{d} ({\bar{x}}^{k}, {\bar{x}}^{k + 1})}{\partial μ} = \int_{t^{k}}^{t^{k + 1}} \sum_{i = 1}^{n} (\frac{\partial L}{\partial x_{i}} - \frac{d}{d t} \frac{\partial L}{\partial {\dot{x}}_{i}} - \sum_{j = 1}^{p} λ_{j} \frac{\partial f_{j}}{\partial x_{i}}) \frac{\partial x_{i}}{\partial μ} d t \\ - \int_{t^{k}}^{t^{k + 1}} \sum_{j = 1}^{p} \frac{\partial λ_{j}}{\partial μ} f_{j} d t + \sum_{i = 1}^{n} \frac{\partial L}{\partial x_{i}} {\frac{\partial x_{i}}{\partial μ}|}_{t^{k}}^{t^{k + 1}} \end{array}

记

\begin{array}{l} S_{1} = \int_{t^{k}}^{t^{k + 1}} \sum_{i = 1}^{n} (\frac{\partial L}{\partial x_{i}} - \frac{d}{d t} \frac{\partial L}{\partial {\dot{x}}_{i}} - \sum_{j = 1}^{p} λ_{j} \frac{\partial f_{j}}{\partial x_{i}}) \frac{\partial x_{i}}{\partial μ} d t, \\ S_{2} = \int_{t^{k}}^{t^{k + 1}} \sum_{j = 1}^{p} \frac{\partial λ_{j}}{\partial μ} f_{j} d t \end{array}

如果假设在离散结点 ${t^{k} + c_{α} τ| c_{α} = α / m, α = 1, \dots, m - 1}$ 处，约束系统的运动方程（4）成立，那么基于这一组离散结点数值计算积分 $S_{1}$ 和 $S_{2}$ 则有 $S_{1} = S_{2} = 0$ ，于是

\frac{\partial {\bar{L}}_{d} ({\bar{x}}^{k}, {\bar{x}}^{k + 1})}{\partial μ} = \sum_{i = 1}^{n} \frac{\partial L}{\partial x_{i}} {\frac{\partial x_{i}}{\partial μ}|}_{t^{k}}^{t^{k + 1}}

平行地

\frac{\partial {\bar{L}}_{d} ({\bar{x}}^{k - 1}, {\bar{x}}^{k})}{\partial μ} = \sum_{i = 1}^{n} \frac{\partial L}{\partial x_{i}} {\frac{\partial x_{i}}{\partial μ}|}_{t^{k - 1}}^{t^{k}}

在时间区间 $[t^{k}, t^{k + 1}]$ 内，如果局部轨迹 $x (t)$ 和 $λ (t)$ 都按照Lagrange插值多项式进行拟合，则最终可得

\begin{array}{l} \frac{\partial {\bar{L}}_{d} ({\bar{x}}^{k}, {\bar{x}}^{k + 1})}{\partial x_{i}^{k}} = - \frac{\partial L (x (t^{k}), \dot{x} (t^{k}))}{\partial {\dot{x}}_{i}}, \\ \frac{\partial {\bar{L}}_{d} ({\bar{x}}^{k - 1}, {\bar{x}}^{k})}{\partial x_{i}^{k}} = \frac{\partial L (x (t^{k}), \dot{x} (t^{k}))}{\partial {\dot{x}}_{i}}, \\ \frac{\partial {\bar{L}}_{d} ({\bar{x}}^{k}, {\bar{x}}^{k + 1})}{\partial λ_{i}^{k}} = 0, \\ \frac{\partial {\bar{L}}_{d} ({\bar{x}}^{k - 1}, {\bar{x}}^{k})}{\partial λ_{i}^{k}} = 0, \end{array}

而方程（14）则随之变为

\frac{\partial {\bar{L}}_{d} ({\bar{x}}^{k}, {\bar{x}}^{k + 1})}{\partial x_{i}^{k}} + \frac{\partial {\bar{L}}_{d} ({\bar{x}}^{k - 1}, {\bar{x}}^{k})}{\partial x_{i}^{k}} = - \frac{\partial L (x (t^{k}), \dot{x} (t^{k}))}{\partial {\dot{x}}_{i}} + \frac{\partial L (x (t^{k}), \dot{x} (t^{k}))}{\partial {\dot{x}}_{i}} = 0

再考虑到所作假设—方程（4）在离散结点 ${t^{k} + c_{α} τ| c_{α} = α / m, α = 1, \dots, m - 1}$ 处成立，结合在一起就得到了完整约束系统的变分积分子.

3 数值算例

3.1　算例1—简谐振子

考虑简谐振子方程

\ddot{x} (t) + x (t) = 0

（15）

取Lagrange函数 $L = \frac{1}{2} {\dot{x}}^{2} - \frac{1}{2} x^{2}$ ，则微分方程等价于Euler-Lagrange方程 $\frac{\partial L}{\partial x} - \frac{d}{d t} \frac{\partial L}{\partial \dot{x}} = 0$ .

在时间区间 $[t^{k}, t^{k + 1}]$ 内，令

x (t) = \frac{t - t^{k + \frac{1}{2}}}{t^{k} - t^{k + \frac{1}{2}}} \cdot \frac{t - t^{k + 1}}{t^{k} - t^{k + 1}} x^{k} + \frac{t - t^{k}}{t^{k + \frac{1}{2}} - t^{k}} \cdot \frac{t - t^{k + 1}}{t^{k + \frac{1}{2}} - t^{k + 1}} η^{1} + \frac{t - t^{k}}{t^{k + 1} - t^{k}} \cdot \frac{t - t^{k + \frac{1}{2}}}{t^{k + 1} - t^{k + \frac{1}{2}}} x^{k + 1}

则求导可得

\begin{array}{l} \dot{x} (t) = \frac{4 t - 2 (t^{k + \frac{1}{2}} + t^{k + 1})}{τ^{2}} x^{k} - \\ \frac{8 t - 4 (t^{k} + t^{k + 1})}{τ^{2}} η^{1} + \frac{4 t - 2 (t^{k} + t^{k + \frac{1}{2}})}{τ^{2}} x^{k + 1} \end{array}

\ddot{x} (t) = \frac{4}{τ^{2}} x^{k} - \frac{8}{τ^{2}} η^{1} + \frac{4}{τ^{2}} x^{k + 1}

在 $t = t^{k} + τ / 2$ 处，由Euler-Lagrange方程，即系统运动方程（15）成立可得

η^{1} = \frac{4}{8 - τ^{2}} (x^{k} + x^{k + 1})

将 $η^{1}$ 带入 $\dot{x} (t)$ 相应可得

\partial_{1} L_{d} (x^{k}, x^{k + 1}) = - {\frac{\partial L}{\partial \dot{x}}|}_{t = t^{k}} = - \dot{x} (t^{k}) = - \frac{3 τ^{2} - 8}{τ (8 - τ^{2})} x^{k} - \frac{τ^{2} + 8}{τ (8 - τ^{2})} x^{k + 1}

\partial_{2} L_{d} (x^{k - 1}, x^{k}) = {\frac{\partial L}{\partial \dot{x}}|}_{t = t^{k}} = \dot{x} (t^{k}) = - \frac{τ^{2} + 8}{τ (8 - τ^{2})} x^{k - 1} - \frac{3 τ^{2} - 8}{τ (8 - τ^{2})} x^{k}

那么，由离散Euler-Lagrange方程（11）可得

x^{k + 1} = \frac{16 - 6 τ^{2}}{8 + τ^{2}} x^{k} - x^{k - 1}

（16）

迭代算法（16）就是利用新的构造方法针对简谐振子方程（15）设计的变分积分子.假定系统（15）的初值为 $x (0) = \dot{x} (0) = 1$ ，以 $τ = 0.1$ 为时间步长，利用差分格式（16）对系统进行数值求解，则对应的数值解如图1所示.显然，该差分格式非常精确地数值模拟了简谐振动，在图示的尺度下，方程的数值解和精确解所对应的曲线完全重合，无任何差别.

图1 方程(15)的精确解与通过差分格式(16)算得的数值解对照

Fig.1 Comparison between the exact solution and the numerical solution computed by algorithm (16) for equation (15)

除了有效性得以验证之外，由图2可以非常直观地看到，差分格式（16）也非常真实地模拟了系统的能量，所计算出的数值能量曲线呈现出长时间的守恒行为，且能量误差始终保持在有界的范围内.保持保守系统的守恒量，这是变分积分子特有的计算优越性，因此，图2中的数值结果也验证了通过局部路径拟合方法构造的变分积分子继承了经典变分积分子优越的计算性能.实际上不仅如此，在同阶精度的情况下，借助新方法得到的变分积分子具有更小的误差.

图2 由差分格式(16)算得的系统(15)的数值能量曲线

Fig.2 The numerical energy curve of system (15) computed by algorithm (16)

为了验证这一事实，采用经典求积公式—中点格式—逼近（2）式中的积分可得与系统（15）对应的离散Lagrange函数

L_{d} (x^{k}, x^{k + 1}) = \frac{τ}{2} {(\frac{x^{k + 1} - x^{k}}{τ})}^{2} - \frac{τ}{2} {(\frac{x^{k + 1} + x^{k}}{2})}^{2}

根据离散Euler-Lagrange方程（3）可得经典变分积分子

x^{k + 1} = \frac{8 - 2 τ^{2}}{4 + τ^{2}} x^{k} - x^{k - 1}

（17）

差分格式（16）和差分格式（17）尽管都是具有二阶精度的数值算法，然而由表1中数据可以直观地对比出差分格式（16）具有更小的局部误差和整体误差.

表1 不同时刻差分格式(16)和差分格式(17)的误差对比

Table 1 Local errors of numerical solutions and energies at different times computed with algorithms (16) and (17) respectively

Time	Algorithm (16) Error of solution/m	Algorithm (16) Error of energy/J	Algorithm (17) Error of solution/m	Algorithm (17) Error of energy/J
10s	0.0007	0.0001	0.0023	0.0186
20s	0.0019	0.0000	0.0086	0.0339
30s	0.0069	0.0006	0.0286	0.0463
40s	0.0116	0.0011	0.0472	0.0116
50s	0.0128	0.0008	0.0518	0.0128

当然，如果单个步长区间 $[t^{k}, t^{k + 1}]$ 内的离散结点取得更加密集，那么可以得到更加精确的差分格式.例如，在时间区间 $[t^{k}, t^{k + 1}]$ 内，如果设

\begin{array}{l} x (t) = \frac{t - t^{k + \frac{1}{3}}}{t^{k} - t^{k + \frac{1}{3}}} \cdot \frac{t - t^{k + \frac{2}{3}}}{t^{k} - t^{k + \frac{2}{3}}} \cdot \frac{t - t^{k + 1}}{t^{k} - t^{k + 1}} x^{k} \\ + \frac{t - t^{k}}{t^{k + \frac{1}{3}} - t^{k}} \cdot \frac{t - t^{k + \frac{2}{3}}}{t^{k + \frac{1}{3}} - t^{k + \frac{2}{3}}} \cdot \frac{t - t^{k + 1}}{t^{k + \frac{1}{3}} - t^{k + 1}} η^{1} \\ + \frac{t - t^{k}}{t^{k + \frac{2}{3}} - t^{k}} \cdot \frac{t - t^{k + \frac{1}{3}}}{t^{k + \frac{2}{3}} - t^{k + \frac{1}{3}}} \cdot \frac{t - t^{k + 1}}{t^{k + \frac{2}{3}} - t^{k + 1}} η^{2} \\ + \frac{t - t^{k}}{t^{k + 1} - t^{k}} \cdot \frac{t - t^{k + \frac{1}{3}}}{t^{k + 1} - t^{k + \frac{1}{3}}} \cdot \frac{t - t^{k + \frac{2}{3}}}{t^{k + 1} - t^{k + \frac{2}{3}}} x^{k + 1} \end{array}

并让系统运动方程在 $t = t^{k} + \frac{1}{3} τ$ 和 $t = t^{k} + \frac{2}{3} τ$ 处成立，则可解得

η^{1} = \frac{81}{(τ^{2} - 9) (τ^{2} - 27)} x^{k + 1} - \frac{9 (τ^{2} - 18)}{(τ^{2} - 9) (τ^{2} - 27)} x^{k}

η^{2} = \frac{81}{(τ^{2} - 9) (τ^{2} - 27)} x^{k} - \frac{9 (τ^{2} - 18)}{(τ^{2} - 9) (τ^{2} - 27)} x^{k + 1}

对应地

\partial_{1} L_{d} (x^{k}, x^{k + 1}) = \frac{11 τ^{4} - 234 τ^{2} + 486}{2 τ (τ^{2} - 9) (τ^{2} - 27)} x^{k} - \frac{2 τ^{4} + 9 τ^{2} + 486}{2 τ (τ^{2} - 9) (τ^{2} - 27)} x^{k + 1}

\partial_{2} L_{d} (x^{k - 1}, x^{k}) = - \frac{2 τ^{4} + 9 τ^{2} + 486}{2 τ (τ^{2} - 9) (τ^{2} - 27)} x^{k - 1} + \frac{11 τ^{4} - 234 τ^{2} + 486}{2 τ (τ^{2} - 9) (τ^{2} - 27)} x^{k}

代入离散Euler-Lagrange方程（11）可得差分格式

x^{k + 1} = \frac{22 τ^{4} - 468 τ^{2} + 972}{2 τ^{4} + 9 τ^{2} + 486} x^{k} - x^{k - 1}

（18）

又例如，如果在 $[t^{k}, t^{k + 1}]$ 内再多取一个结点，即令

x (t) = l_{0} (t) x^{k} + l_{1} (t) η^{1} + l_{2} (t) η^{2} + l_{3} (t) η^{3} + l_{4} (t) x^{k + 1}

则最终可得变分积分子

x^{k + 1} = \frac{- 150 τ^{6} + 9616 τ^{4} - 142848 τ^{2} + 294912}{9 τ^{6} - 184 τ^{4} + 2304 τ^{2} + 147456} x^{k} - x^{k - 1}

（19）

如表2中数据所示，差分格式（18）和差分格式（19）较之差分格式（16）具有更小的误差.内部结点选取得越密集，得到的数值算法越精确，这一规律与经典构造过程中愈加复杂的积分逼近格式诱导愈加精确的变分积分子相一致.

表2 不同时刻差分格式(16-19)的数值解误差对比

Table 2 Local errors of numerical solutions at different times computed with algorithms (16-19) respectively

Time	Algorithm (16)	Algorithm (17)	Algorithm (18)	Algorithm (19)
10s	7.1987 $/ 10^{4}$	2.3042 $/ 10^{3}$	1.3519 $/ 10^{4}$	2.0468 $/ 10^{4}$
20s	1.9163 $/ 10^{3}$	8.5537 $/ 10^{3}$	1.8646 $/ 10^{3}$	3.4417 $/ 10^{4}$
30s	6.9258 $/ 10^{3}$	2.8618 $/ 10^{2}$	5.2588 $/ 10^{3}$	3.7354 $/ 10^{4}$
40s	1.1606 $/ 10^{2}$	4.7230 $/ 10^{2}$	8.2132 $/ 10^{3}$	2.8310 $/ 10^{4}$
50s	1.2754 $/ 10^{2}$	5.1768 $/ 10^{2}$	8.6299 $/ 10^{3}$	1.0158 $/ 10^{4}$

注：对于简谐振子方程(15)，如果基于二项插值(13)对局部轨迹 $x (t)$ 进行插值拟合，那么最终所得到的差分格式与基于Lagrange插值得到的差分格式完全一致，即同样可以对应得到差分格式(16)、(18)和(19)，只不过诱导过程中涉及到的中间变量，例如 $\dot{x} (t)$ 、 $\ddot{x} (t)$ 、 $η^{1}$ 等的表达式并不一致.

3.2　算例2—单摆

单摆运动方程

\ddot{θ} + s i n θ = 0

（20）

也可以等价地表述为Euler-Lagrange方程，对应的Lagrange函数为

L (θ, \dot{θ}) = \frac{1}{2} {(\dot{θ})}^{2} - (1 - c o s θ)

在单个步长区间内如果只选取一个结点，那么对于系统（20）可得如下方程组

\begin{array}{l} \frac{4}{τ^{2}} θ^{k} - \frac{8}{τ^{2}} η^{1} + \frac{4}{τ^{2}} θ^{k + 1} + s i n η^{1} = 0 \\ θ^{k - 1} + 6 θ^{k} + θ^{k + 1} - 4 (η^{1} + {\bar{η}}^{1}) = 0 \\ \frac{4}{τ^{2}} θ^{k - 1} - \frac{8}{τ^{2}} {\bar{η}}^{1} + \frac{4}{τ^{2}} θ^{k} + s i n {\bar{η}}^{1} = 0 \end{array}

（21）

或

\begin{array}{l} \frac{2}{τ^{2}} θ^{k} - \frac{4}{τ^{2}} η^{1} + \frac{2}{τ^{2}} θ^{k + 1} = - s i n (\frac{1}{4} θ^{k} + \frac{1}{2} η^{1} + \frac{1}{4} θ^{k + 1}) \\ η^{1} = 2 θ^{k} - {\bar{η}}^{1} \\ \frac{2}{τ^{2}} θ^{k - 1} - \frac{4}{τ^{2}} {\bar{η}}^{1} + \frac{2}{τ^{2}} θ^{k} = - s i n (\frac{1}{4} θ^{k - 1} + \frac{1}{2} {\bar{η}}^{1} + \frac{1}{4} θ^{k}) \end{array}

（22）

其中方程组（21）是基于Lagrange插值（8）诱导得来，而方程组（22）是基于二项插值（13）得来.隐式地求解方程组（21）或（22）均可得变分积分子 $(θ^{k}, θ^{k + 1}) = Φ (θ^{k - 1}, θ^{k})$ .利用两种差分格式算得系统（20）在相空间中的数值轨迹和数值能量如图3所示.显然，两种差分格式不仅有效地模拟了系统的运动，而且非常真实地反映了系统能量守恒.

图3 由差分格式(21)和(22)分别算得的系统(20)在相空间中的数值轨迹(a)和数值能量曲线(b)

Fig.3 Numerical orbits in the phase space (a) and numerical energy curves (b) of system (20) computed by algorithms (21) and (22) respectively

3.3　算例3—匀速圆周运动

单位质量的质点在二维平面内做匀速圆周运动，其运动方程可以表述为

\{\begin{matrix} {\ddot{x}}_{1} = - 2 λ x_{1}, \\ {\ddot{x}}_{2} = - 2 λ x_{2}, \end{matrix}

其中位形坐标 $(x_{1}, x_{2})$ 满足约束条件 $x_{1}^{2} + x_{2}^{2} = 1$ .基于Lagrange插值多项式（8），可得该完整约束系统的变分积分子

\begin{array}{l} x_{1}^{k - 1} + 6 x_{1}^{k} + x_{1}^{k + 1} - \frac{4 (x_{1}^{k} + x_{1}^{k + 1})}{\sqrt[]{{(x_{1}^{k} + x_{1}^{k + 1})}^{2} + {(x_{2}^{k} + x_{2}^{k + 1})}^{2}}} - \\ \frac{4 (x_{1}^{k - 1} + x_{1}^{k})}{\sqrt[]{{(x_{1}^{k - 1} + x_{1}^{k})}^{2} + {(x_{2}^{k - 1} + x_{2}^{k})}^{2}}} = 0, \\ x_{2}^{k - 1} + 6 x_{2}^{k} + x_{2}^{k + 1} - \frac{4 (x_{2}^{k} + x_{2}^{k + 1})}{\sqrt[]{{(x_{1}^{k} + x_{1}^{k + 1})}^{2} + {(x_{2}^{k} + x_{2}^{k + 1})}^{2}}} - \\ \frac{4 (x_{2}^{k - 1} + x_{2}^{k})}{\sqrt[]{{(x_{1}^{k - 1} + x_{1}^{k})}^{2} + {(x_{2}^{k - 1} + x_{2}^{k})}^{2}}} = 0, \end{array}

（23）

式（23）构成了一个具有二阶精度的隐式差分格式.

4 结论

本文结合局部路径拟合方法提出了一种新的力学系统变分积分子的设计方式.数值结果表明，这种基于局部路径拟合的构造方法不仅让最终得到的数值算法仍然继承了变分积分子特有的优越计算性能，而且显著地提高了算法的精度.另外，这种构造方法不仅适用于自由无约束系统，而且适用于完整约束系统.

参考文献

梅凤翔. 分析力学. 北京：北京理工大学出版社，2013 （Mei F X. Analytical mechanics. Beijing： Beijing Institute of Technology Press，2013 （in Chinese）） [百度学术]

冯康，秦孟兆. 哈密尔顿系统的辛几何算法. 杭州：浙江科技出版社，2003 （Feng K，Qin M Z. Symplectic geometric algorithms for hamiltonian systems. Hangzhou： Zhejiang Science & Technology Press，2003 （in Chinese）） [百度学术]

Hairer E，Lubich C，Wanner G. Geometric numerical integration： structure preserving algorithms for ordinary differential equations. Berlin：Springer，2006 [百度学术]

Wendlandt J，Marsden J E. Mechanical integrators derived from a discrete variational principle. Physica D， 1997， 106： 223~246 [百度学术]

Marsden J E，West M. Discrete mechanics and variational integrators. Acta Numerica， 2001， 10： 357~514 [百度学术]

Kane C，Marsden J E，Ortiz M，et al. Variational integrators and the Newmark algorithm for conservative and dissipative mechanical systems. International Journal for Numerical Methods in Engineering， 2000， 49：1295~1325 [百度学术]

Fetecau R C，Marsden J E，Ortiz M，et al. Nonsmooth Lagrangian mechanics and variational collision integrators. SIAM Journal on Applied Dynamical Systems， 2003， 2： 381~416 [百度学术]

Wenger T，Ober-Blobaum S，Leyendecker S. Construction and analysis of higher order variational integrators for dynamical systems with holonomic constraints. Advances in Computational Mathematics， 2017， 43（5）： 1163~1195 [百度学术]

Ferraro S，Iglesias D，de Diego D M. Momentum and energy preserving integrators for nonholonomic dynamics. Nonlinearity， 2008， 21（8）： 1911~1928 [百度学术]

Ferraro S，Jimenez F，de Diego D M. New developments on the geometric nonholonomic integrator. Nonlinearity， 2015， 28（4）： 871~900 [百度学术]

刘世兴，刘畅，郭永新. Birkhoff意义下Hénon-Heiles方程的离散变分计算. 物理学报， 2011， 60（6）： 064501（Liu S X，Liu C，Guo Y X. Discrete variational calculation of Hénon-Heiles equation in the Birkhoffian sense. Acta Physica Sinica， 2011， 60（6）： 064501 （in Chinese）） [百度学术]

Liu S X，Liu C，Guo Y X. Geometric formulations and variational integrators of discrete autonomous Birkhoff systems. Chinese Physics B， 2011， 20（3）： 034501 [百度学术]

Kong X L，Wu H B，Mei F X. Structure-preserving algorithms for Birkhoffian systems. Journal of Geometry and Physics， 2012， 62（5）： 1157~1166 [百度学术]

Gawlik E S，Mullen P，Pavlov D，et al. Geometric， variational discretization of continuum theories. Physica D， 2011， 240（21）： 1724~1760 [百度学术]

Kraus M，Tassi E，Grasso D. Variational integrators for reduced magneto hydrodynamics. Journal of Computational Physics， 2016， 321： 435~458 [百度学术]

Bou-Rabee N，Owhadi H. Stochastic variational integrators. IMA Journal of Numerical Analysis， 2009， 29（2）： 421~443 [百度学术]

Bou-Rabee N，Owhadi H. Long-run accuracy of variational integrators in the stochastic context. SIAM Journal on Numerical Analysis， 2010， 48（1）： 278~297 [百度学术]

Kosmas O，Vlachos D S. Local path fitting： A new approach to variational integrators. Journal of Computational and Applied Mathematics， 2012， 236（10）： 2632~2642 [百度学术]

Kosmas O，Leyendecker S. Analysis of higher order phase fitted variational integrators. Advances in Computational Mathematics， 2016， 42（3）： 605~619 [百度学术]

基于Lagrange插值多项式拟合的力学系统的变分积分子

摘要

关键词

引言

1 变分积分子

2 基于Lagrange插值多项式拟合的力学系统的变分积分子