A mass spring model for hair simulation

A strand of hair is represented by a set of particles connected by springs.

The velocity for a particular particle is calculated implicitly using the following formula: $\boldsymbol{v}^{n+1/2}=\boldsymbol{v}^{n}+\frac{\Delta t}{2}\boldsymbol{a}(t^{n+1/2},\boldsymbol{x}^{n},\boldsymbol{v}^{n+1/2})$

The force or acceleration ($\boldsymbol{a}$ in the above equation) produced by the spring between two adjacent particles is given by the following: $\boldsymbol{F}^{n+1}=\frac{k}{l_{0}}\left((\boldsymbol{x}_{2}^{n}-\boldsymbol{x}_{1}^{n})^{\mathrm{T}}\hat{\boldsymbol{d}}^{n}-l_{0}\right)\hat{\boldsymbol{d}}^{n}+\Delta t\frac{k}{l_{0}}(\boldsymbol{v}_{2}^{n+1}-\boldsymbol{v}_{1}^{n+1})^{\mathrm{T}}\hat{\boldsymbol{d}}^{n}\hat{\boldsymbol{d}}^{n}$

where $\hat{d}^{n}=(x_{2}^{n}-x_{1}^{n})/\left\Vert x_{2}^{n}-x_{1}^{n}\right\Vert$

My question is, how do you numerically calculate $\boldsymbol{v}_{1}^{n+1}$ and $\boldsymbol{v}_{2}^{n+1}$ in order to calculate the force?

I've tried using Newton's Method but calculating derivative the spring force is just so complicated.

I'm attempting to implement the techniques found in this paper: link

-
+1 for the cool animations on the linked page. – Emil Jeřábek May 23 '11 at 16:10
I'm unsure why this has been deemed off-topic, since it's a legitimate question of research in scientific computing/graphics. Maybe MO is not the right forum for these topics. Anyways. In the paper, the authors suggest an extrapolation step from the $v^{n+1/2}$ stage to $v^{n+1}$, along with a (nebulously defined) post-processing. – Nilima Nigam May 23 '11 at 19:58
@Michael, I did indeed start a discussion on meta about this. – Nilima Nigam May 23 '11 at 20:22
Meta discussion is here: tea.mathoverflow.net/discussion/1052 – S. Carnahan May 24 '11 at 5:01
Sebastian, you should consider editing the question to make it more widely accessible (I don't have enough points to edit it myself). Specify the ODE, the reduction to first order, clarify the notation (especially the $t_{n+1/2}$ and $t_{n+1}$). Specify the method is explicit in the (nonlinear) force update. You should not need to use a Newton update for the velocity, since the authors claim this can be avoided. They are recommending extrapolation. You can ask the question: what is an accurate and efficient extrapolation strategy for the velocity in this algorithm? – Nilima Nigam May 24 '11 at 17:34

Notice that the only unknown quantities in your force equation are the updated velocities $v_i^{n+1}$, and that force depends linearly on these. Everything else is either known from the last time step, e.g. $x_i^n$, or is a simulation parameter/material constant.
Therefore, after substituting for $\mathbf{a}$ in your first equation, you will end up with a set of linear equations in the unknowns $v_i^{n+1/2}$; in particular the system will be of the form $$(\mathbf{I} - c\Delta t^2\ \mathbf{d}\mathbf{d}^T)\mathbf{v}^{n+1/2} = \mathbf{b}$$ for some scalar $c$ and right-hand side $\mathbf{b}$. You can readily solve this system without Newton's method. In particular, for sufficiently small time steps $\Delta t$ the matrix is positive-definite, in which case I recommend using Conjugate Gradients to solve the system iteratively without ever even needing to form the dense matrix $\mathbf{d}\mathbf{d}^T.$