[ODE] Other integrators (was: ODE implementors guide)

[nlin@nlin.net]

> > but this separation is not hard, and after all an integrator is
> > really a small piece of code compared to the complexities of
> > calculating the system matrix.
> 
> is it possible to do that now (time consuming ?).
> It would be amazing to test new features (integrators and hybrids).

After looking through the code, here is my understanding of how the
integration works. I'd appreciate any comments from experts, in particular
those with last names of Smith :-) 

First, the vector of Lagrange multipliers for the constraint forces
is calculated in dInternalStepIsland_x2 with the call to "dSolveLCP",
handling both "normal" constraint forces (unbounded) as well as
LCP constraint forces (contact forces) (but what are the mixed
unbounded+LCP forces?)

Given the vector lambda, the constraint force is calculated in 
dInternalStepIsland_x2 in the code marked, surprisingly, with the
comment "compute the constraint force".

Given the constraint force, it is added to fe (the external force)
to get the total, constraint-satisfying, force on the body and then,
in the section of code marked (still in dInternalStepIsland_x2)
"compute the velocity update", the constraint force is directly integrated
over the entire timestep (since a "time-stepping method" uses integrals of
force over the entire timestep, not instantaneous impulsive forces at the
beginning or end of timesteps) by multiplying it (force) with the stepsize,
dividing it by body mass (since F=m*a => a=F/m), and adding the result
(F/m multiplied by the timestep length) to the body velocity, both linear
and angular. This is an Euler integration step.

With the new linear velocity and angular velocity, the new linear
position and angular orientation are computed in function 
"moveAndRotateBody", where the linear and angular velocites are
multiplied by the stepsize to obtain the new position and orientation.
This is also an Euler integration step.

Given this understanding, it appears that to change this to use a higher
order integrator you would need to:

1. In dInternalStepIsland_x2, don't merely multiply the computed
   acceleration by the stepsize, but instead use the currently
   computed acceleration along with previously computed accelerations,
   saved from previous timesteps, to use a higher order integration
   scheme (BUT!! see note below on time-stepping)
2. In moveAndRotateBody, don't merely multiply the computed velocity
   by the stepsize, but instead use the currently computed velocity
   along with previously computed velocites, saved from previous timesteps,
   to use a higher order integration scheme.

The problem with all of this (assuming the above is correct), and where I
hope some experts will comment, is the interplay of all of this with
the time-stepping scheme. The accelerations computed by ODE (via the call
to dSolveLCP), as far as I understand, are not INSTANTANEOUS accelerations,
but are instead accelerations which must be applied (integrated) over an entire
timestep. I think that the method of solving for the acceleration depends on
the fact that it will be applied (integrated) over the timestep with an
Euler integration method.

Is the above a problem, or not?

If it is a problem, then exactly how and where does the solution of the
accelerations depend on the Euler integration step, and how would one need
to modify it to work under the assumption of a different integrator?
Do any similar assumptions underlie the positional update? In general,
how does the time-stepping formulation depend on the integrator?

Hope this helps some; I'm still trying to learn ODE myself. If you are
feeling bold you could simply dive right in and replace the mentioned code
with higher-order integrator; but I have no idea about the repercussions
this would have on the stability/correctness of the system, because higher
order integrators might violate some of the assumptions used in the
time-stepping formulation. David Stewart's articles talk a lot about
time-stepping, but they require some time to digest. In the article I
quoted in my last message ("Time stepping methods and the mathematics
of rigid body dynamics", see last message for link), on p. 19, the
the rigid body problem over one timestep is formulated as an LCP where
the Euler integration assumption is already contained within the
formulation (see Equation 1.26). This appears to imply that the LCP
which must be solved at each timestep depends on the integration scheme
which is later used to apply the computed accelerations/velocities; a
different integrator appears to need a different formulation of the
LCP.

-Norman