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

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

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

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

弹簧质点系统的简单介绍
隐式方法的推导
Jacobi矩阵的推导
数值方法的简单介绍
程序验证
总结

一、弹簧质点系统

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

Fs=−K(∣xij−r∣)x^ij(1)\boldsymbol{F}_s=-K(|\boldsymbol{x}_{ij}-r|)\hat{\boldsymbol{x}}_{ij}\tag{1} Fs=−K(∣xij−r∣)x^ij(1)

其中xij=xi−xj\boldsymbol{x}_{ij}=\boldsymbol{x}_{i}-\boldsymbol{x}_{j}xij=xi−xj 。系统中两个质点受力相反。x^ij\hat{\boldsymbol{x}}_{ij}x^ij为单位向量。

二、隐式方法的推导

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

下面进行隐式方法的推导：
xt+1=xt+△tvt+1(2)\boldsymbol{x}_{t+1}=\boldsymbol{x}_t+\triangle{t}\boldsymbol{v}_{t+1}\tag{2} xt+1=xt+△tvt+1(2)

vt+1=vt+△tM−1f(xt+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=[x1,...,xn]T\boldsymbol{x}={[\boldsymbol{x}_{1},...,\boldsymbol{x}_{n}]}^T x=[x1,...,xn]T

v=[v1,...,vn]T\boldsymbol{v}={[\boldsymbol{v}_{1},...,\boldsymbol{v}_{n}]}^T v=[v1,...,vn]T

f(x)=[f1(x1,...,xn)T,...,fn(x1,...,xn)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）得：
vt+1=vt+△tM−1f(xt+△tvt+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）式是非线性的，所以对参数进行一阶泰勒展开，得：
vt+1=vt+△tM−1[f(xt)+∂f∂x(xt)△tvt+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)
整理一下，得：
Mvt+1=Mvt+△t[f(xt)+∂f∂x(xt)△tvt+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(xt)\boldsymbol{J}=\frac{\partial{\boldsymbol{f}}}{\partial{\boldsymbol{x}}}(\boldsymbol{x}_t)J=∂x∂f(xt),也就是雅可比矩阵。

整理后，可以得到：
[M−△t2J]vt+1=Mvt+△tf(xt)(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)

Avt+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∂xi∂f∂xj∂f∂xk](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=(∂fi∂xi∂fi∂xj∂fi∂xk∂fj∂xi∂fj∂xj∂fj∂xk∂fk∂xi∂fk∂xj∂fk∂xk)(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）入手，推导弹力的雅可比矩阵：

这里只考虑二维情况，首先考虑一个系统里只有两个质点iii和jjj，那么按照式(10)的定义，我们求导后会分别得到Jii\boldsymbol{J}_{ii}Jii，Jij\boldsymbol{J}_{ij}Jij，Jji\boldsymbol{J}_{ji}Jji，Jjj\boldsymbol{J}_{jj}Jjj四个小矩阵，以计算Jii\boldsymbol{J}_{ii}Jii为例进行推导：
∂Fs∂xi=Jii=−k[(xij−r)∂x^ijxi+x^ij∂∣xij∣−rxi](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)
Fs\boldsymbol{F}_sFs为质点iii所受弹力。

利用公式（11），得：
∂Fs∂xi=Jii=−k[(∣xij−r∣)(I−x^ij⋅x^ijT∣xij∣)+x^ij⋅x^ijT](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)
变一下形：
∂Fs∂xi=Jii=−k[(1−r∣xij∣)(I−x^ij⋅x^ijT)+x^ij⋅x^ijT](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)
将右侧圆括号展开后变形得：
∂Fs∂xi=Jii=−k[I−r∣xij∣(I−x^ij⋅x^ijT)](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)
同理，可以推出其余三个小矩阵：
Jii=−Jij=Jjj−Jji(16)\boldsymbol{J}_{ii}=-\boldsymbol{J}_{ij}=\boldsymbol{J}_{jj}-\boldsymbol{J}_{ji}\tag{16} Jii=−Jij=Jjj−Jji(16)
至此，雅可比矩阵推导完毕，我们需要在每一时步都组装这个矩阵，得到：

∂f∂x=(JiiJijJjiJjj)(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做弹簧的可视化，实现了弹簧质点系统的显式方法，代码地址点这里

在线demo点这里

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

六、总结

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

2D弹簧质点系统的隐式求解相关推荐

网格弹簧质点系统模拟（Spring-Mass System by Verlet Integration）附源码
模拟物体变形最简单的方法就是采用弹簧质点系统(Spring-Mass System),由于模型简单并且实用,它已被广泛应用于服饰.毛发以及弹性固体的动态模拟.对于三角网格而言,弹簧质点系统将网格中的顶 ...
多级弹簧-质量系统瞬态分析（基于Newmark）
该程序主要用于实现多自由度下动态系统的隐式分析,计算模型参见多级弹簧-质量系统的瞬态计算(基于中心差分法).下面是具体的Python程序. #coding:utf-8 ""&quo ...
面向程序员的数据挖掘指南-----第三章：隐式评价和基于物品的过滤算法
本章会从用户的评价类型开始讨论,包括显式评价(赞一下.踩一脚.五星评价等等)和隐式评价(比如在亚马逊上购买了MP3,我们可以认为他喜欢这个产品). 内容: 显式评价隐式评价哪种评价方式更准确? 基 ...
【JY】结构动力学之显隐式
因你精彩即刻关注 "通俗讲动力数值算法" 随着科学研究和工程技术的不断发展,出现了诸如航天飞机.空间站.海洋石油钻井平台等大型或超大型的复杂结构,它们不但自由度高,而且还含有非线 ...
图形学笔记（二十）粒子、刚体、流体的模拟—— 欧拉方法、Errors 和 Instability、中点法、自适应步长、隐式欧拉方法、Runge-Kutta方法、刚体与流体模拟（质点法、网格法、MPM）
图形学笔记(十九)粒子.刚体.流体的模拟-- 欧拉方法.Errors 和 Instability.中点法.自适应步长.隐式欧拉方法.Runge-Kutta方法.刚体与流体模拟(质点法.网格法.MPM) ...
matlab将求解sin隐式解,Matlab隐式符号方程求解和赋值
近日处理了一个隐式方程的求解,由于方程含有较多的未知数,而且这些参数均是跟实验相关的一些参数,所以,必须得到需要求解的解与这些参数之间的一个表达式.之前是考虑用的Maple推导求解了该隐私方程,求解 ...
使用隐式Intent打开系统浏览器的百度网页
使用隐式Intent,我们不仅可以启动自己程序内的活动,还可以启动其它程序的活动,这使得Android多个应用程序之间的功能共享成为了可能.比如说你的应用程序中需要展示一个网页,这时你没有必要自己去实 ...
Android 通过 “隐式意图” 打开系统的浏览器访问百度页面
在MainActivity中,通过"隐式意图"打开系统的浏览器访问百度页面: MainActivity页面: package cn.lwx.openbrowser;import a ...

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

2D弹簧质点系统的隐式求解相关推荐

最新文章

热门文章