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](http://physbam.stanford.edu/~aselle/papers/9/)



A mass spring model for hair simulation
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.$


