Lagrange-multiplier constraint for collinear beads in 3d?

[vputz]

Homework Statement

I'm trying to write a Lagrange-multiplier constraint simulation along the lines of those described by Andrew Witkin here.

The basic idea for this simulation is that you have a number of objects represented by three collinear beads joined by massless rods; one of the two rods changes length, and the other has a fixed length (the beads are in liquid interacting via rotne-prager tensor interactions, but the physics is not the gist of this question). The simulation works by calculating the Jacobian matrix of the constraints, using it to calculate Lagrange multipliers, finding constraint forces, and proceeding.

So we have two constraints:

1) the two beads of the fixed-length leg must remain the same distance apart, and
2) the three beads must remain collinear.

Here's the rub: I have the thing working wonderfully in 2d, but cannot for the life of me get it to behave in 3d.

For terminology, we'll label the beads 1-3 and say that leg 1-2 is changing length and 2-3 is remaining fixed.

Homework Equations

In 2d, the first constraint is that $$r_23 = \delta$$, which we can rewrite as $$C=r_{23}^2 - \delta^2 = 0$$, or

$$C=((x_3-x_2)^2+(y_3-y_2)^2=0$$

$$C = x_{3}^{2} - 2x_{2}x_{3} + x_{2}^{2} + y_{3}^2 - 2y_{2}y_{3} + y_2{^2} - \delta^2 = 0$$

So the Jacobian $$J_i = \frac{\partial C}{\partial q_{i}}$$ is

$$J=[ 0, 0, 2x_{2}-2x_{3}, 2y_{2}-2y{3}, 2x_{3}-2x_{2}, 2y_{3}-2y{2} ]$$
and we can remove the spurious 2's to get

$$J=[ 0, 0, x_{2}-x_{3}, y_{2}-y_{3}, x_{3}-x_{2}, y_{3}-y_{2} ]$$


For the second constraint in 2d (collinear), I used equivalent slopes:

$$\frac{y_{\alpha}}{x_{\alpha}} = \frac{y_{\beta}}{x_{\beta}}$$

$$C=y_{\alpha}x_{\beta}-y_{\beta}x_{\alpha} = 0$$

$$C=(y_{3}-y_{2})(x_{2}-x_{1}) - (y_{2}-y_{1})(x_{3}-x_{2})$$

$$C=y_{3}x_{2}-y_{3}x_{1}-y_{2}x_{2}+y_{2}x_{1}-y_{2}x_{3}+y_{2}x_{2}+y_{1}x_{3}+y_{1}x_{2}$$

and so the Jacobian is

$$J=[y_{2}-y_{3}, x_{3}-x_{2}, y_{3}-y_{1}, x_{1}-x_{3}, y_{1}-y_{2}, x_{2}-x_{1}]$$

(hopefully I've transcribed all that correctly). Okay--so far in 2d, all this works beautifully.

...(continued)

 

[vputz]

The Attempt at a Solution

In 3d, things start going wrong. We can basically extend constraint (1) fairly trivially:

$$J=[ 0, 0, 0, x_{2}-x_{3}, y_{2}-y_{3}, z_{2}-z_{3}, x_{3}-x_2, y_3-y_2, z_3-z_2]$$

...but the collinear constraint is trickier (to me). So far my attempts have centred around treating the problem as if we were constraining the spherical coordinates $$\theta$$ and $$\phi$$, since we can use the 2d "slope" constraint as our $$\theta$$ constraint, since it works and doesn't depend on z at all:

$$J=[y_{2}-y_{3}, x_{3}-x_{2}, 0, y_{3}-y_{1}, x_{1}-x_{3}, 0, y_{1}-y_{2}, x_{2}-x_{1}, 0]$$

Things go wrong when I attempt to add a third constraint for $$\phi$$. Attempts to do variations of $$z_{23}/r_{23}=z_{12}/r_{12}$$ fail--at best merely failing to constrain $$\phi$$ properly (at which point the two rods begin flexing like mad, which does interesting but nonphysical things to the simulation), and at worst interfering with the other constraints (in one case forcing the rods to elongate beyond the bounds of the simulation).

I'm not certain I understand what's going on. My best guess is that attempting to constrain it with separate constraints for $$\theta$$ and $$\phi$$ in cartesian coordinates is overconstraining the cartesian coords and going unstable. I'd like to find a nice way to either phrase it as one constraint equation or as three fairly independent ones. But I can't seem to manage it; I keep thinking in circles. Any good ideas?

Since this is an iterative computer simulation, keeping everything in cartesian is pretty much a must (there's a bunch of matrix math going on). Been hammering at this a while with no good results. The 3d version works just fine as long as everything is coplanar in z, but any variation in z makes it go berserk.

I even briefly tried leveraging $$A\cdotB=AB\cos\theta$$, setting $$\theta=0$$ to try and put everything on one constraint, which was almost interesting but since the Jacobian evaluated to zeros during the simulation turned out to be unuseful.