2D弹簧质点系统的隐式求解

2D弹簧质点系统的隐式求解

本文主要记录如何用隐式方法求解弹簧质点系统，也是对Games201课程作业一的总结。

本文主要分为以下几个部分：

• 弹簧质点系统的简单介绍

• 隐式方法的推导

• Jacobi矩阵的推导

• 数值方法的简单介绍

• 程序验证

• 总结

一、弹簧质点系统

弹簧质点系统在柔体、弹性体中有广泛的应用。根据胡克定律对弹簧系统建模：

$${\mathbf{F}}_s=-K(|{\mathbf{x}}_{ij}-r|){\hat{\mathbf{x}}}_{ij}\text{\hspace{1em}(1)}$$

其中${\mathbf{x}}_{ij}={\mathbf{x}}_i-{\mathbf{x}}_j$ 。系统中两个质点受力相反。${\hat{\mathbf{x}}}_{ij}$为单位向量。

二、隐式方法的推导

相对于显式的方法，隐式方法更稳定，允许较大时步。关于显式、隐式方法的概述，可参考这个course note。

下面进行隐式方法的推导：
$${\mathbf{x}}_{t+1}={\mathbf{x}}_t+△t{\mathbf{v}}_{t+1}\text{\hspace{1em}(2)}$$

$${\mathbf{v}}_{t+1}={\mathbf{v}}_t+△t{\mathbf{M}}^{-1}\mathbf{f}({\mathbf{x}}_{t+1})\text{\hspace{1em}(3)}$$

其中：
$$\mathbf{x}={[{\mathbf{x}}_1,...,{\mathbf{x}}_n]}^T$$

$$\mathbf{v}={[{\mathbf{v}}_1,...,{\mathbf{v}}_n]}^T$$

$$\mathbf{f}(\mathbf{x})={[{\mathbf{f}}_1{({\mathbf{x}}_1,...,{\mathbf{x}}_n)}^T,...,{\mathbf{f}}_n{({\mathbf{x}}_1,...,{\mathbf{x}}_n)}^T]}^T$$

将（2）代入（3）得：
$${\mathbf{v}}_{t+1}={\mathbf{v}}_t+△t{\mathbf{M}}^{-1}\mathbf{f}({\mathbf{x}}_t+△t{\mathbf{v}}_{t+1})\text{\hspace{1em}(4)}$$
注意到（4）式右侧$\mathbf{f}$里的参数，（4）式是非线性的，所以对参数进行一阶泰勒展开，得：
$${\mathbf{v}}_{t+1}={\mathbf{v}}_t+△t{\mathbf{M}}^{-1}[f({\mathbf{x}}_t)+\frac{\partial \mathbf{f}}{\partial \mathbf{x}}({\mathbf{x}}_t)△t{\mathbf{v}}_{t+1}]\text{\hspace{1em}(5)}$$
整理一下，得：
$$\mathbf{M}{\mathbf{v}}_{t+1}=\mathbf{M}{\mathbf{v}}_t+△t[f({\mathbf{x}}_t)+\frac{\partial \mathbf{f}}{\partial \mathbf{x}}({\mathbf{x}}_t)△t{\mathbf{v}}_{t+1}]\text{\hspace{1em}(6)}$$
令$\mathbf{J}=\frac{\partial \mathbf{f}}{\partial \mathbf{x}}({\mathbf{x}}_t)$,也就是雅可比矩阵。

整理后，可以得到：
$$[\mathbf{M}-{△t}^2\mathbf{J}]{\mathbf{v}}_{t+1}=\mathbf{M}{\mathbf{v}}_t+△t\mathbf{f}({\mathbf{x}}_t)\text{\hspace{1em}(7)}$$

$$\mathbf{A}{\mathbf{v}}_{t+1}=\mathbf{b}\text{\hspace{1em}(8)}$$

解此方程组，得到结果。

三、Jacobi矩阵的推导

上一部分我们得到了公式（7），解这个方程组之前我们要先求出$\mathbf{J}$。首先罗列一些基本公式：

标量函数对矢量求导：
$$\frac{\partial f}{\partial \mathbf{x}}=[\frac{\partial f}{\partial x_i}\frac{\partial f}{\partial x_j}\frac{\partial f}{\partial x_k}]\text{\hspace{1em}(9)}$$
矢量函数对矢量求导：
$$\frac{\partial \mathbf{f}}{\partial \mathbf{x}}=\left(\begin{array}{ccc}\frac{\partial f_i}{\partial x_i}& \frac{\partial f_i}{\partial x_j}& \frac{\partial f_i}{\partial x_k}\\ \frac{\partial f_j}{\partial x_i}& \frac{\partial f_j}{\partial x_j}& \frac{\partial f_j}{\partial x_k}\\ \frac{\partial f_k}{\partial x_i}& \frac{\partial f_k}{\partial x_j}& \frac{\partial f_k}{\partial x_k}\end{array}\right)\text{\hspace{1em}(10)}$$

单位向量对向量求导：
$$\frac{\partial \hat{\mathbf{x}}}{\partial \mathbf{x}}=\frac{\mathbf{I}-\hat{\mathbf{x}}\cdot {\hat{\mathbf{x}}}^T}{|\mathbf{x}|}\text{\hspace{1em}(11)}$$
其中，$\mathbf{I}$是单位矩阵。

对于弹力，从胡克定律着手进行雅可比矩阵推导；对于弹簧衰减力、重力等与位置无关的力，这里使用显式的方法去计算。

下面从式（1）入手，推导弹力的雅可比矩阵：

这里只考虑二维情况，首先考虑一个系统里只有两个质点$i$和$j$，那么按照式(10)的定义，我们求导后会分别得到${\mathbf{J}}_{ii}$，${\mathbf{J}}_{ij}$，${\mathbf{J}}_{ji}$，${\mathbf{J}}_{jj}$四个小矩阵，以计算${\mathbf{J}}_{ii}$为例进行推导：
$$\frac{\partial {\mathbf{F}}_s}{\partial {\mathbf{x}}_i}={\mathbf{J}}_{ii}=-k[({\mathbf{x}}_{ij}-r)\frac{\partial {\hat{\mathbf{x}}}_{ij}}{{\mathbf{x}}_i}+{\hat{\mathbf{x}}}_{ij}\frac{\partial |{\mathbf{x}}_{ij}|-r}{{\mathbf{x}}_i}]\text{\hspace{1em}(12)}$$
${\mathbf{F}}_s$为质点$i$所受弹力。

利用公式（11），得：
$$\frac{\partial {\mathbf{F}}_s}{\partial {\mathbf{x}}_i}={\mathbf{J}}_{ii}=-k[(|{\mathbf{x}}_{ij}-r|)(\frac{\mathbf{I}-{\hat{\mathbf{x}}}_{ij}\cdot {\hat{\mathbf{x}}}_{ij}^T}{|{\mathbf{x}}_{ij}|})+{\hat{\mathbf{x}}}_{ij}\cdot {\hat{\mathbf{x}}}_{ij}^T]\text{\hspace{1em}(13)}$$
变一下形：
$$\frac{\partial {\mathbf{F}}_s}{\partial {\mathbf{x}}_i}={\mathbf{J}}_{ii}=-k[(1-\frac{r}{|{\mathbf{x}}_{ij}|})(\mathbf{I}-{\hat{\mathbf{x}}}_{ij}\cdot {\hat{\mathbf{x}}}_{ij}^T)+{\hat{\mathbf{x}}}_{ij}\cdot {\hat{\mathbf{x}}}_{ij}^T]\text{\hspace{1em}(14)}$$
将右侧圆括号展开后变形得：
$$\frac{\partial {\mathbf{F}}_s}{\partial {\mathbf{x}}_i}={\mathbf{J}}_{ii}=-k[\mathbf{I}-\frac{r}{|{\mathbf{x}}_{ij}|}(\mathbf{I}-{\hat{\mathbf{x}}}_{ij}\cdot {\hat{\mathbf{x}}}_{ij}^T)]\text{\hspace{1em}(15)}$$
同理，可以推出其余三个小矩阵：
$${\mathbf{J}}_{ii}=-{\mathbf{J}}_{ij}={\mathbf{J}}_{jj}-{\mathbf{J}}_{ji}\text{\hspace{1em}(16)}$$
至此，雅可比矩阵推导完毕，我们需要在每一时步都组装这个矩阵，得到：

$$\frac{\partial \mathbf{f}}{\partial \mathbf{x}}=\left(\begin{array}{cc}{\mathbf{J}}_{ii}& {\mathbf{J}}_{ij}\\ {\mathbf{J}}_{ji}& {\mathbf{J}}_{jj}\end{array}\right)\text{\hspace{1em}(10)}$$

四、数值方法的简单介绍

至此，我们就可以着手去解这个线性方程组了。本人在数值计算方法方面的学习一拖再拖，还没学多少，只知道可以用Jacobi method、Gauss–Seidel Method、Gradient Descent Method等。关于这些算法网上有大量资料，这里就懒得写了。

五、程序验证

用pixi.js做弹簧的可视化，实现了弹簧质点系统的显式方法，代码地址点这里

在线demo点这里

本文推导的隐式积分法的程序正在coding中…

六、总结

本文主要参考Games201第二讲、Miles Macklin博客里的文章《implict springs》、Matthias Muller等人的《Real Time Physics Class Notes》,其中貌似发现《Real Time Physics Class Notes》里 的一个公式推导的错误，第17页的公式3.32多除了一个$l^2$。