弹簧系统 Spring System

使用变分力学分析弹簧质点系统

单弹簧

假设存在空间中的两个质点

$$\mathbf{x}_0 = \begin{bmatrix} x_0 \\ y_0 \\ z_0 \end{bmatrix} \quad \mathbf{x}_1 = \begin{bmatrix} x_1 \\ y_1 \\ z_1 \end{bmatrix}$$

它们由静止长度为 $l_0$ 的弹簧相连

显然，这个 3 维弹簧物理系统与两个质点的坐标有关，所以 广义坐标 是

$$\mathbf{q} = \begin{bmatrix} \mathbf{x}_0 \\ \mathbf{x}_1 \end{bmatrix}$$

它是一个 6 维向量，我们求它的时间导数得到 广义速度

$$\mathbf{\dot{q}} = \begin{bmatrix} \mathbf{\dot{x}}_0 \\ \mathbf{\dot{x}}_1 \end{bmatrix} = \begin{bmatrix} \mathbf{v}_0 \\ \mathbf{v}_1 \end{bmatrix}$$

我们想要构建这个系统的 拉格朗日量

$$L = T - V$$

所以我们要求系统的 动能 $T$ 和 势能 $V$

系统的 动能 是两个质点的动能之和

$$\begin{align*} T &= \sum_{i=0}^{1}\frac{1}{2}m_i \| \mathbf{v}_i \|_2^2 = \sum_{i=0}^{1}\frac{1}{2}m_i \mathbf{v}_i^T\mathbf{v}_i \\ &= \sum_{i=0}^{1}\frac{1}{2} \mathbf{v}_i^T \underbrace{ \begin{bmatrix} m_x & 0 & 0 \\ 0 & m_y & 0 \\ 0 & 0 & m_z \end{bmatrix}_i }_{\mathbf{M}_i} \mathbf{v}_i \end{align*}$$

• 其中速度向量的二范数平方可以表示为速度向量与自身的内积
• 其中质点的质量可以写成矩阵形式，质点的质量可以是 各向异性 的，即在不同方向上有不同的质量分量

我们可以使用刚刚计算的 广义速度 进一步简化动能等式

$$\begin{align*} \mathbf{M} &= \begin{bmatrix} \mathbf{M}_0 & 0 \\ 0 & \mathbf{M}_1 \end{bmatrix} \\ T &= \frac{1}{2}\mathbf{\dot{q}}^T\mathbf{M}\mathbf{\dot{q}} \end{align*}$$

这利用 广义速度 将两个质点的速度写成一个向量的紧凑形式，将两个质点的质量也写为一个紧凑的质量对角矩阵

系统的 势能 是存储在弹簧中的 弹性势能，弹性势能 以某种形式和弹簧的形变程度相关，总的来说，它需要满足以下条件

1. 在没有外力的情况下，应当趋于其初始状态，例如弹簧在不受外力的情况下应该倾向于恢复到静止长度，也就是说能量的最低值应该由静止状态 (未形变状态) 定义
2. 刚性变换，例如位移、旋转不应该改变 弹性势能
3. 弹性势能 应当仅依赖于质点位置

我们首先将弹簧的 形变程度 定义弹簧的 应变 Strain

$$\text{Strain} = l - l_0$$

• $l$ 是弹簧的形变长度，$l_0$ 是弹簧的静止长度

然后我们可以根据 应变 定义弹簧的 弹性势能 为

$$\text{Potential Energy} = \frac{1}{2}k(l-l_0)^2$$

• $k$ 是弹簧的 刚度参数 Stiffness Parameter

我们需要将弹簧的形变长度 $l$ 定义为质点的坐标形式，首先计算 $\Delta \mathbf{x}$

$$\begin{align*} \Delta \mathbf{x} &= \mathbf{x}_1 - \mathbf{x}_0 \\ &= \begin{bmatrix} -\mathbf{I} & \mathbf{I} \end{bmatrix} \begin{bmatrix} \mathbf{x}_0 \\ \mathbf{x}_1 \end{bmatrix} \end{align*}$$

我们可以看到 $\begin{bmatrix} \mathbf{x}_0 \\ \mathbf{x}_1 \end{bmatrix}$ 实际上是我们之前定义的 广义坐标 $\mathbf{q}$，若我们做如下定义

$$\mathbf{B} = \begin{bmatrix} -\mathbf{I} & \mathbf{I} \end{bmatrix} \quad \mathbf{q} = \begin{bmatrix} \mathbf{x}_0 \\ \mathbf{x}_1 \end{bmatrix}$$

那么可以将 $\Delta \mathbf{x}$ 写成由广义坐标定义

$$\Delta \mathbf{x} = \mathbf{B}\mathbf{q}$$

$\Delta \mathbf{x}$ 向量的模长就是弹簧的形变长度 $l$

$$l = \sqrt{\Delta \mathbf{x}^T \Delta \mathbf{x}} = \sqrt{\mathbf{q}^T\mathbf{B}^T\mathbf{B}\mathbf{q}}$$

我们将这个广义坐标形式的形变长度代入，得到最终的弹性势能函数

$$V = \frac{1}{2}k(\sqrt{\mathbf{q}^T\mathbf{B}^T\mathbf{B}\mathbf{q}} - l_0)^2$$

由此，我们可以得到单弹簧质点系统的 拉格朗日量

$$\begin{align*} L &= T - V \\ &= \frac{1}{2}\mathbf{\dot{q}}^T\mathbf{M}\mathbf{\dot{q}} - \frac{1}{2}k(\sqrt{\mathbf{q}^T\mathbf{B}^T\mathbf{B}\mathbf{q}} - l_0)^2 \end{align*}$$

接着，我们可以考虑 欧拉-拉格朗日方程 Euler-Lagrange Equation

$$\frac{d}{dt} \frac{\partial L}{\partial\mathbf{\dot{q}}} = \frac{\partial L}{\partial\mathbf{q}}$$

我们注意到，拉格朗日量 对 广义坐标 的偏导数仅和 势能函数 有关

$$\begin{align*} \frac{\partial L}{\partial\mathbf{q}} &= \frac{\partial}{\partial\mathbf{q}} (T - V) \\ &= \frac{\partial T}{\partial\mathbf{q}} - \frac{\partial V}{\partial\mathbf{q}} \\ &= 0 - \frac{\partial V}{\partial\mathbf{q}} \\ &= -\frac{\partial V}{\partial\mathbf{q}} \end{align*}$$

代回 欧拉-拉格朗日方程，可得

$$\frac{d}{dt} \frac{\partial L}{\partial\mathbf{\dot{q}}} = -\frac{\partial V}{\partial\mathbf{q}}$$

我们将 $-\dfrac{\partial V}{\partial\mathbf{q}}$ 定义为 广义力 Generalized Forces，也就是说力 $\mathbf{f}$ 是势能场 $V$ 的负梯度，这表示力的方向总是朝向势能下降的方向，例如，在重力场中，物体会自发地沿着势能减小的方向移动

$$\mathbf{f} = -\nabla V = -\frac{\partial V}{\partial\mathbf{q}}$$

我们接着分析等式的左边

$$\frac{d}{dt}\frac{\partial L}{\partial\mathbf{\dot{q}}} = \frac{d}{dt}\frac{\partial}{\partial\mathbf{\dot{q}}}\left(\frac{1}{2}\mathbf{\dot{q}}^T\mathbf{M}\mathbf{\dot{q}} \right)$$

其中 $\frac{1}{2}\mathbf{\dot{q}}^T\mathbf{M}\mathbf{\dot{q}}$ 可以认为是 $\mathbf{\dot{q}}$ 的二次项，所以有

$$\begin{align*} \frac{d}{dt}\frac{\partial}{\partial\mathbf{\dot{q}}}\left(\frac{1}{2}\mathbf{\dot{q}}^T\mathbf{M}\mathbf{\dot{q}} \right) &= \frac{d}{dt}\mathbf{M}\mathbf{\dot{q}} \\ &= \mathbf{M}\mathbf{\ddot{q}} \end{align*}$$

将结果代回 欧拉-拉格朗日方程，我们获得 Equation of Motion 运动方程

$$\mathbf{M}\mathbf{\ddot{q}} = -\frac{\partial V}{\partial\mathbf{q}}$$

我们通过 时间积分 对这个方程进行仿真

多弹簧系统

我们可以将物体离散为多个弹簧质点组成的整体系统 假设系统中有 $n$ 个 质点，那么我们的广义坐标是

$$\mathbf{q} = \begin{bmatrix} \mathbf{x}_0 \\ \mathbf{x}_1 \\ \mathbf{x}_2 \\ \vdots \\ \mathbf{x}_n \end{bmatrix}$$

若其中每个 $\mathbf{x}_{\cdots}$ 都是一个三维坐标，那么 $\mathbf{q} \in \mathbb{R}^{3n \times 1}$

同样的，我们可以得出广义速度

$$\mathbf{\dot{q}} = \begin{bmatrix} \mathbf{v}_0 \\ \mathbf{v}_1 \\ \mathbf{v}_2 \\ \vdots \\ \mathbf{v}_n \end{bmatrix} \in \mathbb{R}^{3n \times 1}$$

接着，我们计算系统的 动能，系统的总动能是所有质点的动能之和

$$T = \sum_{i=0}^{n - 1}\frac{1}{2} \mathbf{v}_i^T \underbrace{ \begin{bmatrix} m_x & 0 & 0 \\ 0 & m_y & 0 \\ 0 & 0 & m_z \end{bmatrix}_i }_{\mathbf{M}_i} \mathbf{v}_i$$

和之前一样，我们可以通过使用 广义速度 以及将所有质点质量矩阵放到一个巨大的对角矩阵内，来简化掉求和符号

$$\begin{align*} \mathbf{M} &= \begin{bmatrix} \mathbf{M}_0 & \cdots & 0 \\ \vdots & \ddots &\vdots \\ 0 & \cdots & \mathbf{M}_n \end{bmatrix} \in \mathbb{R}^{3n \times 3n} \\ T &= \frac{1}{2}\mathbf{\dot{q}}^T\mathbf{M}\mathbf{\dot{q}} \end{align*}$$

与计算 动能 类似，系统的总 势能 是所有弹簧的势能之和，假设系统中有 $m$ 个弹簧

$$V = \sum_{j=0}^{m-1} V_j(\mathbf{x}_A,\mathbf{x}_B)$$

由于每个弹簧的弹性势能都和其端点的两个质点坐标 $\mathbf{x}_A,\mathbf{x}_B$ 相关，所以我们需要从广义坐标 $\mathbf{q}$ 中选出当前弹簧所对应的端点坐标，我们引入 坐标选择矩阵

$$\begin{bmatrix} \mathbf{A}_0 & \mathbf{A}_1 & \cdots & \mathbf{A}_n \end{bmatrix} \begin{bmatrix} \mathbf{x}_0 \\ \mathbf{x}_1 \\ \mathbf{x}_2 \\ \vdots \\ \mathbf{x}_n \end{bmatrix} = \sum_{i=1}^{n}\mathbf{A}_i\mathbf{x}_i$$

若我们想要选取第 $i$ 个质点的坐标，那么只需要将 $\mathbf{A}_i = \mathbf{I}$，其余所有都设置为零 $\mathbf{A}_{j\neq i} = \mathbf{0}$，例如选择第二个质点坐标 $\mathbf{x}_1$ 是

$$\underbrace{ \begin{bmatrix} \mathbf{0} & \mathbf{I} & \cdots & \mathbf{0} \end{bmatrix} }_{\mathbf{S}_1} \begin{bmatrix} \mathbf{x}_0 \\ \mathbf{x}_1 \\ \mathbf{x}_2 \\ \vdots \\ \mathbf{x}_n \end{bmatrix} = \mathbf{x}_1$$

我们将 坐标选择矩阵 记为 $\mathbf{S}$，有了它，我们就可以从完整的 广义坐标 $\mathbf{q}$ 中选择出我们想要的第 $j$ 个弹簧的两端质点坐标 $\mathbf{x}_A$ 和 $\mathbf{x}_B$ 组成的广义坐标 $\mathbf{q}_j$

$$\mathbf{q}_j = \begin{bmatrix} \mathbf{x}_A \\ \mathbf{x}_B \end{bmatrix} = \underbrace{ \begin{bmatrix} \mathbf{S}_A \\ \mathbf{S}_B \end{bmatrix} }_{\mathbf{E}_j} \mathbf{q}$$

我们将 $\mathbf{S}_A$ 和 $\mathbf{S}_B$ 组成的 坐标选择矩阵 记为 $\mathbf{E}_j$ 所以原总势能也变为可以由广义坐标来表述

$$V = \sum_{j=0}^{m-1} V_j(\mathbf{E}_j\mathbf{q})$$

然后，我们可以构建多弹簧质点系统的 拉格朗日量

$$L = T - V = \frac{1}{2}\mathbf{\dot{q}}^T\mathbf{M}\mathbf{\dot{q}} - \sum_{j=0}^{m-1} V_j(\mathbf{E}_j\mathbf{q})$$

然后我们分析 广义力 形式的 欧拉-拉格朗日方程

$$\frac{d}{dt} \frac{\partial L}{\partial\mathbf{\dot{q}}} = -\frac{\partial V}{\partial\mathbf{q}}$$

我们首先分析等式的左边

$$\begin{align*} \frac{d}{dt}\frac{\partial L}{\partial\mathbf{\dot{q}}} &= \frac{d}{dt}\frac{\partial}{\partial\mathbf{\dot{q}}}\left(\frac{1}{2}\mathbf{\dot{q}}^T\mathbf{M}\mathbf{\dot{q}} \right) \\ &= \frac{d}{dt}\mathbf{M}\mathbf{\dot{q}} \\ &= \mathbf{M}\mathbf{\ddot{q}} \end{align*}$$

可以发现这一部分的形式没有改变，我们同样得到 运动方程

$$\mathbf{M}\mathbf{\ddot{q}} = -\frac{\partial V}{\partial\mathbf{q}}$$

当然，这次的 质量矩阵 $\mathbf{M}$ 和 广义坐标 $\mathbf{q}$ 包含了多个弹簧质点，所以要比之前大很多

我们接着分析等式右边的 广义力 项

$$-\frac{\partial V}{\partial\mathbf{q}} = -\frac{\partial}{\partial\mathbf{q}} \sum_{j=0}^{m-1} V_j(\mathbf{E}_j\mathbf{q})$$

由于梯度算子 $\frac{\partial}{\partial\mathbf{q}}$ 是线性算子，所以可以将它移入求和

$$\begin{align*} -\frac{\partial}{\partial\mathbf{q}} \sum_{j=0}^{m-1} V_j(\mathbf{E}_j\mathbf{q}) &= -\sum_{j=0}^{m-1} \frac{\partial}{\partial\mathbf{q}} V_j(\mathbf{E}_j\mathbf{q}) \\ &= -\sum_{j=0}^{m-1} \mathbf{E}^T_j \frac{\partial V_j}{\partial\mathbf{q}_j} (\mathbf{q}_j) \end{align*}$$

可以发现我们需要求每个弹簧的 弹性势能函数 $V_j$ 对其两端点质点广义坐标 $\mathbf{q}_j$ 的梯度，其中 弹性势能函数 一般和单个弹簧定义相同，即

$$V = \frac{1}{2}k(\sqrt{\mathbf{q}^T\mathbf{B}^T\mathbf{B}\mathbf{q}} - l_0)^2$$

可以看到，每个 弹性势能函数 对其两端点 广义坐标 的偏导数 (梯度) 就是该弹簧产生的 广义力，所以我们也可以写为如下形式

$$\begin{align*} \frac{\partial V_j}{\partial\mathbf{q}_j} (\mathbf{q}_j) &= \mathbf{f}_j (\mathbf{q}_j) \\[3mm] -\frac{\partial V}{\partial\mathbf{q}} = -\sum_{j=0}^{m-1} \mathbf{E}^T_j \frac{\partial V_j}{\partial\mathbf{q}_j} (\mathbf{q}_j) &= -\sum_{j=0}^{m-1} \mathbf{E}^T_j \mathbf{f}_j (\mathbf{q}_j) \end{align*}$$

其中 $\mathbf{f}_j (\mathbf{q}_j)$ 是第 $j$ 个弹簧产生的 广义力，它是弹簧两端质点的力组成的向量，在三维情况下，它是

$$\mathbf{f}_j = \begin{bmatrix} \mathbf{f}_A \\ \mathbf{f}_B \end{bmatrix} \in \mathbb{R}^6$$

它通过对 第 $j$ 个弹簧的弹性势能求对 广义坐标 的偏导 (梯度) 得到

$$\frac{\partial V_j}{\partial\mathbf{q}_j} (\mathbf{q}_j) = \mathbf{f}_j (\mathbf{q}_j)$$

然后我们通过对应的 第 $j$ 个 选择矩阵 的 转置 $\mathbf{E}^T_j$ 将 $\mathbf{f}_j$ 放到最终大型 广义力 向量 $\mathbf{f}$ 的对应位置 (每个弹簧的力只作用在广义坐标的一部分分量上)，也就是说，这一步得到是多个只在某些分量有值的向量，求和以后才是每个分量都有值的广义力向量

$$-\sum_{j=0}^{m-1} \mathbf{E}^T_j \mathbf{f}_j (\mathbf{q}_j) = \mathbf{f} = \begin{bmatrix} \mathbf{f}_0 \\ \mathbf{f}_1 \\ \vdots \\ \mathbf{f}_n \end{bmatrix} \in \mathbb{R}^{3n \times 1}$$

• 注意负号，这里的负号是必要的，因为我们求的是势能的负梯度，代表的是从势能中导出的广义力，力指向势能减少的方向

注意，每个质点可能连接多个弹簧，所以质点最终的力是所有连接到该质点的弹簧对其产生的合力，也就是说，每个弹簧都会对不同质点产生分力，即对不同的坐标分量有贡献，而合力是每个弹簧力对不同分量的总贡献

有了 广义力 向量，我们就可以进行 时间积分 了，并且我们导出了最终 运动方程

$$\mathbf{M}\mathbf{\ddot{q}} = \mathbf{f}$$

我们发现这实际上就是牛顿的 $f = ma$，根据 拉格朗日-欧拉方程 推导出 运动方程 以后，我们可以通过计算势能函数的梯度得到 广义力 $\mathbf{f}$，然后根据 运动方程 进行 时间积分 得到 广义加速度 $\mathbf{\ddot{q}}$，然后再进行 时间积分 得到 广义速度 $\mathbf{\dot{q}}$，最后进行 时间积分 得到 广义位置 $\mathbf{q}$

时间积分 Time Integration

我们采用 反向欧拉 Backward Euler 作为时间积分方法，这是一种隐式积分方法，它通过在当前积分步骤中使用 “未来的位置” 来保持数值积分的稳定性

我们首先进行一次积分以得到 广义速度

$$\mathbf{M} \mathbf{\dot{q}}^{t+1} = \mathbf{M} \mathbf{\dot{q}}^{t} + \Delta t \ \mathbf{f}(\mathbf{q}^{t+1})$$

• $\mathbf{M}$ 是质量矩阵
• $\mathbf{\dot{q}}^{t+1}$ 是下一步的 广义速度
• $\mathbf{f}(\mathbf{q}^{t+1})$ 是下一步位置的 广义力
• $\Delta t$ 是时间步长

然后再对 广义速度 积分，得到 广义位置

$$\mathbf{q}^{t+1} = \mathbf{q}^{t} + \Delta t \ \mathbf{\dot{q}}^{t+1}$$

• $\mathbf{q}^{t+1}$ 是下一步的 广义位置
• $\mathbf{q}^{t}$ 是当前步的 广义位置
• $\mathbf{\dot{q}}^{t+1}$ 是下一步的 广义速度
• $\Delta t$ 是时间步长

我们发现 广义力 $\mathbf{f}(\mathbf{q}^{t+1})$ 是非线性函数，直接求解比较困难，所以我们需要进行线性化近似，我们使用泰勒展开对其进行一阶近似

$$\mathbf{f}(\mathbf{q}^{t+1}) \approx \mathbf{f}(\mathbf{q}^{t}) + \frac{\partial\mathbf{f}}{\partial\mathbf{q}}\bigg|_{\mathbf{q}^t} (\mathbf{q}^{t+1} - \mathbf{q}^t)$$

• $\mathbf{f}(\mathbf{q}^{t})$ 是当前时刻 $t$ 的 广义力
• $\frac{\partial\mathbf{f}}{\partial\mathbf{q}}\bigg|_{\mathbf{q}^t}$ 是当前时刻 广义力 对 广义坐标 的导数，也就是力梯度

根据之前的 广义位置 积分等式，可得

$$\begin{align*} \mathbf{q}^{t+1} - \mathbf{q}^{t} &= \Delta t \ \mathbf{\dot{q}}^{t+1} \\ \mathbf{f}(\mathbf{q}^{t+1}) &= \mathbf{f}(\mathbf{q}^{t}) + \frac{\partial\mathbf{f}}{\partial\mathbf{q}}\bigg|_{\mathbf{q}^t} \Delta t \ \mathbf{\dot{q}}^{t+1} \end{align*}$$

我们将这个结果代回 广义速度 积分公式

$$\begin{align*} \mathbf{M} \mathbf{\dot{q}}^{t+1} &= \mathbf{M} \mathbf{\dot{q}}^{t} + \Delta t \ \left(\mathbf{f}(\mathbf{q}^{t}) + \frac{\partial\mathbf{f}}{\partial\mathbf{q}}\bigg|_{\mathbf{q}^t} \Delta t \ \mathbf{\dot{q}}^{t+1} \right) \\ &= \mathbf{M} \mathbf{\dot{q}}^{t} + \Delta t \ \mathbf{f}(\mathbf{q}^{t}) + \frac{\partial\mathbf{f}}{\partial\mathbf{q}}\bigg|_{\mathbf{q}^t} \Delta t^2 \ \mathbf{\dot{q}}^{t+1} \end{align*}$$

这里我们特别关注一下 $\frac{\partial\mathbf{f}}{\partial\mathbf{q}}\bigg|_{\mathbf{q}^t}$，它是力梯度，也被称为 刚度矩阵 Stiffness Matrix，记为 $\mathbf{K}$

$$\mathbf{K} = \frac{\partial\mathbf{f}}{\partial\mathbf{q}}\bigg|_{\mathbf{q}^t}$$

同时，由于 Hessian 矩阵 是某个标量函数 (例如我们的势能函数 $V(\mathbf{q})$ 就是一个标量函数) 的二阶导数矩阵，描述了该标量函数的曲率

$$\mathbf{H} = \frac{\partial^2 V}{\partial\mathbf{q}^2}$$

而我们的力函数 $\mathbf{f}(\mathbf{q})$ 是根据势能函数的 负梯度 (负一阶偏导) 定义的，所以它也是 势能的 Hessian (差了个负号)

$$\mathbf{K} = \frac{\partial\mathbf{f}}{\partial\mathbf{q}} = \frac{\partial}{\partial\mathbf{q}}\left(-\frac{\partial V}{\partial\mathbf{q}} \right) = -\frac{\partial^2 V}{\partial\mathbf{q}^2} = -\mathbf{H}$$

如此，我们的积分等式变为

$$\mathbf{M} \mathbf{\dot{q}}^{t+1} = \mathbf{M} \mathbf{\dot{q}}^{t} + \Delta t \ \mathbf{f}(\mathbf{q}^{t}) + \mathbf{K} \Delta t^2 \ \mathbf{\dot{q}}^{t+1}$$

重新整理，我们得到

$$\begin{align*} (\mathbf{M} - \mathbf{K} \Delta t^2) \mathbf{\dot{q}}^{t+1} &= \mathbf{M} \mathbf{\dot{q}}^{t} + \Delta t \ \mathbf{f}(\mathbf{q}^{t}) \\ \mathbf{q}^{t+1} &= \mathbf{q}^t + \Delta t \mathbf{\dot{q}}^{t+1} \end{align*}$$

其中也可以用 $-\mathbf{H}$ 取代 $\mathbf{K}$，具体来说，也就是

$$\begin{align*} \mathbf{f} &= -\frac{\partial}{\partial\mathbf{q}} \sum_{j=0}^{m-1} V_j(\mathbf{E}_j\mathbf{q}) \\ \mathbf{K} &= \frac{\partial\mathbf{f}}{\partial\mathbf{q}} \\ \mathbf{K} &= -\frac{\partial^2}{\partial\mathbf{q}^2} \sum_{j=0}^{m-1} V_j(\mathbf{E}_j\mathbf{q}) \\ &= -\sum_{j=0}^{m-1}\frac{\partial^2}{\partial\mathbf{q}^2} V_j(\mathbf{E}_j\mathbf{q}) \\ &= \sum_{j=0}^{m-1} \left( -\mathbf{E}^T_j \frac{\partial^2 V_j}{\partial\mathbf{q}_j^2} \mathbf{E}_j \right) \\ &= \sum_{j=0}^{m-1} \left( \mathbf{E}^T_j (-\mathbf{H}_j) \mathbf{E}_j \right) \\ &= \sum_{j=0}^{m-1} \left( \mathbf{E}^T_j \mathbf{K}_j \mathbf{E}_j \right) \end{align*}$$

• 引入了逐弹簧的 $-\mathbf{H}_j$ 和 $\mathbf{K}_j$

和之前一样，所有单弹簧的矩阵加起来合成了大的 $\mathbf{K}$ 矩阵