中心差分法：简要总结

中心差分方法是一种经典的显式时间积分算法，广泛应用于结构动力学、有限元分析及分子动力学模拟中。

1. 核心原理与公式推导

中心差分法基于泰勒级数展开，利用当前时刻 t 前后两个时刻（ t-\Delta t 和 t+\Delta t ）的位移来近似计算 t 时刻的速度和加速度。

速度与加速度近似公式：通过前差分和后差分公式相减与相加，忽略高阶小量 O(\Delta t^2) ，得到：
- 速度： \dot{u}_t \approx \frac{1}{2\Delta t}(u_{t+\Delta t} - u_{t-\Delta t}) ，
- 加速度： \ddot{u}_t \approx \frac{1}{\Delta t^2}(u_{t+\Delta t} - 2u_t + u_{t-\Delta t}) ，
运动方程离散化：将上述近似式代入 t 时刻的运动方程 M\ddot{u}_t + C\dot{u}_t + Ku_t = Q_t ，整理后得到求解 t+\Delta t 时刻位移的递推公式：\hat{M} u_{t+\Delta t} = \hat{Q}_t其中：
- 有效质量矩阵： \hat{M} = \frac{1}{\Delta t^2}M + \frac{1}{2\Delta t}C ，
- 有效载荷向量： \hat{Q}_t = Q_t - (K - \frac{2}{\Delta t^2}M)u_t - (\frac{1}{\Delta t^2}M - \frac{1}{2\Delta t}C)u_{t-\Delta t} ，

2. 算法特性：显式积分 (Explicit Integration)

中心差分法最显著的特征是其显式性质，这带来了巨大的计算优势：

无需刚度矩阵求逆/分解：在递推公式中，刚度矩阵 K 不出现在待求解矩阵（有效质量矩阵 \hat{M} ）的左侧。因此，不需要对 K 进行三角分解或求逆。
无需组装总体刚度矩阵：由于 Ku_t 项可以通过单元级别的弹性力累加得到（ F_E = \sum K_e u_e ），在有限元分析中可以完全不组装全局刚度矩阵，极大地节省了内存存储量。
对角化加速：如果质量矩阵 M 是对角阵（集中质量矩阵）且阻尼 C 可忽略或也是对角阵，则有效质量矩阵 \hat{M} 也是对角阵。此时求解 u_{t+\Delta t} 不需要矩阵运算，直接对每个自由度进行标量除法即可：u_{t+\Delta t}^{(i)} = \frac{\hat{Q}_t^{(i)}}{\hat{M}_{ii}}

这使得该方法非常适合大规模并行计算。

3. 起步程序 (Starting Procedure)

由于公式需要 u_{t-\Delta t} 来计算 u_{t+\Delta t} ，而在初始时刻 t=0 时只有 u_0 和 \dot{u}_0 ，因此需要特殊处理来启动计算：

虚拟位移计算：利用后差分公式反推 t=-\Delta t 时刻的位移 u_{-\Delta t} ：
u_{-\Delta t} = u_0 - \Delta t \dot{u}_0 + \frac{\Delta t^2}{2} \ddot{u}_0
初始加速度求解：其中的 \ddot{u}_0 需通过 t=0 时刻的运动方程求得： M\ddot{u}_0 = Q_0 - C\dot{u}_0 - Ku_0 。

4. 常用变异格式

文档介绍了三种常见的变体，适用于不同场景：

Verlet 格式：
- 直接利用 u_{t+\Delta t} = 2u_t - u_{t-\Delta t} + \ddot{u}_t \Delta t^2 。
- 特点：位移截断误差为 O(\Delta t^4) ，精度较高；但不直接输出速度，速度需额外计算且误差为 O(\Delta t^2) 。
蛙跳格式 (Leapfrog)：
- 位移和速度在时间轴上交错求解（速度在 t+\Delta t/2 ，位移在 t+\Delta t ）。
- 特点：同样不直接给出同一时刻的速度和位移。
速度 Verlet 格式 (Velocity Verlet / Leapfrog Verlet)：
- 结合了蛙跳的思想，但通过两步更新速度，使得能在同一时刻获得位移、速度和加速度。
- 应用：是显式有限元程序（如 LS-DYNA）和分子动力学程序中最常用的格式。
- 步骤：
  1. 更新半步速度；
  2. 更新全步位移；
  3. 计算新时刻加速度；
  4. 更新另半步速度得到全步速度。

5. 稳定性分析 (Stability Analysis)

中心差分法是条件稳定的，其稳定性取决于时间步长 \Delta t 与系统固有特性的关系。

无阻尼系统临界步长：为了保证递推矩阵特征值的模 |\lambda| \le 1 ，时间步长必须满足：\Delta t \le \Delta t_{cr} = \frac{T_n}{\pi}其中 T_n 是系统的最小固有周期（对应最高频率）。
有阻尼系统：阻尼的存在会略微改变临界步长公式，但基本原则不变： \Delta t 必须小于由系统最高频率决定的临界值。
数值阻尼特性：在稳定范围内（ \Delta t \le \Delta t_{cr} ），中心差分法的谱半径 \rho(D)=1 ，意味着它没有数值阻尼（人工阻尼）。这对于保留低频精度是有利的，但也意味着无法自动滤除由离散化引起的高频虚假振荡。

6. 实际应用中的限制与注意事项

时间步长极小：由于 \Delta t 受限于最小单元的尺寸（ T_n 很小），对于精细网格的有限元模型， \Delta t 可能非常小，导致需要成千上万个时间步才能完成模拟，计算成本高。
质量矩阵要求：要求质量矩阵的所有对角元素 m_{ii} > 0 。如果某些节点质量为 0 或接近 0，会导致 \Delta t_{cr} \to 0 ，算法失效。
网格均匀性：在有限元建模时，应避免个别单元尺寸远小于其他单元，否则整个模型的 \Delta t 将被这些微小单元“拖累”，大幅降低计算效率。

7. 算法流程总结

根据文档提供的 Fortran 代码逻辑，标准实施步骤为：

初始化：形成 M, C, K （或仅 M, C ），给定初值，计算 u_{-\Delta t} ，计算积分常数 c_0, c_1, c_2, c_3 ，形成并分解有效质量矩阵 \hat{M} 。
循环迭代 (每步 \Delta t )：
- 计算 t 时刻的有效载荷 \hat{Q}_t 。
- 求解 \hat{M} u_{t+\Delta t} = \hat{Q}_t 得到新位移。
- （可选）回代计算 t 时刻的速度和加速度。
- 更新变量 ( u_{t-\Delta t} \leftarrow u_t , u_t \leftarrow u_{t+\Delta t} )，时间推进。

中心差分法以其显式、无需迭代、内存占用低的特点，成为解决高频、瞬态、非线性及接触问题（如冲击、爆炸）的首选方法。然而，其条件稳定性强制要求极小的时间步长，使其不适用于低频主导或长时程的动力学问题（此类问题通常选用 Newmark- \beta 法等隐式算法）。在实际工程中，往往需要在网格精度（影响最高频率）和计算成本之间寻找平衡。

中心差分法：算法流程

输入:

\mathbf{M}, \mathbf{C}, \mathbf{K} : 质量、阻尼、刚度矩阵
\mathbf{u}_0, \dot{\mathbf{u}}_0 : 初始位移和速度
\mathbf{Q}(t) : 外部载荷函数
\Delta t : 时间步长
N_{steps} : 总时间步数

输出:

\mathbf{u}_t, \dot{\mathbf{u}}_t, \ddot{\mathbf{u}}_t : 各时刻的位移、速度、加速度

1. 初始化阶段

计算积分常数: c_0 = 1 / \Delta t^2 c_1 = 1 / (2\Delta t) c_2 = 2 \times c_0 c_3 = \Delta t^2 / 2
计算初始加速度 ( \ddot{\mathbf{u}}_0 ):
- \ddot{\mathbf{u}}_0 = \mathbf{M}^{-1} (\mathbf{Q}_0 - \mathbf{C}\dot{\mathbf{u}}_0 - \mathbf{K}\mathbf{u}_0) ,
计算虚拟起始位移 ( \mathbf{u}_{-\Delta t} ):
- \mathbf{u}_{-\Delta t} = \mathbf{u}_0 - \Delta t \cdot \dot{\mathbf{u}}_0 + c_3 \cdot \ddot{\mathbf{u}}_0 ,
设置迭代变量:
- \mathbf{u}_{prev} \leftarrow \mathbf{u}_{-\Delta t} (代表 t-\Delta t )
- \mathbf{u}_{curr} \leftarrow \mathbf{u}_0 (代表 t )
构建并处理有效质量矩阵 ( \hat{\mathbf{M}} ):
- \hat{\mathbf{M}} = c_0 \mathbf{M} + c_1 \mathbf{C} ,
- 如果 \hat{\mathbf{M}} 不是对角矩阵: - 对 \hat{\mathbf{M}} 进行三角分解 (例如 LDL^T 分解)
- 否则 (若为对角阵): - 直接存储其对角元素的倒数用于后续除法运算

2. 时间步循环迭代

对于 k = 1 到 N_{steps} 执行:

更新当前时间:
- t = (k-1) \times \Delta t ,
获取当前载荷:
- \mathbf{Q}_t = \mathbf{Q}(t) ,
计算有效载荷向量 ( \hat{\mathbf{Q}}_t ):
- \hat{\mathbf{Q}}_t = \mathbf{Q}_t - (\mathbf{K} - c_2 \mathbf{M})\mathbf{u}_{curr} - (c_0 \mathbf{M} - c_1 \mathbf{C})\mathbf{u}_{prev} ,
- 注: \mathbf{K}\mathbf{u}_{curr} 项可通过单元内力累加计算，无需组装全局 \mathbf{K}
求解新位移 ( \mathbf{u}_{next} , 即 \mathbf{u}_{t+\Delta t} ):
- \hat{\mathbf{M}} \mathbf{u}_{next} = \hat{\mathbf{Q}}_t ,
- 若已分解，使用前代/回代求解；若为对角阵，直接逐元素相除。
计算当前时刻 ( t ) 的速度和加速度 (可选，用于输出):
- \ddot{\mathbf{u}}_t = c_0 (\mathbf{u}_{prev} - 2\mathbf{u}_{curr} + \mathbf{u}_{next}) ,
- \dot{\mathbf{u}}_t = c_1 (-\mathbf{u}_{prev} + \mathbf{u}_{next}) ,
- 操作: 输出或存储 t 时刻的 \mathbf{u}_{curr}, \dot{\mathbf{u}}_t, \ddot{\mathbf{u}}_t ,
更新变量以推进时间:
- \mathbf{u}_{prev} \leftarrow \mathbf{u}_{curr} ,
- \mathbf{u}_{curr} \leftarrow \mathbf{u}_{next} ,

3. 结束处理

输出最后一步结果:
- 计算 t_{end} = N_{steps} \times \Delta t 时刻的速度和加速度并输出。

关键说明

显式特性: 整个循环中，刚度矩阵 \mathbf{K} 从未出现在待求解矩阵的左侧，因此不需要对 \mathbf{K} 求逆或分解。
起步技巧: 通过公式 \mathbf{u}_{-\Delta t} = \mathbf{u}_0 - \Delta t \dot{\mathbf{u}}_0 + \frac{\Delta t^2}{2} \ddot{\mathbf{u}}_0 巧妙解决了中心差分法需要两个历史时间点才能启动的问题。
效率优化: 若使用集中质量矩阵（对角阵），步骤 4 退化为简单的标量除法，计算速度极快且内存占用极低。

中心差分法：算法推导

从泰勒级数展开出发，严格构建离散化格式，处理初始条件起步问题，并深入分析其稳定性准则。

1. 问题定义与控制方程

考虑一个 n 自由度线性动力学系统，其运动微分方程由下式给出：

\mathbf{M} \ddot{\mathbf{u}}(t) + \mathbf{C} \dot{\mathbf{u}}(t) + \mathbf{K} \mathbf{u}(t) = \mathbf{Q}(t) \quad \text{......(1)}

其中：

t ：时间变量。
\mathbf{u}(t) \in \mathbb{R}^n ：位移向量。
\dot{\mathbf{u}}(t) \in \mathbb{R}^n ：速度向量（位移的一阶导数）。
\ddot{\mathbf{u}}(t) \in \mathbb{R}^n ：加速度向量（位移的二阶导数）。
\mathbf{M} \in \mathbb{R}^{n \times n} ：质量矩阵（通常为正定对称矩阵）。
\mathbf{C} \in \mathbb{R}^{n \times n} ：阻尼矩阵。
\mathbf{K} \in \mathbb{R}^{n \times n} ：刚度矩阵。
\mathbf{Q}(t) \in \mathbb{R}^n ：外部载荷向量。

目标：已知 t=0 时刻的初始位移 \mathbf{u}_0 和初始速度 \dot{\mathbf{u}}_0 ，求解离散时间点 t_k = k \Delta t ( k=0, 1, 2, \dots ) 上的位移响应 \mathbf{u}_k 。其中 \Delta t 为时间步长。

2. 基于泰勒级数的差分近似推导

中心差分法的核心思想是利用当前时刻 t 前后两个时刻的位移值来近似该时刻的速度和加速度。

2.1 泰勒级数展开

将位移函数 \mathbf{u}(t) 在时刻 t 处分别向 t+\Delta t 和 t-\Delta t 进行泰勒级数展开：

前向展开 (Forward Expansion):

\mathbf{u}_{t+\Delta t} = \mathbf{u}_t + \Delta t \dot{\mathbf{u}}_t + \frac{\Delta t^2}{2!} \ddot{\mathbf{u}}_t + \frac{\Delta t^3}{3!} \dddot{\mathbf{u}}_t + O(\Delta t^4) \quad \text{......(2)}

后向展开 (Backward Expansion):

\mathbf{u}_{t-\Delta t} = \mathbf{u}_t - \Delta t \dot{\mathbf{u}}_t + \frac{\Delta t^2}{2!} \ddot{\mathbf{u}}_t - \frac{\Delta t^3}{3!} \dddot{\mathbf{u}}_t + O(\Delta t^4) \quad \text{......(3)}

其中 O(\Delta t^4) 表示截断误差项，其量级为 \Delta t^4 。

2.2 速度与加速度的中心差分格式

为了消去高阶项并解出 \dot{\mathbf{u}}_t 和 \ddot{\mathbf{u}}_t ，我们对上述两式进行线性组合。

推导速度公式： 将式 (2) 减去式 (3)

\mathbf{u}_{t+\Delta t} - \mathbf{u}_{t-\Delta t} = 2\Delta t \dot{\mathbf{u}}_t + \frac{2\Delta t^3}{3!} \dddot{\mathbf{u}}_t + O(\Delta t^5)

忽略 O(\Delta t^3) 及更高阶小量，整理得 t 时刻速度的中心差分近似：

\dot{\mathbf{u}}_t \approx \frac{1}{2\Delta t} (\mathbf{u}_{t+\Delta t} - \mathbf{u}_{t-\Delta t}) \quad \text{......(4)}

注：该近似的截断误差为 O(\Delta t^2) 。 推导加速度公式： 将式 (2) 加上式 (3)：

\mathbf{u}_{t+\Delta t} + \mathbf{u}_{t-\Delta t} = 2\mathbf{u}_t + \Delta t^2 \ddot{\mathbf{u}}_t + O(\Delta t^4)

忽略 O(\Delta t^4) 项，整理得 t 时刻加速度的中心差分近似：

\ddot{\mathbf{u}}_t \approx \frac{1}{\Delta t^2} (\mathbf{u}_{t+\Delta t} - 2\mathbf{u}_t + \mathbf{u}_{t-\Delta t}) \quad \text{......(5)}

注：该近似的截断误差亦为 O(\Delta t^2) 。

3. 运动方程的离散化与递推公式

将推导出的近似公式 (4) 和 (5) 代入原始运动方程 (1) 中，且在 t 时刻建立平衡关系:

\mathbf{M} \left[ \frac{1}{\Delta t^2} (\mathbf{u}_{t+\Delta t} - 2\mathbf{u}_t + \mathbf{u}_{t-\Delta t}) \right] + \mathbf{C} \left[ \frac{1}{2\Delta t} (\mathbf{u}_{t+\Delta t} - \mathbf{u}_{t-\Delta t}) \right] + \mathbf{K} \mathbf{u}_t = \mathbf{Q}_t \quad \text{......(6)}

我们的目标是求解未知量 \mathbf{u}_{t+\Delta t} 。因此，需要将含有 \mathbf{u}_{t+\Delta t} 的项移至等式左侧，其余已知项（ \mathbf{u}_t, \mathbf{u}_{t-\Delta t} ）移至右侧。

3.1 合并同类项

提取 \mathbf{u}_{t+\Delta t} 的系数：

\left( \frac{1}{\Delta t^2} \mathbf{M} + \frac{1}{2\Delta t} \mathbf{C} \right) \mathbf{u}_{t+\Delta t} = \mathbf{Q}_t - \mathbf{K} \mathbf{u}_t + \mathbf{M} \left( \frac{2}{\Delta t^2} \mathbf{u}_t - \frac{1}{\Delta t^2} \mathbf{u}_{t-\Delta t} \right) + \mathbf{C} \left( \frac{1}{2\Delta t} \mathbf{u}_{t-\Delta t} \right)

整理右侧关于 \mathbf{u}_t 和 \mathbf{u}_{t-\Delta t} 的项：

\begin{aligned} \text{Right Hand Side} &= \mathbf{Q}_t - \left( \mathbf{K} - \frac{2}{\Delta t^2} \mathbf{M} \right) \mathbf{u}_t - \left( \frac{1}{\Delta t^2} \mathbf{M} - \frac{1}{2\Delta t} \mathbf{C} \right) \mathbf{u}_{t-\Delta t} \end{aligned}

3.2 定义有效矩阵与载荷向量

为了简化表达，引入以下记号（对应 P1 中的积分常数）：令 c_0 = \frac{1}{\Delta t^2}, \quad c_1 = \frac{1}{2\Delta t} 。则 c_2 = 2c_0 = \frac{2}{\Delta t^2} 。

定义 有效质量矩阵 (Effective Mass Matrix) \hat{\mathbf{M}} ：

\hat{\mathbf{M}} = c_0 \mathbf{M} + c_1 \mathbf{C} = \frac{1}{\Delta t^2} \mathbf{M} + \frac{1}{2\Delta t} \mathbf{C} \quad \text{......(7)}

定义 ** t 时刻的有效载荷向量 (Effective Load Vector)** \hat{\mathbf{Q}}_t ：

\hat{\mathbf{Q}}_t = \mathbf{Q}_t - \left( \mathbf{K} - c_2 \mathbf{M} \right) \mathbf{u}_t - \left( c_0 \mathbf{M} - c_1 \mathbf{C} \right) \mathbf{u}_{t-\Delta t} \quad \text{......(8)}

于是，离散化的运动方程简化为线性方程组：

\hat{\mathbf{M}} \mathbf{u}_{t+\Delta t} = \hat{\mathbf{Q}}_t \quad \text{......(9)}

3.3 显式积分特性分析

观察式 (9) 左端的 \hat{\mathbf{M}} ：

刚度矩阵缺席： \mathbf{K} 仅出现在右端载荷项中，不参与左端系数矩阵的构成。这意味着不需要对 \mathbf{K} 进行三角分解或求逆。
对角化优势：若采用 集中质量矩阵 (Lumped Mass Matrix)，则 \mathbf{M} 为对角阵。若阻尼 \mathbf{C} 可忽略或也为对角阵（如瑞利阻尼中的质量比例部分），则 \hat{\mathbf{M}} 也是对角阵。此时，式 (9) 退化为 n 个独立的标量方程，第 i 个自由度的解为：u_{t+\Delta t}^{(i)} = \frac{\hat{Q}_t^{(i)}}{\hat{M}_{ii}} = \frac{\hat{Q}_t^{(i)}}{c_0 M_{ii} + c_1 C_{ii}} \quad \text{......(10)}这避免了大型矩阵求逆运算，计算效率极高，称为 显式积分 (Explicit Integration)。

4. 起步程序 (Starting Procedure)

中心差分法是两步法（Two-step method），计算 \mathbf{u}_{t+\Delta t} 需要依赖 \mathbf{u}_t 和 \mathbf{u}_{t-\Delta t} 。在初始时刻 t=0 ，已知 \mathbf{u}_0 和 \dot{\mathbf{u}}_0 ，但 \mathbf{u}_{-\Delta t} 未知。必须构造虚拟的起始位移。

4.1 初始加速度计算

首先利用 t=0 时刻的运动方程 (1) 求出初始加速度 \ddot{\mathbf{u}}_0 ：

\mathbf{M} \ddot{\mathbf{u}}_0 = \mathbf{Q}_0 - \mathbf{C} \dot{\mathbf{u}}_0 - \mathbf{K} \mathbf{u}_0 \implies \ddot{\mathbf{u}}_0 = \mathbf{M}^{-1} (\mathbf{Q}_0 - \mathbf{C} \dot{\mathbf{u}}_0 - \mathbf{K} \mathbf{u}_0) \quad \text{......(11)}

4.2 反向泰勒展开推导 \mathbf{u}_{-\Delta t}

利用后差分公式 (3) 的截断形式，但在 t=0 处展开至二阶精度以保证启动精度：

\mathbf{u}_{-\Delta t} = \mathbf{u}_0 - \Delta t \dot{\mathbf{u}}_0 + \frac{\Delta t^2}{2} \ddot{\mathbf{u}}_0 + O(\Delta t^3)

忽略高阶项，得到起步公式：

\mathbf{u}_{-\Delta t} = \mathbf{u}_0 - \Delta t \dot{\mathbf{u}}_0 + \frac{\Delta t^2}{2} \ddot{\mathbf{u}}_0 \quad \text{......(12)}

注：P1 中使用积分常数 c_3 = 1/c_2 = \Delta t^2/2 ，故写作 \mathbf{u}_{-\Delta t} = \mathbf{u}_0 - \Delta t \dot{\mathbf{u}}_0 + c_3 \ddot{\mathbf{u}}_0 。

至此， \mathbf{u}_{-\Delta t} 已求得，算法可以进入 t=0 步，求解 \mathbf{u}_{\Delta t} 。

5. 稳定性分析 (Stability Analysis)

中心差分法是 条件稳定 (Conditionally Stable) 的。其稳定性取决于时间步长 \Delta t 与系统最高频率的关系。

5.1 单自由度无阻尼系统模型

考虑单自由度无阻尼自由振动方程：

\ddot{u} + \omega^2 u = 0

其中 \omega 为固有圆频率。将中心差分格式 (5) 代入（此时 \mathbf{C}=0, \mathbf{Q}=0 ）：

\frac{u_{t+\Delta t} - 2u_t + u_{t-\Delta t}}{\Delta t^2} + \omega^2 u_t = 0

整理得递推关系：

u_{t+\Delta t} - (2 - \omega^2 \Delta t^2) u_t + u_{t-\Delta t} = 0 \quad \text{......(13)}

5.2 特征方程与谱半径

设解的形式为 u_t = \lambda^t u_0 （或更严谨地，写成状态空间形式 \mathbf{X}_{t+\Delta t} = \mathbf{D} \mathbf{X}_t ，其中 \mathbf{X}_t = [u_t, u_{t-\Delta t}]^T ）。代入式 (13) 得到特征方程：

\lambda^2 - (2 - \omega^2 \Delta t^2) \lambda + 1 = 0 \quad \text{......(14)}

该二次方程的判别式为：

\Delta = (2 - \omega^2 \Delta t^2)^2 - 4

根据 P1 中的稳定性理论，算法稳定的充要条件是递推矩阵 \mathbf{D} 的谱半径 \rho(\mathbf{D}) = \max |\lambda_i| \le 1 。对于重根情况需严格小于 1，对于非重根可等于 1。

情形讨论：

**若 (2 - \omega^2 \Delta t^2)^2 < 4 **：此时 \Delta < 0 ，特征根为一对共轭复数。由韦达定理，两根之积 \lambda_1 \lambda_2 = 1 。因为是共轭复数，模长相等，故 |\lambda_1| = |\lambda_2| = 1 。此时算法稳定，且无数值阻尼（能量守恒）。条件等价于：
-2 < 2 - \omega^2 \Delta t^2 < 2
右不等式 2 - \omega^2 \Delta t^2 < 2 \implies -\omega^2 \Delta t^2 < 0 恒成立。左不等式 -2 < 2 - \omega^2 \Delta t^2 \implies \omega^2 \Delta t^2 < 4 \implies \omega \Delta t < 2 。
若 (2 - \omega^2 \Delta t^2)^2 \ge 4 **：此时 \Delta \ge 0 ，特征根为实数。由于两根之积为 1，必有一个根的模 |\lambda| \ge 1 。若 \Delta > 0 ，则必有一根 |\lambda| > 1 ，导致解随时间指数增长，算法 **不稳定。

5.3 临界时间步长

综上所述，稳定性条件为：

\omega \Delta t \le 2

即：

\Delta t \le \frac{2}{\omega} = \frac{2}{2\pi/T} = \frac{T}{\pi} \quad \text{......(15)}

其中 T 为系统周期。

对于多自由度系统，稳定性由 最高频率 \omega_{\max} （或最小周期 T_{\min} ）控制：

\Delta t \le \Delta t_{cr} = \frac{2}{\omega_{\max}} = \frac{T_{\min}}{\pi} \quad \text{......(16)}

结论：中心差分法没有数值阻尼（在稳定区内谱半径恒为 1），且时间步长必须小于临界步长 \Delta t_{cr} 。如果网格划分导致某些单元极小， \omega_{\max} 会极大，从而迫使 \Delta t 极小，显著增加计算成本。

6. 常用变异格式 (Variants)

P1 中还提到了几种基于中心差分思想的变体，主要用于不同场景（如分子动力学）：

6.1 Verlet 格式

直接利用位移展开式相加的结果：

\mathbf{u}_{t+\Delta t} = 2\mathbf{u}_t - \mathbf{u}_{t-\Delta t} + \ddot{\mathbf{u}}_t \Delta t^2

特点：位移精度 O(\Delta t^4) ，但不直接输出速度。速度需事后通过 \dot{\mathbf{u}}_t \approx (\mathbf{u}_{t+\Delta t} - \mathbf{u}_{t-\Delta t}) / 2\Delta t 计算，精度降为 O(\Delta t^2) 。

6.2 蛙跳格式 (Leapfrog)

速度与位移在时间轴上交错半步：

\begin{aligned} \mathbf{v}_{t+\frac{1}{2}\Delta t} &= \mathbf{v}_{t-\frac{1}{2}\Delta t} + \ddot{\mathbf{u}}_t \Delta t \\ \mathbf{u}_{t+\Delta t} &= \mathbf{u}_t + \mathbf{v}_{t+\frac{1}{2}\Delta t} \Delta t \end{aligned}

特点：同样不直接给出同一时刻的速度和位移。

6.3 速度 Verlet 格式 (Velocity Verlet)

结合了 Leapfrog 的优点，能在同一时刻获得位移、速度和加速度，步骤如下：

1: \mathbf{v}_{t+\frac{1}{2}\Delta t} = \mathbf{v}_t + \frac{1}{2} \ddot{\mathbf{u}}_t \Delta t ,

2: \mathbf{u}_{t+\Delta t} = \mathbf{u}_t + \mathbf{v}_{t+\frac{1}{2}\Delta t} \Delta t ,

3: 计算新加速度 \ddot{\mathbf{u}}_{t+\Delta t} = \mathbf{M}^{-1}(\mathbf{Q}_{t+\Delta t} - \mathbf{K}\mathbf{u}_{t+\Delta t}) ,

4: \mathbf{v}_{t+\Delta t} = \mathbf{v}_{t+\frac{1}{2}\Delta t} + \frac{1}{2} \ddot{\mathbf{u}}_{t+\Delta t} \Delta t ,

此格式广泛应用于 LS-DYNA 等显式有限元软件及分子动力学模拟。该方法是解决短时、高频、强非线性（如冲击、碰撞）动力学问题的首选数值工具。

返回主页

目录

中心差分法：简要总结

1. 核心原理与公式推导

2. 算法特性：显式积分 (Explicit Integration)

3. 起步程序 (Starting Procedure)

4. 常用变异格式

5. 稳定性分析 (Stability Analysis)

6. 实际应用中的限制与注意事项

7. 算法流程总结

中心差分法：算法流程

1. 初始化阶段

2. 时间步循环迭代

3. 结束处理

关键说明

中心差分法：算法推导

1. 问题定义与控制方程

2. 基于泰勒级数的差分近似推导

2.1 泰勒级数展开

2.2 速度与加速度的中心差分格式

3. 运动方程的离散化与递推公式

3.1 合并同类项

3.2 定义有效矩阵与载荷向量

3.3 显式积分特性分析

4. 起步程序 (Starting Procedure)

4.1 初始加速度计算

4.2 反向泰勒展开推导 \mathbf{u}_{-\Delta t}

5. 稳定性分析 (Stability Analysis)

5.1 单自由度无阻尼系统模型

5.2 特征方程与谱半径

5.3 临界时间步长

6. 常用变异格式 (Variants)

6.1 Verlet 格式

6.2 蛙跳格式 (Leapfrog)

6.3 速度 Verlet 格式 (Velocity Verlet)