2D弹簧质点系统的隐式求解
本文主要记录如何用隐式方法求解弹簧质点系统,也是对Games201课程作业一的总结。
本文主要分为以下几个部分:
-
弹簧质点系统的简单介绍
-
隐式方法的推导
-
Jacobi矩阵的推导
-
数值方法的简单介绍
-
程序验证
-
总结
一、弹簧质点系统
弹簧质点系统在柔体、弹性体中有广泛的应用。根据胡克定律对弹簧系统建模:
F s = − K ( ∣ x i j − r ∣ ) x ^ i j (1) \boldsymbol{F}_s=-K(|\boldsymbol{x}_{ij}-r|)\hat{\boldsymbol{x}}_{ij}\tag{1} Fs=−K(∣xij−r∣)x^ij(1)
其中 x i j = x i − x j \boldsymbol{x}_{ij}=\boldsymbol{x}_{i}-\boldsymbol{x}_{j} xij=xi−xj 。系统中两个质点受力相反。 x ^ i j \hat{\boldsymbol{x}}_{ij} x^ij为单位向量。
二、隐式方法的推导
相对于显式的方法,隐式方法更稳定,允许较大时步。关于显式、隐式方法的概述,可参考这个course note。
下面进行隐式方法的推导:
x
t
+
1
=
x
t
+
△
t
v
t
+
1
(2)
\boldsymbol{x}_{t+1}=\boldsymbol{x}_t+\triangle{t}\boldsymbol{v}_{t+1}\tag{2}
xt+1=xt+△tvt+1(2)
v t + 1 = v t + △ t M − 1 f ( x t + 1 ) (3) \boldsymbol{v}_{t+1}=\boldsymbol{v}_t+\triangle{t}\boldsymbol{M}^{-1}\boldsymbol{f}(\boldsymbol{x}_{t+1})\tag{3} vt+1=vt+△tM−1f(xt+1)(3)
其中:
x
=
[
x
1
,
.
.
.
,
x
n
]
T
\boldsymbol{x}={[\boldsymbol{x}_{1},...,\boldsymbol{x}_{n}]}^T
x=[x1,...,xn]T
v = [ v 1 , . . . , v n ] T \boldsymbol{v}={[\boldsymbol{v}_{1},...,\boldsymbol{v}_{n}]}^T v=[v1,...,vn]T
f ( x ) = [ f 1 ( x 1 , . . . , x n ) T , . . . , f n ( x 1 , . . . , x n ) T ] T \boldsymbol{f(\boldsymbol{x})}={[\boldsymbol{f}_1{(\boldsymbol{x}_1,...,\boldsymbol{x}_n)}^T,...,\boldsymbol{f}_n{(\boldsymbol{x}_1,...,\boldsymbol{x}_n)}^T]}^T f(x)=[f1(x1,...,xn)T,...,fn(x1,...,xn)T]T
将(2)代入(3)得:
v
t
+
1
=
v
t
+
△
t
M
−
1
f
(
x
t
+
△
t
v
t
+
1
)
(4)
\boldsymbol{v}_{t+1}=\boldsymbol{v}_t+\triangle{t}\boldsymbol{M}^{-1}\boldsymbol{f}(\boldsymbol{x}_{t}+\triangle{t}\boldsymbol{v}_{t+1})\tag{4}
vt+1=vt+△tM−1f(xt+△tvt+1)(4)
注意到(4)式右侧
f
\boldsymbol{f}
f里的参数,(4)式是非线性的,所以对参数进行一阶泰勒展开,得:
v
t
+
1
=
v
t
+
△
t
M
−
1
[
f
(
x
t
)
+
∂
f
∂
x
(
x
t
)
△
t
v
t
+
1
]
(5)
\boldsymbol{v}_{t+1}=\boldsymbol{v}_t+\triangle{t}\boldsymbol{M}^{-1}[f(\boldsymbol{x}_{t})+\frac{\partial\boldsymbol{f}}{\partial\boldsymbol{x}}(\boldsymbol{x}_t)\triangle{t}\boldsymbol{v}_{t+1}]\tag{5}
vt+1=vt+△tM−1[f(xt)+∂x∂f(xt)△tvt+1](5)
整理一下,得:
M
v
t
+
1
=
M
v
t
+
△
t
[
f
(
x
t
)
+
∂
f
∂
x
(
x
t
)
△
t
v
t
+
1
]
(6)
\boldsymbol{M}\boldsymbol{v}_{t+1}=\boldsymbol{M}\boldsymbol{v}_t+\triangle{t}[f(\boldsymbol{x}_{t})+\frac{\partial\boldsymbol{f}}{\partial\boldsymbol{x}}(\boldsymbol{x}_t)\triangle{t}\boldsymbol{v}_{t+1}]\tag{6}
Mvt+1=Mvt+△t[f(xt)+∂x∂f(xt)△tvt+1](6)
令
J
=
∂
f
∂
x
(
x
t
)
\boldsymbol{J}=\frac{\partial{\boldsymbol{f}}}{\partial{\boldsymbol{x}}}(\boldsymbol{x}_t)
J=∂x∂f(xt),也就是雅可比矩阵。
整理后,可以得到:
[
M
−
△
t
2
J
]
v
t
+
1
=
M
v
t
+
△
t
f
(
x
t
)
(7)
[\boldsymbol{M}-{\triangle{t}}^2\boldsymbol{J}]\boldsymbol{v}_{t+1}=\boldsymbol{M}\boldsymbol{v}_{t}+\triangle{t}\boldsymbol{f}(\boldsymbol{x}_t)\tag{7}
[M−△t2J]vt+1=Mvt+△tf(xt)(7)
A v t + 1 = b (8) \boldsymbol{A}\boldsymbol{v}_{t+1}=\boldsymbol{b}\tag{8} Avt+1=b(8)
解此方程组,得到结果。
三、Jacobi矩阵的推导
上一部分我们得到了公式(7),解这个方程组之前我们要先求出 J \boldsymbol{J} J。首先罗列一些基本公式:
标量函数对矢量求导:
∂
f
∂
x
=
[
∂
f
∂
x
i
∂
f
∂
x
j
∂
f
∂
x
k
]
(9)
\frac{\partial f}{\partial\boldsymbol{x}} = [\frac{\partial f}{\partial x_i}\ \frac{\partial f} {\partial x_j}\ \frac{\partial f} {\partial x_k}]\tag{9}
∂x∂f=[∂xi∂f ∂xj∂f ∂xk∂f](9)
矢量函数对矢量求导:
∂
f
∂
x
=
(
∂
f
i
∂
x
i
∂
f
i
∂
x
j
∂
f
i
∂
x
k
∂
f
j
∂
x
i
∂
f
j
∂
x
j
∂
f
j
∂
x
k
∂
f
k
∂
x
i
∂
f
k
∂
x
j
∂
f
k
∂
x
k
)
(10)
\frac{\partial\boldsymbol{f}}{\partial\boldsymbol{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)\tag{10}
∂x∂f=⎝⎜⎛∂xi∂fi∂xi∂fj∂xi∂fk∂xj∂fi∂xj∂fj∂xj∂fk∂xk∂fi∂xk∂fj∂xk∂fk⎠⎟⎞(10)
单位向量对向量求导:
∂
x
^
∂
x
=
I
−
x
^
⋅
x
^
T
∣
x
∣
(11)
\frac{\partial \hat{\boldsymbol{x}}}{\partial \boldsymbol{x} } = \frac{\boldsymbol{I}-\hat{\boldsymbol{x}}\cdot{\hat{\boldsymbol{x}}^T}}{|\boldsymbol{x}|}\tag{11}
∂x∂x^=∣x∣I−x^⋅x^T(11)
其中,
I
\boldsymbol{I}
I是单位矩阵。
对于弹力,从胡克定律着手进行雅可比矩阵推导;对于弹簧衰减力、重力等与位置无关的力,这里使用显式的方法去计算。
下面从式(1)入手,推导弹力的雅可比矩阵:
这里只考虑二维情况,首先考虑一个系统里只有两个质点
i
i
i和
j
j
j,那么按照式(10)的定义,我们求导后会分别得到
J
i
i
\boldsymbol{J}_{ii}
Jii,
J
i
j
\boldsymbol{J}_{ij}
Jij,
J
j
i
\boldsymbol{J}_{ji}
Jji,
J
j
j
\boldsymbol{J}_{jj}
Jjj四个小矩阵,以计算
J
i
i
\boldsymbol{J}_{ii}
Jii为例进行推导:
∂
F
s
∂
x
i
=
J
i
i
=
−
k
[
(
x
i
j
−
r
)
∂
x
^
i
j
x
i
+
x
^
i
j
∂
∣
x
i
j
∣
−
r
x
i
]
(12)
\frac{\partial\boldsymbol{F}_{s}}{\partial\boldsymbol{x}_{i}}=\boldsymbol{J}_{ii}=-k[(\boldsymbol{x}_{ij}-r)\frac{\partial\hat{\boldsymbol{x}}_{ij}}{\boldsymbol{x}_i}+\hat{\boldsymbol{x}}_{ij}\frac{\partial|\boldsymbol{x}_{ij}|-r}{\boldsymbol{x}_i}]\tag{12}
∂xi∂Fs=Jii=−k[(xij−r)xi∂x^ij+x^ijxi∂∣xij∣−r](12)
F
s
\boldsymbol{F}_s
Fs为质点
i
i
i所受弹力。
利用公式(11),得:
∂
F
s
∂
x
i
=
J
i
i
=
−
k
[
(
∣
x
i
j
−
r
∣
)
(
I
−
x
^
i
j
⋅
x
^
i
j
T
∣
x
i
j
∣
)
+
x
^
i
j
⋅
x
^
i
j
T
]
(13)
\frac{\partial\boldsymbol{F}_{s}}{\partial\boldsymbol{x}_{i}}=\boldsymbol{J}_{ii}=-k[(|\boldsymbol{x}_{ij}-r|)(\frac{\boldsymbol{I}-\hat{\boldsymbol{x}}_{ij}\cdot{\hat{\boldsymbol{x}}_{ij}^T}}{|\boldsymbol{x}_{ij}|})+\hat{\boldsymbol{x}}_{ij}\cdot{\hat{\boldsymbol{x}}_{ij}^T}]\tag{13}
∂xi∂Fs=Jii=−k[(∣xij−r∣)(∣xij∣I−x^ij⋅x^ijT)+x^ij⋅x^ijT](13)
变一下形:
∂
F
s
∂
x
i
=
J
i
i
=
−
k
[
(
1
−
r
∣
x
i
j
∣
)
(
I
−
x
^
i
j
⋅
x
^
i
j
T
)
+
x
^
i
j
⋅
x
^
i
j
T
]
(14)
\frac{\partial\boldsymbol{F}_{s}}{\partial\boldsymbol{x}_{i}}=\boldsymbol{J}_{ii}=-k[(1-\frac{r}{|\boldsymbol{x}_{ij}|})({\boldsymbol{I}-\hat{\boldsymbol{x}}_{ij}\cdot{\hat{\boldsymbol{x}}_{ij}^T}})+\hat{\boldsymbol{x}}_{ij}\cdot{\hat{\boldsymbol{x}}_{ij}^T}]\tag{14}
∂xi∂Fs=Jii=−k[(1−∣xij∣r)(I−x^ij⋅x^ijT)+x^ij⋅x^ijT](14)
将右侧圆括号展开后变形得:
∂
F
s
∂
x
i
=
J
i
i
=
−
k
[
I
−
r
∣
x
i
j
∣
(
I
−
x
^
i
j
⋅
x
^
i
j
T
)
]
(15)
\frac{\partial\boldsymbol{F}_{s}}{\partial\boldsymbol{x}_{i}}=\boldsymbol{J}_{ii}=-k[\boldsymbol{I}-\frac{r}{|\boldsymbol{x}_{ij}|}(\boldsymbol{I}-\hat{\boldsymbol{x}}_{ij}\cdot{\hat{\boldsymbol{x}}_{ij}^T})]\tag{15}
∂xi∂Fs=Jii=−k[I−∣xij∣r(I−x^ij⋅x^ijT)](15)
同理,可以推出其余三个小矩阵:
J
i
i
=
−
J
i
j
=
J
j
j
−
J
j
i
(16)
\boldsymbol{J}_{ii}=-\boldsymbol{J}_{ij}=\boldsymbol{J}_{jj}-\boldsymbol{J}_{ji}\tag{16}
Jii=−Jij=Jjj−Jji(16)
至此,雅可比矩阵推导完毕,我们需要在每一时步都组装这个矩阵,得到:
∂ f ∂ x = ( J i i J i j J j i J j j ) (10) \frac{\partial\boldsymbol{f}}{\partial\boldsymbol{x}}=\left( \begin{array}{ccc} \boldsymbol{J}_{ii} &\boldsymbol{J}_{ij}\\ \boldsymbol{J}_{ji} & \boldsymbol{J}_{jj}\\ \end{array} \right)\tag{10} ∂x∂f=(JiiJjiJijJjj)(10)
四、数值方法的简单介绍
至此,我们就可以着手去解这个线性方程组了。本人在数值计算方法方面的学习一拖再拖,还没学多少,只知道可以用Jacobi method、Gauss–Seidel Method、Gradient Descent Method等。关于这些算法网上有大量资料,这里就懒得写了。
五、程序验证
用pixi.js做弹簧的可视化,实现了弹簧质点系统的显式方法,代码地址点这里
本文推导的隐式积分法的程序正在coding中…
六、总结
本文主要参考Games201第二讲、Miles Macklin博客里的文章《implict springs》、Matthias Muller等人的《Real Time Physics Class Notes》,其中貌似发现《Real Time Physics Class Notes》里 的一个公式推导的错误,第17页的公式3.32多除了一个 l 2 l^2 l2。