ODE Solvers

[Mario Marengo]

Hi Jens,

I didn't want to hijack the other thread, so I'm starting a new one...

// mario: in case you read this, I would be intrested if you've implemented/worked with a verlet ODE solver; I wrote not too long ago a rigid body engine and for realtime simulations higher-order runge-kutta solvers were simply too slow. I menaged to speed up the collision detection / contact forces algos quite a bit, but I ended with euler as ODE-solver for the realtime part. (Though runge-kutta or other solvers with adaptive step control gave more predictable results and the solution is far from perfect). You know roughly how fast verlet is compared to runge-kutta of 3rd order ?!

Yes; Verlet is faster than RK4, which is why it is used extensively for modeling molecular dynamics -- i.e: lots'o'particles. But it's a little "apples and oranges" since RK is an explicit method. Like every other solver though, it's not a perfect solution for all occasions; it's got it's quirks. It is efficient if the system can be described without velocities: X''(t) = f(t,X(t))

I've used Verlet for a cloth simulation, and it worked really well (at the time, I compared it with Baraff's Midpoint, and decided to go with Verlet). The good thing about plain Verlet is that you can modify positional information directly (as a result of collisions, say) and the solver will adjust predictably (no need to manage velocities). But if some part of the system needs accurate velocities, you need to try something else (Midpoint or "Velocity Verlet" would be my choices).

Here are some of the comments from my (now 2-year-old ) implementation of these various solvers. I don't mind sharing the actual code, but I'd have to clean things up a bit, because a lot of it is tied to the rest of the library (they use generic containers, get built through a factory class, etc...)

This is Plain Verlet:

//-----------------------------------------------------------------------------
//                                  Verlet
//-----------------------------------------------------------------------------
//           Implements the Verlet method, whose update rule is:
//                X[t+h] = 2*X[t] - X[t-h] + A[t] * h^2
//           Where:
//           X    (vector) -> position
//           A    (vector) -> acceleration (F/m, or X'')
//           h    (scalar) -> delta-time (timestep)
//
//           Similar to RK2 or Midpoint.
//           WARNING! Due to the fact that velocity does not enter into
//           the equation explicitly, it should not be used in systems
//           that depend on an accurate V. In this implementation, V is
//           computed as the average velocity in interval [n-1,n+1]:
//                V[t] = (V[t+h]-V[t-h])/(2*h)
//           Which of course means that velocites are lagging behind
//           by one whole timestep!!
//
//           See VelVerlet or Leapfrog for a (slightly more
//           expensive) version that does a much better job with
//           velocities.
//
//           OTOH, if velocities aren't crucial, this method rocks! :)
//
//-----------------------------------------------------------------------------

Note that it needs the position from the previous step. This influences implementation of course, but it's not too much of a headache. Another neat feature of this solver is that it is reversible -- i.e: it can solve forward or backward in time.

Here's a variation that *does* use velocity. I think it's usually referred to as "Velocity Verlet", but I'm not sure:

//-----------------------------------------------------------------------------
//                             Velocity Verlet
//-----------------------------------------------------------------------------
//
//         Implements the Velocity Verlet method, whose update rule is:
//                  X[t+h] = X[t] + h*V[t] + 0.5*h^2 * A[t]
//                  V[t+h] = V[t] + 0.5*h * (A[t] + A[t+h])
//         Where:
//         X    (vector) -> position
//         V    (vector) -> velocity (or X')
//         F    (vector) -> force
//         A    (vector) -> acceleration (or F/m, or X'')
//         m    (scalar) -> mass
//         h    (scalar) -> delta-time (timestep)
//
//         Unlike the vanilla Verlet, this one can more accurately
//         compute V[t+h]. However; in order to do this, it needs
//         to evaluate F at X[t+h], which makes it more expensive
//         to compute. Exactly "how much" more expensive, depends
//         on the cost of evaluating F.
//         <snip>
//-----------------------------------------------------------------------------

As you can see, this one doesn't require the client to provide old positions, so the implementation is simpler, and uses less memory, etc.

Let me know if you want the actual implementation, and I'll put something together.

Actually; looking through this stuff, I notice I have two implementations: one as an ODESolver object hirearchy, and another as free-floating function templates (which I could extract much more easily).

Cheers!

[edit] Ooops; forgot to mention: as you know, there are litterally thousands of references for ODEs out there, but I came across one in particular, where they do a more in-depth comparison analysis. Check it out:

http://www.gris.uni-tuebingen.de/~mhauth/t...29lNumerics.pdf

[TheDunadan]

Whoa, didn't expect such detailed answer, but glad you found the time

Unless I've missed something crucial the verlet solver looks pretty easy to implement and it shouldn't take much time to integrate it in my little rigidbody engine. If I should run into any troubles I'll come back to your offer for the source code, but I wanna give it a try myself first.

Once I cleaned up the code for the rb engine I intend to make it public and anyone intrested can have a look at it. Now it's in a little bowling simulation game and the players can drop in all sorts forces and see how the physics change. (And there are 7 different ODE solvers to choose from: the verlet one will be #8 )

there are litterally thousands of references for ODEs out there

Yep, and the hardest part is to check wich ones are worth reading and which not . Concerning the paper you linked here, doesn't surprise me at all that Ron Fedkiw has his hands in here, lots of intresting papers on his site. I hope to find some time implementing his collision propergation algorithm for rigid body engines.... but there is always this contant lack of time

In magnis et voluisse sat est

Jens

btw. have I ever mentioned I like the odforce smilies

[Mario Marengo]

Unless I've missed something crucial the verlet solver looks pretty easy to implement and it shouldn't take much time to integrate it in my little rigidbody engine. If I should run into any troubles I'll come back to your offer for the source code, but I wanna give it a try myself first.

Nope; I don't think you've missed anything; the update itself is almost as trivial as Euler. The only slight complication comes with mantaining the cache for the "previous sample" and the logistics associated with that. Beyond that, it's all about the data structures of the system itself: are the particles stored as a single array of N-dimentional vectors; or as separate dynamically sized arrays of arbitrary attributes; or... the usual million-and-one conciderations. But the ODE solver itself is simple

.....but there is always this contant lack of time  

Amen.

In magnis et voluisse sat est

You made me go and look that up!

Spanish, no problem; Latin.....hmmmmnot so much...

[TheDunadan]

You made me go and look that up!

I just tried to find something witty and by accident found my old latin book Our teacher loved those latin sayings and I grew fond of them.

Jens

Felix qui potuit rerum cognoscere causas *couldn't resist *