[ODE] More accuracy ???

[Gary R. Van Sickle]

> Nguyen Binh wrote:
> >     So, the code is modified to
> >
> >     ------Begin code fragment
> >     ......
> >         // handle linear velocity
> >         dReal h2 = h*h/2;
> >         for (j = 0; j < 3; j++)
> >                 b->pos[j] += h * b->lvel[j] + h2*b->facc[j];
> >     .....
> >     ------End
> >     So our ODE has O(h3) accuracy compare to just O(h2) with the old
> >     one.
>
> I've given it a try-out and it doesn't seem to hurt matters at all,
> though it's hard to say that it helps, and I'm not qualified to speak
> of its correctness.

I am ;-).

Assuming facc[] is constant over time interval h (which we have to assume, since
there's no other info to go on for the forces which get added into facc), we
have:

  a = v' = facc/Mass ==> constant
  v = last_v + integral of v' over time period h
    = last_v + (facc/Mass)*h
  p = last_p + p'
  p'= integral of v over time period h
    = last_v*h + ((facc/Mass)*h^2)/2


So, from what I can tell, Nguyen's code is missing a b->invMass, and should read
something like this:

     ------Begin code fragment
     ......
         // handle linear velocity
         dReal h2 = h*h/2;
         for (j = 0; j < 3; j++)
                 b->pos[j] += h * b->lvel[j] + h2*(b->facc[j]*b->invMass);
     .....
     ------End

There's two other problems though:

1.  At the point in the StepIsland code where moveAndRotateBody() gets called,
it's not clear to me what's been done to lvel[]; in particular, I see this
earlier in the code:

  for (j=0; j<3; j++) body[i]->lvel[j] += body_invMass * cforce[i*8+j];

cforce looks to be the constraint force.  So, lvel in the presence of a
constraint is != last_v in the above derivation.  What's that do to us?  I don't
know, it's too late, my brain shut down a half hour ago.  If there is no
constraint force though, this is the *exact* solution to the given problem, not
merely another Taylor-series factor.  Well, it is, but it's the *last* non-zero
factor.

2.  The forces acting on a body in the real world we're trying to model are
almost never (always never?) constant over a finite time period.  It would be
nice to get externally-controlled forces (e.g. springs, interbody gravitation,
etc.) inside the integration loop.  This is of course well beyond the scope of
Nguyen's patch, I just thought I'd mention it.

>  I find that rotational velocity is much less
> stable (I don't know about 'accurate', though the issues are
> probably related) than linear velocity -- is there a similar trick
> for rotational velocity?

Probably, but I think it would involve a bit of matrix and/or quaternion math.
I'm just not up to it at this hour. ;-)

--
Gary R. Van Sickle
Brewer.  Patriot.