This is a quick post to document some work I did while writing a mass spring simulation using an implicit integrator. Implicit, or backward Euler integration is well described in David Baraff's Physically Based Modelling SIGGRAPH course and this post assumes some familiarity with it.

Springs are a workhorse in physical simulation, once you have unconditionally stable springs you can use them to model just about anything, from rigid bodies to drool and snot. For example, Industrial Light & Magic used a tetrahedral mesh with edge and altitude springs to model the damage to ships in Avatar (see Avatar: Bending Rigid Bodies).

If you sit down and try and implement an implicit integrator one of the first things you need is the Jacobian of the particle forces with respect to the particle positions and velocities. The rest of this post shows how to derive these Jacobians for a basic Hookean spring in a form ready to be plugged into a linear system solver (I use a hand-rolled conjugate gradient solver, see Jonathon Shewchuk's painless introduction for details, it is all of about 20 lines of code to implement).

\[\renewcommand{\v}[1]{\mathbf{#1}} \newcommand{\uv}[1]{\mathbf{\widehat{#1}}} \newcommand\ddx[1]{\frac{\partial#1}{\partial \v{x} }} \newcommand\dd[2]{\frac{\partial#1}{\partial #2}}\]

Vector Calculus Basics

In order to calculate the force Jacobians we first need to know how to calculate the derivatives of some basic geometric quantities with respect to a vector.

In general the derivative of a scalar valued function with respect to a vector is defined as the following row vector of partial derivatives:

\[ \ddx{f} = \begin{bmatrix} \dd{f}{x_i} & \dd{f}{x_j} & \dd{f}{x_k} \end{bmatrix}\]

And for a vector valued function with respect to a vector:

\[\ddx{\v{f}} = \begin{bmatrix} \dd{f_i}{x_i} & \dd{f_i}{x_j} & \dd{f_i}{x_k} \\ \dd{f_j}{x_i} & \dd{f_j}{x_j} & \dd{f_j}{x_k} \\ \dd{f_k}{x_i} & \dd{f_k}{x_j} & \dd{f_k}{x_k} \end{bmatrix}\]

Applying the first definition to the dot product of two vectors we can calculate the derivative with respect to one of the vectors:

\[\ddx{\v{x}^T \cdot \v{y}} = \v{y}^T \]

Note that I'll explicitly keep track of whether vectors are row or column vectors as it will help keep things straight later on.

The derivative of a vector magnitude with respect to the vector, gives the normalized vector transposed:

\[\ddx{|\v{x}|} = \left(\frac{\v{x}}{|\v{x}|}\right)^T = \uv{x}^T \]

The derivative of a normalized vector \(\v{\widehat{x}} = \frac{\v{x}}{|\v{x}|} \) can be obtained using the quotient rule:

\[\ddx{\uv{x}} = \frac{\v{I}|\v{x}| - \v{x}\cdot\uv{x}^T}{|\v{x}|^2}\]

Where \(\v{I}\) is the \(n\) x \(n\) identity matrix and n is the dimension of \(x\). The product of a column vector and a row vector \(\uv{x}\cdot\uv{x}^T\) is the outer product which is a \(n\) x \(n\) matrix that can be constructed using standard matrix multiplication definition.

Dividing through by \(|\v{x}|\) we have:

\[\ddx{\uv{x}} = \frac{\v{I} - \uv{x}\cdot\uv{x}^T}{\v{|x|}}\]

Jacobian of Stretch Force

Now we are ready to compute the Jacobian of the spring forces. Recall the equation for the elastic force on a particle \(i\) due to an undamped Hookean spring:

\[\v{F_s} = -k_s(|\v{x}_{ij}| - r)\uv{x}_{ij}\]

Where \(\v{x}_{ij} = \v{x}_i - \v{x}_j\) is the vector between the two connected particle positions, \(r\) is the rest length and \(k_s\) is the stiffness coefficient.

The Jacobian of this force with respect to particle \(i\)'s position is obtained by using the product rule for the two \(\v{x}_i\) dependent terms in \(\v{F_s}\):

\[\dd{\v{F_s}}{\v{x}_i} = -ks\left[(|\v{x}_{ij}| - r)\dd{\uv{x}_{ij}}{\v{x}_i} + \uv{x}_{ij}\dd{(|\v{x}_{ij}| - r)}{\v{x}_i}\right]\]

Using the previously derived formulas for the derivative of a vector magnitude and normalized vector we have:

\[\dd{\v{F_s}}{\v{x}_i} = -ks\left[(|\v{x}_{ij}| - r)\left(\frac{\v{I} - \uv{x}_{ij}\cdot \uv{x}_{ij}^T}{|\v{x}_{ij}|}\right) + \uv{x}_{ij}\cdot\uv{x}_{ij}^T\right]\]

Dividing the first two terms through by \(|\v{x}_{ij}|\):

\[\dd{\v{F_s}}{\v{x}_i} = -ks\left[(1 - \frac{r}{|\v{x}_{ij}|})\left(\v{I} - \uv{x}_{ij}\cdot \uv{x}_{ij}^T\right) + \uv{x}_{ij}\cdot \uv{x}_{ij}^T\right]\]

Due to the symmetry in the definition of \(\v{x}_{ij}\) we have the following force derivative with respect to the opposite particle:

\[\dd{\v{F_s}}{\v{x}_j} = -\dd{\v{F_s}}{\v{x}_i}\]

Jacobian of Damping Force

The equation for the damping force on a particle \(i\) due to a spring is:

\[\v{F_d} = -k_d\cdot\uv{x}(\v{v}_{ij}\cdot \uv{x}_{ij})\]

Where \(\v{v}_{ij} = \v{v}_i-\v{v}_j\) is the relative velocities of the two particles. This is the preferred formulation because it damps only relative velocity along the spring axis.

Taking the derivative with respect to \(\v{v}_i\):

\[\dd{\v{F_d}}{\v{v}_i} = -k_d\cdot\uv{x}\cdot\uv{x}^T\]

As with stretching, the force on the opposite particle is simply negated:

\[\dd{\v{F_d}}{\v{v}_j} = -\dd{\v{F_d}}{\v{v}_i} \]

Note that implicit integration introduces it's own artificial damping so you might find it's not necessary to add as much additional damping as you would with an explicit integration scheme.

I'll be going into more detail about implicit methods and FEM in subsequent posts, stay tuned!

Refs

[Baraff Witkin] - Physically Based Modelling, SIGGRAPH course
[Baraff Witkin] - Large Steps in Cloth Simulation
[N. Joubert] - An Introduction to Simulation
[D Prichard] - Implementing Baraff and Witkin's Cloth Simulation
[Choi] - Stable but Responsive Cloth
Numerical Recipes, 3rd edition 2007 - ch17.5