参考资料：

Fast Simulation of Mass-Spring Systems in Rust

论文阅读：Fast Simulation of Mass-Spring Systems

【论文精读】讲解刘天添2013年的fast simulation of mass spring system(Projective Dynamics最早的论文)

Projective Dynamics笔记(一)——Fast Mass Spring

概述

这篇论文通过引入弹簧方向的辅助变量，并采用块坐标下降法解决了传统隐式欧拉法中速度慢的问题，成功实现了弹簧质点系统的快速仿真。该方法特别适用于实时应用，并且在低迭代次数下也能得到较好的视觉效果。

流程概述：

1.前置知识

1.1 运动方程（牛顿第二定律）

在物理仿真中，系统的运动通常由二阶微分方程描述，例如：$M\frac{d^{2}q}{d{t}^2}=f(q)$
其中：

• ( M ) 是质量矩阵，表示系统中各质点的质量。
• ( q(t) ) 是位置向量，表示质点的位置。
• ( f(q) ) 是作用在质点上的力，通常是位移 ( q ) 的函数。

1.2 二阶导数的离散化

使用中心差分法，二阶导数 ( $\frac{d^{2}q}{d{t}^2}$ ) 的离散形式为：
$$\frac{d^{2}q}{d{t}^2}\approx \frac{q_{n+1}-2q_n+q_{n-1}}{h^2}$$
其中：

• ( q_n ) 是时间 ( t_n ) 时的位置。
• ( q_{n+1} ) 是时间 ( t_{n+1} = t_n + h ) 时的位置。
• ( q_{n-1} ) 是时间 ( t_{n-1} = t_n - h ) 时的位置。
• ( h ) 是时间步长

1.3 代入运动方程

将二阶导数的离散化形式代入运动方程 ( $M\frac{d^{2}q}{d{t}^2}=f(q)$ )：
$$M\frac{q_{n+1}-2q_n+q_{n-1}}{h^2}=f(q_{n+1})$$
最后我们得到隐式欧拉法的公式：
$$M(q_{n+1}-2q_n+q_{n-1})=h^{2}f(q_{n+1})$$
其中，( $q_{n+1}$ ) 是需要通过求解获得的未知位置，( $f(q_{n+1})$ ) 依赖于 ( $q_{n+1}$ )，因此这是一个隐式公式

1.4 物理意义

2. 将隐式积分问题转化为一个优化问题

2.1 要解的是隐式积分问题是：

$$M\frac{q_{n+1}-2q_n+q_{n-1}}{h^2}=f(q_{n+1})$$
其中：
n是时间迭代步
q是所有粒子的位置向量
M是粒子质量对角矩阵
h是时间步长
f是整个系统的保守力
$q_{n+1}$是未知状态量，设为x。$q_n$和$q_{n-1}$是已知量，设为$y=2q_n-q_{n-1}$
所以得到式子$M(x-y)=h^{2}f(x)$
求解这个方程，等效于求解下面这个方程的极小值:
(令g(x)求导为0得到上式)
$$g(x)=\frac{1}{2}(x-y{)}^{T}M(x-y)+h^{2}E(x)$$

其中，( E ) 为系统的势能（因为 ( $\nabla E=-f$ )，因此 ( $\nabla g=0$ ) 等效于公式 (7)）。

按照胡克定律，弹簧的弹性势能为：

$$E=\frac{1}{2}k(\Vert p_1-p_2\Vert -r{)}^2$$

其中：

• ( $p_1,p_2$ ) 为两个粒子的位置，
• ( r ) 为弹簧的静止长度（rest length）。

但如果直接采用这个形式，上面的 ( g(x) ) 极值问题就不太好解。为了将上式变形为一个方便求解的形式，作者引入了一个辅助变量 ( $d\in {\mathbb{R}}^3$ )。

2.2 引入辅助变量d

公式如下：
假设d是一个未知的三位向量，那么：
$$(\Vert p_1-p_2\Vert -r{)}^2=\underset{\Vert d\Vert =r}{\min }\Vert p_1-p_2-d{\Vert }^2$$

1. 左边公式的物理意义：

左边的公式 ( $(\Vert p_1-p_2\Vert -r{)}^2$ ) 是弹簧势能的表示形式，依据胡克定律，它表示弹簧当前长度 ( |p_1 - p_2| ) 与静止长度 ( r ) 之间的差的平方。这里 ( $p_1$ ) 和 ( $p_2$ ) 是弹簧两端质点的位置

2. 右边公式的几何解释：

右边的公式引入了一个辅助向量 ( d )，并对其施加约束 ( |d| = r )。这个向量表示固定长度为 ( r ) 的向量，但方向可以自由变化。优化问题为：
$$\underset{\Vert d\Vert =r}{\min }\Vert p_1-p_2-d{\Vert }^2$$
这个问题的意思是寻找一个向量 ( d )，使得 ( $p_1-p_2$ ) 与 ( d ) 的差最小，即让 ( $\Vert p_1-p_2-d\Vert$ ) 最小化

$$显然当d=r\frac{p_1-p_2}{\Vert p_1-p_2\Vert }时取极小值$$

重新定义弹簧的弹性势能E(x,d)

矩阵L和J的推导

第一行的式子为什么可以等于L和J，证明如下：
忽略$d_i^T{d}_i$是关于 d 的平方项，但这项对 x 的优化没有直接影响，因此我们暂时忽略它，只关注 x 和 d 之间的关系

$S_i^T$是一个选择矩阵，是一个单位向量(标准基），用于选择第 𝑖个弹簧的位移变量 $d_i$
​
在这里，$d=S_i^{T}d$ 可以理解为，向量 $d$ 中的每个分量 $d_i$ 对应一个特定的弹簧的偏移量。通过 $S_i^{T}d$，我们提取了 $d$ 中与第 $i$ 个弹簧相关的那部分位移。

换句话说，$S_i^T$ 确保我们从总的 $d$ 向量中只选取第 $i$ 个元素（因为 $S_i^T$ 是一个标准基，其他位置上的元素都会被置为 0）。

因此，对于每个弹簧 $i$，我们有：

$$d_i=S_i^{T}d$$

这表示 $d_i$ 是通过选择矩阵 $S_i^T$ 从 $d$ 中提取出来的。

外力为什么放在$X^T（）$里面

外力 $f_{ext}$被视作对系统的额外作用力，所以在总的能量函数中，它会影响线性部分，即外力对位移 x 的作用是线性的。因此，外力项自然可以放入这个线性项中
乘以 $h^2$体现了外力在整个时间步长内的累积效应
力和加速度是直接相关的。加速度是位移 𝑥对时间 𝑡的二阶导数。而当我们离散化这个二阶导数时，时间步长 ℎ被平方了

2.3 Block Coordinate Descent（块坐标下降法）

Global Step部分