[ODE] ODE Source code vs. Baraff

[david@csworkbench.com]

Russ (or anyone else who has a clue),

I've been reading up a bit more on this algorithm, trying to figure out
why my latest round of optimizations didn't work as expected, and I
noticed something I wasn't sure about..... It probably has something to do
with the fact that it's 2 AM here, but I'll ask it anyway in case the
cobwebs haven't cleared by morning (they usually don't).  According to
Baraff, the whole of the algorithm can be summed up as:

A*lamda = b (plus some stuff for the LCP solver to simulate friction and
limits..)
where:
invM = inverse of Mass matrix, Fext = external forces (and torques)
A = J * invM * J'  (I got this one!)
b = -(J * invM * Fext + c) (here's my problem)

What I see in the code is:

b = c/stepsize - (J * (invM * Fext + vel/stepsize))

Ignoring the difference in sign (I'm sure I just missed it somewhere), if
I understand correctly, this is the difference between yours and Baraff's
implementations.  You calculate the amount of force to be applied over the
step, and he calculates instantaneous forces.

What I was hoping was, if I only updated the Fext accumulators inside the
iteration loop, then I would only have to recalculate b (rhs) each
iteration, and lamda would become progressively smaller.  I'm
interpretting lamda (well, J * lamda) as the force required to move the
system from the state implied in rhs to a state consistent with the
constraints defined by A.  What is happening is, lamda is not changing
(much), and iterating is only giving me that many times as much force as
should have been added.  Am I misinterpretting lamda, or is its definition
slightly different in your implementation?  If so, how can I reformulate
the problem so that:

for i = 0 to iterations
  for each joint in island (world for now)
    calculate rhs = f(body->facc, body->tacc)...
    dSolveLCP(A, lamda, rhs, lo, hi, other stuff)
    add J * lamda to accumulators
  next joint
next i

assuming A has been precalculated for each joint separately, so that the
accumulators would eventually converge on the global max.  All I really
need is a formula for rhs that makes lamda mean the amount of force that
needs to be added to the accumulators (that were used to calculate rhs in
the first place....) to bring the two bodies this joint connects to a
consistent state.  It sounds like, from Baraff's linear time dynamics
paper, that if rhs only computes instantaneous forces for Fext, that the
resulting Ax=b equation can be iteratively solved, but I didn't get that
far before my brain ran out of steam.

If what I'm looking for is mathematically impossible (or slow), I've still
got a copy of the first (well, second) algorithm that I tried and works
fairly well.  But this one is bound to be faster, if it works, just
because you only have to calculate A once per step, not once per
iteration.  (I tried using the first A generated as an approximation for
all the iterations, and had to bump the number of iterations up to keep
the same stability -- no net speed increase.)

Any thoughts would be helpful,
David Whittaker

P.S. I haven't taken Diff. Eq. yet, and my Matrix Theory professor
struggled with English....  If you answered my question in one of those
last couple emails, I missed it.