基于一阶后向欧拉法(BDF1)求解非线性动力学系统的牛顿迭代推导

在结构动力学和有限元分析中，求解二阶非线性常微分方程是常见需求。本文将详细推导利用 一阶后向欧拉法（BDF1/Backward Euler） 配合 直接牛顿-拉夫逊（Newton-Raphson）迭代法 求解该类动力学系统的完整公式推导过程。

---

1. 动力学控制方程与状态降阶

考虑一般的二阶非线性动力学平衡方程：
$$M\ddot{u}+K(u)u=F$$

其中：

• $M$ 为质量矩阵；
• $K(u)$ 为依赖于位移的非线性刚度矩阵；
• $u,\dot{u},\ddot{u}$ 分别为位移、速度与加速度向量；
• $F$ 为外载荷向量。

为了使用一阶隐式积分算法（后向欧拉），引入速度变量 $v=\dot{u}$，将二阶系统降阶为一阶耦合微分方程组：
$$\{\begin{array}{l}\dot{u}=v\\ M\dot{v}+K(u)u=F\end{array}$$

---

2. 时间离散化（后向欧拉法 BDF1）

设时间步长为 $\Delta t$（简记为 $dt$），在离散时刻 $t_{n+1}$，采用一阶后向差分近似时间导数：

• 速度与位移的关系：
  $${\dot{u}}_{n+1}\approx \frac{u_{n+1}-u_n}{dt}=v_{n+1}$$
• 加速度与速度的关系：
  $${\dot{v}}_{n+1}\approx \frac{v_{n+1}-v_n}{dt}$$

将其代入动力学方程组，在 $t_{n+1}$ 时刻的理想平衡状态应满足：
$$\{\begin{array}{l}\frac{u_{n+1}-u_n}{dt}-v_{n+1}=0\\ \frac{M}{dt}(v_{n+1}-v_n)+K(u_{n+1})u_{n+1}-F_{n+1}=0\end{array}$$

---

3. 构造残差方程（Residuals）

在非线性求解过程中，第 $k$ 次牛顿迭代给出的状态估计值记为 $(u_{n+1}^k,v_{n+1}^k)$。此时系统尚未完全平衡，定义如下残差向量：

3.1 运动学位移残差

$$R_u=\frac{u_{n+1}^k-u_n}{dt}-v_{n+1}^k$$

3.2 动力学平衡残差

$$R_v=\frac{M}{dt}(v_{n+1}^k-v_n)+K(u_{n+1}^k)u_{n+1}^k-F_{n+1}$$

(当残差范数 $\Vert R_u\Vert \to 0$ 且 $\Vert R_v\Vert \to 0$ 时，方程收敛)

---

4. 多元牛顿法一阶泰勒展开（线性化）

寻求增量步 $\Delta u$ 和 $\Delta v$，令更新后的残差近似为零：
$$\{\begin{array}{l}{R}_u+\frac{\partial R_u}{\partial u}\Delta u+\frac{\partial R_u}{\partial v}\Delta v=0\\ R_v+\frac{\partial R_v}{\partial u}\Delta u+\frac{\partial R_v}{\partial v}\Delta v=0\end{array}$$

计算 Jacobian 偏导数矩阵：

1. 对 $R_u$ 求导：
   $$\frac{\partial R_u}{\partial u}=\frac{1}{dt},\frac{\partial R_u}{\partial v}=-I$$
2. 对 $R_v$ 求导：
   $$\frac{\partial R_v}{\partial u}=\frac{\partial (K(u)u)}{\partial u}=K(u_{n+1}^k)+\frac{\partial K}{\partial u}{u}_{n+1}^k\triangleq K_T(\text{切线刚度矩阵 Tangent Stiffness})$$
   $$\frac{\partial R_v}{\partial v}=\frac{M}{dt}$$

代入得到线性方程组：
$$\{\begin{array}{ll}\frac{1}{dt}\Delta u-\Delta v=-R_u& \cdots (1)\\ K_T\Delta u+\frac{M}{dt}\Delta v=-R_v& \cdots (2)\end{array}$$

---

5. 变量消元与核心方程推导

由方程 $(1)$ 可将速度增量 $\Delta v$ 显式表达为位移增量 $\Delta u$ 的形式：
$$\Delta v=\frac{\Delta u}{dt}+R_u\cdots (3)$$

将方程 $(3)$ 代入方程 $(2)$ 中：
$$K_T\Delta u+\frac{M}{dt}\left(\frac{\Delta u}{dt}+R_u\right)=-R_v$$

展开整理：
$$K_T\Delta u+\frac{M}{d{t}^2}\Delta u+\frac{M}{dt}{R}_u=-R_v$$

合并同类项，得到核心求解方程：
$$\left(K_T+\frac{M}{d{t}^2}\right)\Delta u=-R_v-\frac{M}{dt}{R}_u$$

核心符号定义：

• 等效切线刚度矩阵：
  $$\hat{K}=K_T+\frac{M}{d{t}^2}$$
• 等效右端残差项：
  $$\hat{R}=-R_v-\frac{M}{dt}{R}_u$$

---

6. 状态更新与算法流程

算法步骤总结：

1. 求解增量：
   $$\Delta u={\hat{K}}^{-1}\hat{R}$$
2. 更新位移：
   $$u_{n+1}^{k+1}=u_{n+1}^k+\Delta u$$
3. 更新速度：
   $$v_{n+1}^{k+1}=v_{n+1}^k+\frac{\Delta u}{dt}+R_u$$
4. 收敛判断：
   检查 $\Vert \hat{R}\Vert \le \text{tol}$ 或 $\Vert \Delta u\Vert \le \text{tol}$。若未达到收敛准则，令 $k=k+1$ 重复上述步骤；若收敛，则进入下一个时间步 $t_{n+2}$。