🔒 Locked

Thoughts on velocity verlet

Started by Numsgil Feb 10, 2006 at 4:32 PM 22 replies 21.8k views

Original Post

Feb 10, 2006 at 4:32 PM

I've been implementing a somewhat simple physics engine for a project I've been working on. Since about last June I've been researching different ways of doing it. I started with Euler, then moved on to Verlet. Both have problems, though with different things. Verlet has issues with things like velocity limiting and collisions (specifically, making them respect proper physics rules, like conservation of momentum). Euler has issues with rigid constraints, and alot of instability can occur if forces become rather large. Dissatisfied, I decided to look for more schemes to use. I found references to "leap frog" and "velocity" versions of verlet in the Jakobsen article. Leap frog didn't really appeal to me, but velocity verlet certainly did. This short article describes the method Since you have an explicit velocity term, velocity limiting and collision response become as easy as they were in Euler. Since it's a verlet scheme, it inherits all the stability and rigid constraint advantages of the more regular verlet method. I've already done the simpler aspects of my physics demo (collision response between spheres, rigid in-bounds constraints, etc.) in velocity verlet, and I'm quite pleased. I simply don't see anything that would indicate it's inferior in any way to either the regular verlet or Euler. Why then can't I find any game-related literature on it? Every reference I find to it is related to particle simulations in legitimate physics research. Things like explosion simulators and atomic interaction simulations. Is there some hiccup to using velocity verlet I'm not seeing? Maybe everyone is just following that Jakobsen article without doing any further research?

ozmic66

Feb 11, 2006 at 10:52 PM

looks like nobody has an answer :0

i'm wodering what ppl have to say too

Sneftel

Feb 11, 2006 at 11:09 PM

"Velocity Verlet" integration is essentially identical to regular Verlet integration. The only difference is that instead of keeping x(t) and x(t-1) around, you keep x(t) and x(t)-x(t-1) around (as v(t)). It makes velocity-limiting clearer, not easier.

Have you considered high-order Runge-Kutta integration? Tastes like Euler, kicks like Verlet.

Numsgil

Feb 12, 2006 at 12:07 AM

That's a bold assertion Sneftel. I disagree.

They are mathematically equivellant, and have similar orders of convergence. However, in practice they work quite differently, as my experience at the very least can attest.

Take for instance collisions. In regular verlet, since you don't have explicit control over velocity, correct collision response becomes difficult. Anyone whose ever tried to use verlet as their main integration method for their physics engine can probably attest to this, I found such problems discussed quite often during my research.

Verlet is absolutely beautiful at handling projections and other constraints. Push the ball off the wall, keep the stick rigid, etc. It fails, however, at tasks that are less rigid.

One account in particular (which I can't seem to find at the moment, but I may be able to look up if anyone is interested) describes how a programmer became enamored with verlet after reading Jakobsen's article, and suggested using it to his boss. They ran into serious problems down the road and had to use several hacks to get the physics to look right in their project, which the author was most displeased with. I was having similar problems.

Velocity verlet has several advantages over vanilla verlet, addressing the pitfall I described above:

Source 1
Source 2
Source 3

note that none of the above really seem to be game-related, which is what I really find curious. Game programmers seem to have really dropped the ball IMO. I would expect to find even a reference somewhere, something to show that someone involved with video games tried it even once.

To sum it all the sources up, velocity verlet allows for more precise treatments of velocities (absolutely a necessity for proper collision response) and updates all variables from the previous time step only, as opposed to vanilla verlet which updates from the previous two time steps. It does require more vectors to be stored than vanilla verlet, and is slightly more computationally expensive.

edit:
Oh, I forgot to respond to the second part of your post.

I did try using the Runge-Kutta method. I found it overkill for my purposes, and too slow. At least compared to Euler and Verlet, etc. I could see it being useful if you're running a very precise simulation, where numerical accuracy is paramount, but I think that's the exception instead of the rule for game physics.

2nd edit:
This article describes the problems of using regular verlet

Apparently the velocity term is far more inaccurate than the position term.

[Edited by - Numsgil on February 12, 2006 12:07:51 AM]

John Schultz

Feb 12, 2006 at 12:47 AM

Quote:
Original post by Numsgil
Take for instance collisions. In regular verlet, since you don't have explicit control over velocity, correct collision response becomes difficult. Anyone whose ever tried to use verlet as their main integration method for their physics engine can probably attest to this, I found such problems discussed quite often during my research.

Try this test program. The first version is really Newton-Stormer-Verlet (NSV), AKA "Semi-implicit" Euler, as velocity is computed first and used to update position (traditional Euler has the order reversed).

All of those methods compute the same result (within the limits of float precision). When using velocity-less Verlet, updating velocity explicitly is straightforward: back project the old position to directly update the effective velocity.

I would agree that NSV with stability clamping for spring-like forces and potentially unstable impulses is superior to RK2/RK4 for typical game applications. However, RK and similar higher order methods provide better behavior for high angular velocities (and for some periodic force systems with larger time steps). Also keep in mind that the higher order methods are not always more accurate (they can lose energy, and may not be symplectic), though they are generally more stable at a particular time step.

As you may find, there are many formulations for Verlet; some with RK-like properties while maintaining symplectic behavior (important for some applications). The higher order methods will require more derivative evaluations and more float ops: it's a tradeoff.

Numsgil

Feb 12, 2006 at 1:13 AM

Quote:

All of those methods compute the same result (within the limits of float precision). When using velocity-less Verlet, updating velocity explicitly is straightforward: back project the old position to directly update the effective velocity.

Yes, however there is significantly more error in this velocity term than in the position term, as my second edit above explains. I don't know if this error term is significant enough to matter in games necessarily, but you add to it the conceptual overhead of not having an explicit velocity term, and I don't see any significant advantage at all.

It's all a matter of your needs of course. Most games, however, generally follow similar needs. One of those needs is having a velocity term for several different calculations (collisions in particular). You'll probably be saving the velocity term to a variable instead of recalculating it several times if you use vanilla verlet.

If that's the case, vanilla verlet doesn't save you any vectors in storage. Even if that wasn't the case, a measly one extra vector to give you explicit control over velocity seems worth it to me at least. Admitedly, the implementation is slightly more complex, but I doubt it's above the head of anyone who understands verlet.

There are a whole family of verlet-esque algorithms that don't really seem explored, which I guess is what I have issue with. They all have strengths/weaknesses, (though some are clearly superior to others in all but extreme cases) but any articles game-related on physics only makes a passing reference to leap-frog and velocity versions without exploring why the author picked vanilla verlet.

Case in point:

Article on Molecular Dynamics Simulations

Near the bottom, there is a small passage on integration techniques, exploring the pros/cons of each method. It also lists something called "Beeman's algorithm" which I haven't found mentioned even in passing on any game sites.

This stuff isn't hard to find, I'm just wondering why no one's made a survey of it for game physics specifically.

John Schultz

Feb 12, 2006 at 1:34 AM

Quote:
Original post by Numsgil
Quote:

All of those methods compute the same result (within the limits of float precision). When using velocity-less Verlet, updating velocity explicitly is straightforward: back project the old position to directly update the effective velocity.

Yes, however there is significantly more error in this velocity term than in the position term, as my second edit above explains. I don't know if this error term is significant enough to matter in games necessarily, but you add to it the conceptual overhead of not having an explicit velocity term, and I don't see any significant advantage at all.

From experience, the following methods are simple and fast, with no major issues for simple game applications when used with force/acceleration clamping (including very stiff springs for vehicle suspensions). Generalized constraints are most easily handled using big-matrix (semi-implicit (solved simultaneously) optimization-based: LCP/Jacobi/PGS) formulations (though I have so far gotten away with a pure impulse-with-tricks methodology):

v += a*dt
x += v*dt

or

v += a*dt*dt
x += v // v is prescaled, now a displacement

First best results, use a fixed time step. The following is then possible:

dt2 = dt*dt

...

v += a*dt2
x += v // v is prescaled, now a displacement

dinnyesg

Feb 12, 2006 at 2:20 AM

Numsgil:

I also read the Jacobsen paper and started to experiment with verlet physics. For my simple game collision response is not that difficult because I don't need bounces. Still I am curious how you managed to handle bouncing collisions with the newly introduced velocity term. For a single point particle I know it works. But what do you do if you have a more complex object falling to the ground and penetrating into it only at one point. In this case if you relfect the velocity for bouncing collision this only affets that single point and not the whole body. The result is that the other particles push that particle back and you don't have bouncing. The more other particles the less bouncing. How did you approach this problem?

Numsgil

Feb 12, 2006 at 2:22 AM

Quote:
Original post by John Schultz
Quote:
Original post by Numsgil
Quote:

All of those methods compute the same result (within the limits of float precision). When using velocity-less Verlet, updating velocity explicitly is straightforward: back project the old position to directly update the effective velocity.

Yes, however there is significantly more error in this velocity term than in the position term, as my second edit above explains. I don't know if this error term is significant enough to matter in games necessarily, but you add to it the conceptual overhead of not having an explicit velocity term, and I don't see any significant advantage at all.

From experience, the following methods are simple and fast, with no major issues for simple game applications when used with force/acceleration clamping (including very stiff springs for vehicle suspensions).

There is something to say for straightforwardness and the ability to plop down off the top of your head. That's why bubble-sort still exists (though not without controversy).

However if Euler alone was adequate, I don't think that Jakobsen article would even be around. There are solutions that are moderately more complex and generally more stable, meaning the range of forces you have to clamp is diminished.

We aren't even talking hundreds or even tens of lines of code. The code from my engine using velocity verlet as it stands (it uses a time step of 1):

pos += vel + 0.5f * oldImpulse / mass;
vel += (oldImpulse + Impulse) / (2 * mass);
oldImpulse = Impulse;

The same code would look in Euler like:

vel += Impulse / mass;
pos += vel;

The difference is so minor, I don't see why you wouldn't take the (I'm assuming) better method on a larger project. There is a slightly higher computational overhead, but I can't imagine it actually being a bottleneck for any modern game. Especially since the code in my engine could be further optimized if need be.

Quote:

Generalized constraints are most easily handled using big-matrix (semi-implicit (solved simultaneously) optimization-based: LCP/Jacobi/PGS) formulations

A little off topic, but are there any resources for constraint solving in games other than relaxation shown in the Jakobsen article? My matrix skills are a little rusty, and most of the methods I find described in math articles are matrix based. I'm sure the final implementation of a Jacobi or other technique would involve matrices, but some source is always nice ;)

John Schultz

Feb 12, 2006 at 3:57 AM

Quote:
Original post by Numsgil
pos += vel + 0.5f * oldImpulse / mass;
vel += (oldImpulse + Impulse) / (2 * mass);
oldImpulse = Impulse;

I experimented with similar methods (effectively applies a low-pass filter on the velocity update) while trying to get cheap RK-like (multi-derivative/higher order) methods by staggering the calculations and averaging.

Quote:
Original post by Numsgil
The same code would look in Euler like:

vel += Impulse / mass;
pos += vel;

If you have a simple demo that shows the first method provides better behavior than the second, other developers might find it interesting.

Quote:
Original post by Numsgil
A little off topic, but are there any resources for constraint solving in games other than relaxation shown in the Jakobsen article? My matrix skills are a little rusty, and most of the methods I find described in math articles are matrix based. I'm sure the final implementation of a Jacobi or other technique would involve matrices, but some source is always nice ;)

See the ODE and Doom III source code. See also the Bullet forum (active discussion on LCP, PGS, etc.).

jjd

Feb 12, 2006 at 8:22 AM

Quote:
Original post by Numsgil
They are mathematically equivellant, and have similar orders of convergence. However, in practice they work quite differently, as my experience at the very least can attest.

No, they are definitely not mathematically equivalent. Regular Verlet is higher order in displacement, velocity-Verlet is higher order in velocity. Regular Verlet has the disadvantage ov being multistep (not self starting), and velocity-Verlet requires an additional call to the forcing function each iteration. As John Schultz alluded too, there are pros and cons to using each and which one is superior probably depends upon your application.

Numsgil

Feb 12, 2006 at 10:28 AM

Quote:
Original post by jjd
Quote:
Original post by Numsgil
They are mathematically equivellant, and have similar orders of convergence. However, in practice they work quite differently, as my experience at the very least can attest.

No, they are definitely not mathematically equivalent.

They are mathematically equivelant according to David Kofke of the Department of Chemical Engineering at SUNY Buffalo. See my "source 3" above (which are his lecture notes apparently). It may depend on exactly what you mean by "mathematically equivelant" of course. I'm assuming Kofke meant they are algebraic transformations of one another.

Do you have any sources for your assertion that vanilla verlet has better positional accuracy than velocity verlet?

Quote:

As John Schultz alluded too, there are pros and cons to using each and which one is superior probably depends upon your application.

I agree, but again why hasn't there been any examination (at least that I can find) of the different methods for the purposes of video game physics, which has different needs than molecular dynamics like most of my sources? Games all have very similar physics needs (Except maybe racing games (I just haven't done any looking into racing physics, so I won't care to speculate)), if not, then the Jakobsen article wouldn't be quoted umpteen-million times.

Has anyone used an integrating method in their project other than Euler, vanilla verlet, or Runge-Kutta? Has anyone even thought to look for another method?

Quote:

If you have a simple demo that shows the first method provides better behavior than the second, other developers might find it interesting.

I'm assuming by this you're implying that there isn't a difference. To be honest, I've only really examined velocity verlet vs. vanilla verlet. I was turned off on Euler early on while using it because it was incapable of handling really tiny masses. Something like 1/10000 of the mass of more normal objects caused the simulation to explode.

My simulation does not seem to explode with velocity-verlet. That's hardly a unbiased comparison, however, since the engine isn't as done as the Euler one was, and some things have changed in the interim.

Quote:

See the ODE and Doom III source code. See also the Bullet forum (active discussion on LCP, PGS, etc.).

Thanks, I'll give them a look.

jjd

Feb 12, 2006 at 12:19 PM

Quote:
Original post by Numsgil
They are mathematically equivelant according to David Kofke of the Department of Chemical Engineering at SUNY Buffalo. See my "source 3" above (which are his lecture notes apparently). It may depend on exactly what you mean by "mathematically equivelant" of course. I'm assuming Kofke meant they are algebraic transformations of one another.

Do you have any sources for your assertion that vanilla verlet has better positional accuracy than velocity verlet?

I used to be a member of the faculty in the dept of applied mathematics and the university of colorado in boulder. I'm not too interested in the claims of a chemical engineer on the mathematical equivalence of two numerical integrators.

The Verlet method is a special case of a Stormer method. It is simply derived by adding two Taylor series in the following way. Suppose you have a differential equation of the form x'' = f(t,x(t)), where x(t) is smooth. Using a Taylor series we can write

x(t+h) = x(t) + h x'(t) + h^2 x''(t)/2 + h^3 x'''(t)/6 + O(h^4)

where h > 0. In other words, we take a small step forward in time. We can also take a small step backward in time resulting in the following series

x(t-h) = x(t) - h x'(t) + h^2 x''(t)/2 - h^3 x'''(t)/6 + O(h^4)

It is important to note that the terms with odd-order derivatives have different signs but those of even-order do not. Add these two series together and you get

x(t+h) + x(t-h) = 2 x(t) + h^2 x''(t) + O(h^4)

i.e. all odd-order terms cancel. From the original equation (and with a little rearrangement) we get

x(t+h) = 2 x(t) - x(t-h) + h^2 x''(t) + O(h^4)

Thus Verlet's algorithm has a truncation error of O(h^4). To get an approximation of the velocity, v(t), you can use

v(t) = (x(t) - x(t-h))/h + O(h)

The 'velocity' Verlet algorithm begins with the Taylor series

x(t+h) = x(t) + x'(t) + h^2 x''(t) + O(h^3)

Using the relation v(t) = x'(t), and the original equation, we get

x(t+h) = x(t) + v(t) + h^2 f(t)/2 + O(h^3)

'velocity' Verlet truncates the h^3 and higher-order terms here, so the truncation error scales like O(h^3). Here we see that the displacement in 'velocity' Verlet is one order less than the original.

The velocity is calculated in a way similar to the original Verlet method. We write out the Taylor series for the velocity

v(t+h) = v(t) + h f(t) + O(h^2)

We can also write a Taylor series from v(t+h) back to v(t)

v(t) = v(t+h) - h f(t+h) + O(h^2)

which we can also write as

v(t+h) = v(t) + h f(t+h) + O(h^2)

Combining these equations gives

v(t+h) = v(t) + (h / 2) ( f(t) + f(t+h) ) + O(h^2)

Note that this formula is implicit if the forcing depends on the velocity. Note that the truncation error is O(h^2) or one order higher than in the original Verlet method. When the forcing depends only of the displacement, the force can be calculated before calculating the velocity. Thus we get the steps

x(t+h) = x(t) + h v(t) + h^2 f(t,x(t)) / 2
f(t+h) = f(t+h,x(t+h))
v(t+h) = v(t) + (h / 2) ( f(t) + f(t+h) )

Forgive my abuse of notation, I think the meaning is clear though.

I suspect that the reason only a few methods are commonly used in games reflects the history and requirements that are particular to games:
(1) kickass physics engines are a relatively recent development
(2) physical accuracy is not so important as it is in the physical sciences
(3) the realtime demand on games is a severe constraint, although one that will probably lessen as computers become more powerful.

Games require fast, stable methods. Often the amount of memory used is a constraint too. This limits the methods that people are willing to explore. Implicit methods would be great, if only they didn't come with a high cost to memory. Multi-step methods suffer a similar fate and the additional problem of start-up problems and that changing the size of the timestep is awkward. As such, explicit methods with few function evaluations are currently the best option for most games. I believe that will change in time, and I hope that more use can be made of implicit methods so that my ragdolls will behave more stably [smile]

Oh, almost forgot. You wanted some sources "the art of molecular dynamics simulation" by rapaport; "numerical recipes" by press et al.; "solving ordinary differential equation I and II" by hairer and wanner. I also really like "computer methods for ordinary differential equations and differential-algebraic equations" by ascher and petzold.

John Schultz

Feb 12, 2006 at 12:43 PM

struct Test {  Test() {    int count = 100;    flt x = 0.f;    flt v = 0.f;    flt a = 9.81f;    flt dt = 1.f/flt(count);    int i;    lprintf("Simple Euler:\n");    for (i=0; i < count; i++) {      x += v*dt;      v += a*dt;      lprintf("V: %f X: %f\n",v,x);    } // for    x = v = 0.f;    lprintf("More Accurate Euler:\n");    for (i=0; i < count; i++) {      x += v*dt + .5f*a*dt*dt;      v += a*dt;      lprintf("V: %f X: %f\n",v,x);    } // for    x = v = 0.f;    lprintf("NSV: Newton-Stormer-Verlet ('Semi-implicit' Euler):\n");    for (i=0; i < count; i++) {      v += a*dt;      x += v*dt;      lprintf("V: %f X: %f\n",v,x);    } // for    lprintf("Velocity-less Verlet:\n");    flt xc = 0.f;    flt xo = xc;    for (i=0; i < count; i++) {      v = xc - xo + a*dt*dt;      xo = xc;      xc += v;      lprintf("V: %f X: %f\n",v/dt,xc);    } // for    lprintf("My NSV variant\n");    x = v = 0.f;    flt dt2 = dt*dt;    for (i=0; i < count; i++) {      v += a*dt2;      x += v; // v is prescaled: really a displacement.      lprintf("V: %f X: %f\n",v/dt,x);    } // for    lprintf("Velocity Verlet variant\n");    flt oldAccel = a; //0.f; // May be self-starting issues...    x = v = 0.f;    for (i=0; i < count; i++) {      x += v*dt + .5f*oldAccel*dt*dt;      v +=        .5f*(oldAccel+a)*dt;      oldAccel = a;      lprintf("V: %f X: %f\n",v,x);    } // for  }} test;

Output:

Simple Euler:V: 0.098100 X: 0.000000V: 0.196200 X: 0.000981V: 0.294300 X: 0.002943V: 0.392400 X: 0.005886V: 0.490500 X: 0.009810V: 0.588600 X: 0.014715V: 0.686700 X: 0.020601V: 0.784800 X: 0.027468V: 0.882900 X: 0.035316V: 0.981000 X: 0.044145V: 1.079100 X: 0.053955V: 1.177200 X: 0.064746V: 1.275300 X: 0.076518V: 1.373400 X: 0.089271V: 1.471500 X: 0.103005V: 1.569600 X: 0.117720V: 1.667700 X: 0.133416V: 1.765800 X: 0.150093V: 1.863900 X: 0.167751V: 1.962000 X: 0.186390V: 2.060100 X: 0.206010V: 2.158200 X: 0.226611V: 2.256299 X: 0.248193V: 2.354399 X: 0.270756V: 2.452499 X: 0.294300V: 2.550599 X: 0.318825V: 2.648699 X: 0.344331V: 2.746799 X: 0.370818V: 2.844899 X: 0.398286V: 2.942999 X: 0.426735V: 3.041099 X: 0.456165V: 3.139199 X: 0.486576V: 3.237299 X: 0.517968V: 3.335399 X: 0.550341V: 3.433499 X: 0.583695V: 3.531599 X: 0.618030V: 3.629699 X: 0.653346V: 3.727799 X: 0.689643V: 3.825899 X: 0.726921V: 3.923999 X: 0.765180V: 4.022099 X: 0.804420V: 4.120199 X: 0.844641V: 4.218299 X: 0.885843V: 4.316399 X: 0.928026V: 4.414499 X: 0.971190V: 4.512599 X: 1.015335V: 4.610700 X: 1.060461V: 4.708800 X: 1.106568V: 4.806900 X: 1.153656V: 4.905000 X: 1.201725V: 5.003100 X: 1.250775V: 5.101201 X: 1.300806V: 5.199301 X: 1.351818V: 5.297401 X: 1.403811V: 5.395501 X: 1.456785V: 5.493601 X: 1.510740V: 5.591702 X: 1.565676V: 5.689802 X: 1.621593V: 5.787902 X: 1.678491V: 5.886002 X: 1.736370V: 5.984102 X: 1.795230V: 6.082202 X: 1.855071V: 6.180303 X: 1.915893V: 6.278403 X: 1.977696V: 6.376503 X: 2.040480V: 6.474603 X: 2.104245V: 6.572703 X: 2.168991V: 6.670804 X: 2.234718V: 6.768904 X: 2.301426V: 6.867004 X: 2.369115V: 6.965104 X: 2.437785V: 7.063204 X: 2.507436V: 7.161304 X: 2.578068V: 7.259405 X: 2.649681V: 7.357505 X: 2.722275V: 7.455605 X: 2.795850V: 7.553705 X: 2.870406V: 7.651805 X: 2.945943V: 7.749906 X: 3.022461V: 7.848006 X: 3.099961V: 7.946106 X: 3.178441V: 8.044206 X: 3.257902V: 8.142305 X: 3.338344V: 8.240405 X: 3.419767V: 8.338505 X: 3.502171V: 8.436604 X: 3.585556V: 8.534704 X: 3.669922V: 8.632804 X: 3.755269V: 8.730904 X: 3.841597V: 8.829003 X: 3.928906V: 8.927103 X: 4.017196V: 9.025203 X: 4.106467V: 9.123302 X: 4.196719V: 9.221402 X: 4.287952V: 9.319502 X: 4.380167V: 9.417602 X: 4.473361V: 9.515701 X: 4.567537V: 9.613801 X: 4.662694V: 9.711901 X: 4.758832V: 9.810000 X: 4.855951More Accurate Euler:V: 0.098100 X: 0.000491V: 0.196200 X: 0.001962V: 0.294300 X: 0.004415V: 0.392400 X: 0.007848V: 0.490500 X: 0.012262V: 0.588600 X: 0.017658V: 0.686700 X: 0.024035V: 0.784800 X: 0.031392V: 0.882900 X: 0.039731V: 0.981000 X: 0.049050V: 1.079100 X: 0.059351V: 1.177200 X: 0.070632V: 1.275300 X: 0.082895V: 1.373400 X: 0.096138V: 1.471500 X: 0.110362V: 1.569600 X: 0.125568V: 1.667700 X: 0.141755V: 1.765800 X: 0.158922V: 1.863900 X: 0.177070V: 1.962000 X: 0.196200V: 2.060100 X: 0.216311V: 2.158200 X: 0.237402V: 2.256299 X: 0.259474V: 2.354399 X: 0.282528V: 2.452499 X: 0.306562V: 2.550599 X: 0.331578V: 2.648699 X: 0.357574V: 2.746799 X: 0.384552V: 2.844899 X: 0.412510V: 2.942999 X: 0.441450V: 3.041099 X: 0.471370V: 3.139199 X: 0.502272V: 3.237299 X: 0.534154V: 3.335399 X: 0.567018V: 3.433499 X: 0.600862V: 3.531599 X: 0.635688V: 3.629699 X: 0.671494V: 3.727799 X: 0.708282V: 3.825899 X: 0.746050V: 3.923999 X: 0.784800V: 4.022099 X: 0.824530V: 4.120199 X: 0.865242V: 4.218299 X: 0.906934V: 4.316399 X: 0.949608V: 4.414499 X: 0.993262V: 4.512599 X: 1.037898V: 4.610700 X: 1.083514V: 4.708800 X: 1.130112V: 4.806900 X: 1.177690V: 4.905000 X: 1.226250V: 5.003100 X: 1.275790V: 5.101201 X: 1.326312V: 5.199301 X: 1.377814V: 5.297401 X: 1.430298V: 5.395501 X: 1.483762V: 5.493601 X: 1.538208V: 5.591702 X: 1.593634V: 5.689802 X: 1.650042V: 5.787902 X: 1.707430V: 5.886002 X: 1.765800V: 5.984102 X: 1.825150V: 6.082202 X: 1.885482V: 6.180303 X: 1.946794V: 6.278403 X: 2.009088V: 6.376503 X: 2.072362V: 6.474603 X: 2.136618V: 6.572703 X: 2.201854V: 6.670804 X: 2.268072V: 6.768904 X: 2.335270V: 6.867004 X: 2.403450V: 6.965104 X: 2.472610V: 7.063204 X: 2.542752V: 7.161304 X: 2.613875V: 7.259405 X: 2.685978V: 7.357505 X: 2.759063V: 7.455605 X: 2.833128V: 7.553705 X: 2.908175V: 7.651805 X: 2.984202V: 7.749906 X: 3.061211V: 7.848006 X: 3.139200V: 7.946106 X: 3.218171V: 8.044206 X: 3.298123V: 8.142305 X: 3.379055V: 8.240405 X: 3.460969V: 8.338505 X: 3.543863V: 8.436604 X: 3.627739V: 8.534704 X: 3.712595V: 8.632804 X: 3.798433V: 8.730904 X: 3.885252V: 8.829003 X: 3.973051V: 8.927103 X: 4.061831V: 9.025203 X: 4.151593V: 9.123302 X: 4.242336V: 9.221402 X: 4.334059V: 9.319502 X: 4.426764V: 9.417602 X: 4.520449V: 9.515701 X: 4.615116V: 9.613801 X: 4.710763V: 9.711901 X: 4.807392V: 9.810000 X: 4.905001NSV: Newton-Stormer-Verlet ('Semi-implicit' Euler):V: 0.098100 X: 0.000981V: 0.196200 X: 0.002943V: 0.294300 X: 0.005886V: 0.392400 X: 0.009810V: 0.490500 X: 0.014715V: 0.588600 X: 0.020601V: 0.686700 X: 0.027468V: 0.784800 X: 0.035316V: 0.882900 X: 0.044145V: 0.981000 X: 0.053955V: 1.079100 X: 0.064746V: 1.177200 X: 0.076518V: 1.275300 X: 0.089271V: 1.373400 X: 0.103005V: 1.471500 X: 0.117720V: 1.569600 X: 0.133416V: 1.667700 X: 0.150093V: 1.765800 X: 0.167751V: 1.863900 X: 0.186390V: 1.962000 X: 0.206010V: 2.060100 X: 0.226611V: 2.158200 X: 0.248193V: 2.256299 X: 0.270756V: 2.354399 X: 0.294300V: 2.452499 X: 0.318825V: 2.550599 X: 0.344331V: 2.648699 X: 0.370818V: 2.746799 X: 0.398286V: 2.844899 X: 0.426735V: 2.942999 X: 0.456165V: 3.041099 X: 0.486576V: 3.139199 X: 0.517968V: 3.237299 X: 0.550341V: 3.335399 X: 0.583695V: 3.433499 X: 0.618030V: 3.531599 X: 0.653346V: 3.629699 X: 0.689643V: 3.727799 X: 0.726921V: 3.825899 X: 0.765180V: 3.923999 X: 0.804420V: 4.022099 X: 0.844641V: 4.120199 X: 0.885843V: 4.218299 X: 0.928026V: 4.316399 X: 0.971190V: 4.414499 X: 1.015335V: 4.512599 X: 1.060461V: 4.610700 X: 1.106568V: 4.708800 X: 1.153656V: 4.806900 X: 1.201725V: 4.905000 X: 1.250775V: 5.003100 X: 1.300806V: 5.101201 X: 1.351818V: 5.199301 X: 1.403811V: 5.297401 X: 1.456785V: 5.395501 X: 1.510740V: 5.493601 X: 1.565676V: 5.591702 X: 1.621593V: 5.689802 X: 1.678491V: 5.787902 X: 1.736370V: 5.886002 X: 1.795230V: 5.984102 X: 1.855071V: 6.082202 X: 1.915893V: 6.180303 X: 1.977696V: 6.278403 X: 2.040480V: 6.376503 X: 2.104245V: 6.474603 X: 2.168991V: 6.572703 X: 2.234718V: 6.670804 X: 2.301426V: 6.768904 X: 2.369115V: 6.867004 X: 2.437785V: 6.965104 X: 2.507436V: 7.063204 X: 2.578068V: 7.161304 X: 2.649681V: 7.259405 X: 2.722275V: 7.357505 X: 2.795850V: 7.455605 X: 2.870406V: 7.553705 X: 2.945943V: 7.651805 X: 3.022461V: 7.749906 X: 3.099960V: 7.848006 X: 3.178441V: 7.946106 X: 3.257902V: 8.044206 X: 3.338344V: 8.142305 X: 3.419767V: 8.240405 X: 3.502171V: 8.338505 X: 3.585556V: 8.436604 X: 3.669922V: 8.534704 X: 3.755269V: 8.632804 X: 3.841597V: 8.730904 X: 3.928906V: 8.829003 X: 4.017196V: 8.927103 X: 4.106467V: 9.025203 X: 4.196719V: 9.123302 X: 4.287952V: 9.221402 X: 4.380166V: 9.319502 X: 4.473362V: 9.417602 X: 4.567538V: 9.515701 X: 4.662694V: 9.613801 X: 4.758832V: 9.711901 X: 4.855951V: 9.810000 X: 4.954051Velocity-less Verlet:V: 0.098100 X: 0.000981V: 0.196200 X: 0.002943V: 0.294300 X: 0.005886V: 0.392400 X: 0.009810V: 0.490500 X: 0.014715V: 0.588600 X: 0.020601V: 0.686700 X: 0.027468V: 0.784800 X: 0.035316V: 0.882900 X: 0.044145V: 0.981000 X: 0.053955V: 1.079100 X: 0.064746V: 1.177200 X: 0.076518V: 1.275300 X: 0.089271V: 1.373400 X: 0.103005V: 1.471501 X: 0.117720V: 1.569601 X: 0.133416V: 1.667700 X: 0.150093V: 1.765801 X: 0.167751V: 1.863901 X: 0.186390V: 1.962001 X: 0.206010V: 2.060102 X: 0.226611V: 2.158202 X: 0.248193V: 2.256302 X: 0.270756V: 2.354403 X: 0.294300V: 2.452501 X: 0.318825V: 2.550602 X: 0.344331V: 2.648702 X: 0.370818V: 2.746802 X: 0.398286V: 2.844903 X: 0.426735V: 2.943003 X: 0.456165V: 3.041103 X: 0.486576V: 3.139203 X: 0.517968V: 3.237304 X: 0.550341V: 3.335401 X: 0.583695V: 3.433498 X: 0.618030V: 3.531596 X: 0.653346V: 3.629693 X: 0.689643V: 3.727790 X: 0.726921V: 3.825888 X: 0.765180V: 3.923985 X: 0.804420V: 4.022082 X: 0.844641V: 4.120180 X: 0.885842V: 4.218277 X: 0.928025V: 4.316374 X: 0.971189V: 4.414472 X: 1.015333V: 4.512569 X: 1.060459V: 4.610672 X: 1.106566V: 4.708770 X: 1.153654V: 4.806867 X: 1.201722V: 4.904964 X: 1.250772V: 5.003062 X: 1.300802V: 5.101159 X: 1.351814V: 5.199256 X: 1.403806V: 5.297354 X: 1.456780V: 5.395451 X: 1.510734V: 5.493548 X: 1.565670V: 5.591646 X: 1.621586V: 5.689743 X: 1.678484V: 5.787840 X: 1.736362V: 5.885938 X: 1.795221V: 5.984035 X: 1.855062V: 6.082132 X: 1.915883V: 6.180229 X: 1.977685V: 6.278327 X: 2.040469V: 6.376436 X: 2.104233V: 6.474533 X: 2.168978V: 6.572643 X: 2.234705V: 6.670752 X: 2.301412V: 6.768861 X: 2.369101V: 6.866970 X: 2.437771V: 6.965080 X: 2.507422V: 7.063189 X: 2.578054V: 7.161298 X: 2.649667V: 7.259407 X: 2.722261V: 7.357517 X: 2.795836V: 7.455626 X: 2.870393V: 7.553735 X: 2.945930V: 7.651844 X: 3.022449V: 7.749954 X: 3.099948V: 7.848063 X: 3.178429V: 7.946172 X: 3.257891V: 8.044281 X: 3.338334V: 8.142391 X: 3.419758V: 8.240500 X: 3.502163V: 8.338609 X: 3.585549V: 8.436718 X: 3.669916V: 8.534828 X: 3.755265V: 8.632937 X: 3.841594V: 8.731046 X: 3.928905V: 8.829155 X: 4.017196V: 8.927241 X: 4.106469V: 9.025350 X: 4.196722V: 9.123435 X: 4.287956V: 9.221521 X: 4.380171V: 9.319606 X: 4.473367V: 9.417692 X: 4.567544V: 9.515777 X: 4.662702V: 9.613862 X: 4.758840V: 9.711948 X: 4.855960V: 9.810033 X: 4.954060My NSV variantV: 0.098100 X: 0.000981V: 0.196200 X: 0.002943V: 0.294300 X: 0.005886V: 0.392400 X: 0.009810V: 0.490500 X: 0.014715V: 0.588600 X: 0.020601V: 0.686700 X: 0.027468V: 0.784800 X: 0.035316V: 0.882900 X: 0.044145V: 0.981000 X: 0.053955V: 1.079100 X: 0.064746V: 1.177200 X: 0.076518V: 1.275300 X: 0.089271V: 1.373400 X: 0.103005V: 1.471500 X: 0.117720V: 1.569600 X: 0.133416V: 1.667700 X: 0.150093V: 1.765800 X: 0.167751V: 1.863900 X: 0.186390V: 1.962000 X: 0.206010V: 2.060100 X: 0.226611V: 2.158200 X: 0.248193V: 2.256300 X: 0.270756V: 2.354400 X: 0.294300V: 2.452500 X: 0.318825V: 2.550600 X: 0.344331V: 2.648699 X: 0.370818V: 2.746799 X: 0.398286V: 2.844899 X: 0.426735V: 2.942999 X: 0.456165V: 3.041099 X: 0.486576V: 3.139199 X: 0.517968V: 3.237299 X: 0.550341V: 3.335399 X: 0.583695V: 3.433499 X: 0.618030V: 3.531599 X: 0.653346V: 3.629699 X: 0.689643V: 3.727799 X: 0.726921V: 3.825899 X: 0.765180V: 3.923999 X: 0.804420V: 4.022099 X: 0.844641V: 4.120199 X: 0.885843V: 4.218299 X: 0.928026V: 4.316399 X: 0.971190V: 4.414498 X: 1.015335V: 4.512598 X: 1.060461V: 4.610698 X: 1.106568V: 4.708798 X: 1.153656V: 4.806898 X: 1.201725V: 4.904998 X: 1.250775V: 5.003098 X: 1.300806V: 5.101198 X: 1.351818V: 5.199298 X: 1.403811V: 5.297398 X: 1.456785V: 5.395498 X: 1.510740V: 5.493598 X: 1.565676V: 5.591698 X: 1.621593V: 5.689798 X: 1.678491V: 5.787897 X: 1.736370V: 5.885997 X: 1.795230V: 5.984097 X: 1.855071V: 6.082197 X: 1.915892V: 6.180297 X: 1.977695V: 6.278397 X: 2.040479V: 6.376497 X: 2.104244V: 6.474598 X: 2.168990V: 6.572698 X: 2.234717V: 6.670798 X: 2.301425V: 6.768899 X: 2.369114V: 6.866999 X: 2.437784V: 6.965099 X: 2.507435V: 7.063200 X: 2.578067V: 7.161300 X: 2.649680V: 7.259400 X: 2.722274V: 7.357500 X: 2.795849V: 7.455601 X: 2.870405V: 7.553701 X: 2.945942V: 7.651801 X: 3.022460V: 7.749902 X: 3.099959V: 7.848002 X: 3.178439V: 7.946102 X: 3.257900V: 8.044203 X: 3.338342V: 8.142303 X: 3.419765V: 8.240403 X: 3.502170V: 8.338504 X: 3.585555V: 8.436604 X: 3.669921V: 8.534704 X: 3.755268V: 8.632804 X: 3.841596V: 8.730905 X: 3.928905V: 8.829005 X: 4.017195V: 8.927105 X: 4.106466V: 9.025206 X: 4.196718V: 9.123306 X: 4.287951V: 9.221406 X: 4.380165V: 9.319507 X: 4.473360V: 9.417607 X: 4.567536V: 9.515707 X: 4.662693V: 9.613807 X: 4.758832V: 9.711908 X: 4.855951V: 9.810008 X: 4.954050Velocity Verlet variantV: 0.098100 X: 0.000491V: 0.196200 X: 0.001962V: 0.294300 X: 0.004415V: 0.392400 X: 0.007848V: 0.490500 X: 0.012262V: 0.588600 X: 0.017658V: 0.686700 X: 0.024035V: 0.784800 X: 0.031392V: 0.882900 X: 0.039731V: 0.981000 X: 0.049050V: 1.079100 X: 0.059351V: 1.177200 X: 0.070632V: 1.275300 X: 0.082895V: 1.373400 X: 0.096138V: 1.471500 X: 0.110362V: 1.569600 X: 0.125568V: 1.667700 X: 0.141755V: 1.765800 X: 0.158922V: 1.863900 X: 0.177070V: 1.962000 X: 0.196200V: 2.060100 X: 0.216311V: 2.158200 X: 0.237402V: 2.256299 X: 0.259474V: 2.354399 X: 0.282528V: 2.452499 X: 0.306562V: 2.550599 X: 0.331578V: 2.648699 X: 0.357574V: 2.746799 X: 0.384552V: 2.844899 X: 0.412510V: 2.942999 X: 0.441450V: 3.041099 X: 0.471370V: 3.139199 X: 0.502272V: 3.237299 X: 0.534154V: 3.335399 X: 0.567018V: 3.433499 X: 0.600862V: 3.531599 X: 0.635688V: 3.629699 X: 0.671494V: 3.727799 X: 0.708282V: 3.825899 X: 0.746050V: 3.923999 X: 0.784800V: 4.022099 X: 0.824530V: 4.120199 X: 0.865242V: 4.218299 X: 0.906934V: 4.316399 X: 0.949608V: 4.414499 X: 0.993262V: 4.512599 X: 1.037898V: 4.610700 X: 1.083514V: 4.708800 X: 1.130112V: 4.806900 X: 1.177690V: 4.905000 X: 1.226250V: 5.003100 X: 1.275790V: 5.101201 X: 1.326312V: 5.199301 X: 1.377814V: 5.297401 X: 1.430298V: 5.395501 X: 1.483762V: 5.493601 X: 1.538208V: 5.591702 X: 1.593634V: 5.689802 X: 1.650042V: 5.787902 X: 1.707430V: 5.886002 X: 1.765800V: 5.984102 X: 1.825150V: 6.082202 X: 1.885482V: 6.180303 X: 1.946794V: 6.278403 X: 2.009088V: 6.376503 X: 2.072362V: 6.474603 X: 2.136618V: 6.572703 X: 2.201854V: 6.670804 X: 2.268072V: 6.768904 X: 2.335270V: 6.867004 X: 2.403450V: 6.965104 X: 2.472610V: 7.063204 X: 2.542752V: 7.161304 X: 2.613875V: 7.259405 X: 2.685978V: 7.357505 X: 2.759063V: 7.455605 X: 2.833128V: 7.553705 X: 2.908175V: 7.651805 X: 2.984202V: 7.749906 X: 3.061211V: 7.848006 X: 3.139200V: 7.946106 X: 3.218171V: 8.044206 X: 3.298123V: 8.142305 X: 3.379055V: 8.240405 X: 3.460969V: 8.338505 X: 3.543863V: 8.436604 X: 3.627739V: 8.534704 X: 3.712595V: 8.632804 X: 3.798433V: 8.730904 X: 3.885252V: 8.829003 X: 3.973051V: 8.927103 X: 4.061831V: 9.025203 X: 4.151593V: 9.123302 X: 4.242336V: 9.221402 X: 4.334059V: 9.319502 X: 4.426764V: 9.417602 X: 4.520449V: 9.515701 X: 4.615116V: 9.613801 X: 4.710763V: 9.711901 X: 4.807392V: 9.810000 X: 4.905001

By using the old accel, this Velocity Verlet (VV) variant (other variants exist) only requires one derivative evaluation per step. In this simple example, VV has the same accuracy as "More Accurate Euler". The second most accurate is my NSV variant. However, all of these variants are more than accurate enough for most game applications.

The next step would be to show a demo of a stiff enough spring-mass example that shows when each variant goes unstable at a particular delta time step size. Higher order methods can show their value by allowing for a larger time step, saving iteration cycles, which can offset the extra memory/cache/cpu overhead. Comparing these methods to RK2/RK4 (2 and 4 derivative evaluations per step) would also be interesting for developers.

My staggered/averaging tests (similar to the Velocity Verlet variant above) were inspired by the nature of RK2, using the same trick (by staggering evaluations and using a previous value to prevent requiring the extra derivative computation).

Note: when a rigid impulse (such as with a collision) must be applied, the higher order methods must use a different method for integration (in this example, the oldAccel would end up making the net impulse incorrect for the collision event). This concept is also illustrated in Baraff's papers (which are generalized for most types of integrators), where a 'discontinuous' notification must be used.

Numsgil

Feb 12, 2006 at 12:48 PM

Quote:
Original post by jjd
Quote:
Original post by Numsgil
They are mathematically equivelant according to David Kofke of the Department of Chemical Engineering at SUNY Buffalo. See my "source 3" above (which are his lecture notes apparently). It may depend on exactly what you mean by "mathematically equivelant" of course. I'm assuming Kofke meant they are algebraic transformations of one another.

Do you have any sources for your assertion that vanilla verlet has better positional accuracy than velocity verlet?

I used to be a member of the faculty in the dept of applied mathematics and the university of colorado in boulder. I'm not too interested in the claims of a chemical engineer on the mathematical equivalence of two numerical integrators.

...

Thanks, it's been difficult for me to find error analysis for the different algorithms.

You'll have to forgive me, it's been a few years since my Numerical Analysis classes. Am I right in the following summary?

vanilla verlet:
position error: O(h^4)
velocity error: O(h^2)

velocity verlet:
position error: O(h^3)
velocity error: O(h^3)

If so, then it would seem to me at least that velocity verlet is indeed preferable to vanilla verlet for games in general (since velocities are at least as important as positions).

Quote:

I suspect that the reason only a few methods are commonly used in games reflects the history and requirements that are particular to games:
(1) kickass physics engines are a relatively recent development
(2) physical accuracy is not so important as it is in the physical sciences
(3) the realtime demand on games is a severe constraint, although one that will probably lessen as computers become more powerful.

Games require fast, stable methods. Often the amount of memory used is a constraint too. This limits the methods that people are willing to explore. Implicit methods would be great, if only they didn't come with a high cost to memory. Multi-step methods suffer a similar fate and the additional problem of start-up problems and that changing the size of the timestep is awkward. As such, explicit methods with few function evaluations are currently the best option for most games. I believe that will change in time, and I hope that more use can be made of implicit methods so that my ragdolls will behave more stably [smile]

Someone more knowledgable than I should write a short article that describes different families of integration methods for the purposes of games. I agree that multi-step methods probably aren't the best route for games. Perhaps suggest some future routes that can be explored once computers become more powerful.

Quote:

Oh, almost forgot. You wanted some sources "the art of molecular dynamics simulation" by rapaport; "numerical recipes" by press et al.; "solving ordinary differential equation I and II" by hairer and wanner. I also really like "computer methods for ordinary differential equations and differential-algebraic equations" by ascher and petzold.

Thanks I'll give them a look. Numerical Recipes is one of my favorite books, I find myself coming back to it time and time again. You can find it online at http://library.lanl.gov/numerical/bookcpdf.html

jjd

Feb 12, 2006 at 12:57 PM

Quote:
Original post by Numsgil

You'll have to forgive me, it's been a few years since my Numerical Analysis classes. Am I right in the following summary?

vanilla verlet:
position error: O(h^4)
velocity error: O(h^2)

velocity verlet:
position error: O(h^3)
velocity error: O(h^3)

If so, then it would seem to me at least that velocity verlet is indeed preferable to vanilla verlet for games in general (since velocities are at least as important as positions).

Slight correction.

vanilla verlet:
position error: O(h^4)
velocity error: O(h)

velocity verlet:
position error: O(h^3)
velocity error: O(h^2)

@John Schultz: Comparing the exact solution next to the numerical solutions would be a nice addition to your simulations.

Numsgil

Feb 12, 2006 at 12:59 PM

Quote:
Original post by John Schultz
*** Source Snippet Removed ***

Output:
*** Source Snippet Removed ***

...

The next step would be to show a demo of a stiff enough spring-mass example that shows when each variant goes unstable at a particular delta time step size. Higher order methods can show their value by allowing for a larger time step, saving iteration cycles, which can offset the extra memory/cache/cpu overhead. Comparing these methods to RK2/RK4 (2 and 4 derivative evaluations per step) would also be interesting for developers.

I agree such a demo would be of value.

Quote:

My staggered/averaging tests (similar to the Velocity Verlet variant above) were inspired by the nature of RK2, using the same trick (by staggering evaluations and using a previous value to prevent requiring the extra derivative computation).

Do you have an error term for position and velocity?

I would actually be interested in error term analysis of all the techniques you showed.

Quote:

Note: when a rigid impulse (such as with a collision) must be applied, the higher order methods must use a different method for integration (in this example, the oldAccel would end up making the net impulse incorrect for the collision event). This concept is also illustrated in Baraff's papers (which are generalized for most types of integrators), where a 'discontinuous' notification must be used.

I think I know where you're coming from, but I don't understand what you mean when you say "use a different method for integration".

Numsgil

Feb 12, 2006 at 1:04 PM

Quote:
Original post by jjd
Slight correction.

vanilla verlet:
position error: O(h^4)
velocity error: O(h)

I'm confused then. This article says it's O(h^2). Maybe you're using different methods for calculating velocity?

edit:

Yes, I see now, you are using different methods. The article uses time step t-1 and t+1, whereas you used t-1 and t

John Schultz

Feb 12, 2006 at 1:19 PM

Quote:
Original post by Numsgil
I agree such a demo would be of value.

Quote:
Original post by jjd
Comparing the exact solution next to the numerical solutions would be a nice addition to your simulations.

The source was posted for others to experiment with... Feel free to use it in any way you wish. [wink]

Quote:
Original post by Numsgil
Do you have an error term for position and velocity?

I'm primarily interested with stability and correct handling of total kinetic energy: the motion needs to look correct (given the time steps used, any error has not been an issue).

Quote:
Original post by Numsgil
I think I know where you're coming from, but I don't understand what you mean when you say "use a different method for integration".

What happens to the total kinetic energy when an object/particle hits a rigid surface and must instantaneously change direction in a perfectly elastic collision?

Numsgil

Feb 12, 2006 at 2:14 PM

Quote:
Original post by John Schultz
I'm primarily interested with stability and correct handling of total kinetic energy: the motion needs to look correct (given the time steps used, any error has not been an issue).

So that would be no? :D

Quote:
Original post by Numsgil
I think I know where you're coming from, but I don't understand what you mean when you say "use a different method for integration".

What happens to the total kinetic energy when an object/particle hits a rigid surface and must instantaneously change direction in a perfectly elastic collision?[/quote]

As far as I know, you'd just find the velocity parallel to the normal vector of the collision, then take the opposite sign. That's actually fairly easy to do. KE would remain constant.

Maybe you're talking about how you generally integrate before you apply constraints?

Topic Locked

This topic has been locked by a moderator. New replies are not allowed.