## Newton Was Wrong ?

Please don't post Bullet support questions here, use the above forums instead.
Posts: 7
Joined: Fri Aug 28, 2009 8:54 am

### Newton Was Wrong ?

Hi !

I am trying to make a particle simulator using Newtons law of motion in C++.
I use a Vector class and Atom class which is represented by r - position, v - velocity and a - acceleration (im not using BulletPhysicsLibrary right now).

In Atom class I have Update function :

Code: Select all

Update(float dt)
{
r += (v * dt) + (a * 0.5f * dt * dt); // dr = v*dt + 0.5*a*dt^2
v += (a * dt); // dv/dt = a
}

I set dt for about 0.005f or 0.001f.

Now I try to simulate a normal spring using Hooke's Law.
In main program I simply write something like :

Code: Select all

my_atom->r.y = 2; // starting position, so it will get accelerated
while(...)
{
my_atom->a.y = -2*my_atom->r.y; // a = -(k/m)*x .. (-k/m)==-2
my_atom->Update(dt);
}
But the problem is that every cicle the atom getting his Amplitude increase ... If I would draw it, then you would see as the atom just keep getting enery instead of keeping a consistent one !
I think it happends because my dt isnt small enough (should be 0.00000000001 I guess but its not real to use) ?

Can someone help me solve this problem ?

bone
Posts: 235
Joined: Tue Feb 20, 2007 4:56 pm

### Re: Newton Was Wrong ?

Welcome to the wide world of numerical integration. Your Update() function is integrating time forward, specifically with a method often called "parabolic" integration. This method is perfect if and only if the acceleration is constant throughout the timestep. Otherwise it is somewhat unstable. Since a real spring would be changing the force and hence acceleration in the middle of the timestep, this method is not particularly suitable.

Around these parts, we recommend using the symplectic Euler integration method (a.k.a. semi-implicit and possibly some other names):

v += ( a * dt )
r += ( v * dt )

Notice that we are using the *updated* (next timestep's) velocity for the position integration. This is conditionally stable, and more stable than the parabolic method. If the dt is small enough, you will be able to simulate your spring with pretty good accuracy for quite awhile. Hope this helps.

Posts: 7
Joined: Fri Aug 28, 2009 8:54 am

### Re: Newton Was Wrong ?

Wow, it is surprisingly accurate !
In the first few runs I even got an integer instead of .9998 values.

I also wiki-ed some strange words you used (AKA "numerical integration) and I didn't thought that subject is so wide, I acctualy didnt know thers a subject for my problem at all !

I'll sure have reading matirials now, thanks for the help !

h4tt3n
Posts: 23
Joined: Wed Mar 12, 2008 9:08 am

### Re: Newton Was Wrong ?

An even better integration method is the Velocity Verlet, which is an order of magnitude more accurate than symplectic Euler.

r += v * dt + 0.5 * a * dt * dt
v += 0.5 * a * dt
a = f / m
v += 0.5 * a * dt

cheers,
Mike
Last edited by h4tt3n on Fri Sep 04, 2009 12:00 pm, edited 1 time in total.

bone
Posts: 235
Joined: Tue Feb 20, 2007 4:56 pm

### Re: Newton Was Wrong ?

Even better? In what way? If those equations are correct for the Velocity Verlet, then it has the *exact* same solution for a mass/spring/damper system as the much simpler Symplectic Euler method. The only meaningful difference is that Verlet is wasting CPU time because the equations are more complex to accomplish the same thing.

There is another difference that I consider negligible: the velocity for the Symplectic method is equivalent to the midpoint velocity rather than the final velocity of the Verlet method. This difference *could* be considered important if you are measuring the energy of the system at the end of the timestep, as it might appear to be wobbling even though it really isn't. In any case, one could easily move the bottom two equations of Verlet to the top of the loop, though, and then it becomes perfectly equivalent.

P.S. Unless I'm mistaken, I think the Velocity Verlet method has a bunch of other names including Newton-Stormer-Verlet and possibly Leapfrog as well.

h4tt3n
Posts: 23
Joined: Wed Mar 12, 2008 9:08 am

### Re: Newton Was Wrong ?

Admittedly the difference between the two is very small, but it is there. If you add an energy check to your simulation, you'll see that VV is more accurate than SE. VV is very similar to the leapfrog method, but not identical (vv is slightly more accurate).

Posts: 7
Joined: Fri Aug 28, 2009 8:54 am

### Re: Newton Was Wrong ?

I don't really understand how people got to those equations, but after 3110 iteration with VV i got an INTEGER value when it acctualy needed to be an integer! Thats outstandingly accurate and I didnt even use a rounding-error algorithem.

EV gave me an integer too but at the 3109 AND 3110.. I dont really know whats the difference but i'll test it more when I use it in the acctual graphic simulation (all the above were used in command line to check).

And just to point out, Newton laws gives me an error of .198 at this stage

EDIT:
I did now a test to check the Energy of the system and you are right, Euler is abit bluffing (+-0.5%) while Verlet keeping it energy at 24.999+-0.001% or something like that (I guess on dt and other parameters as well).

bone
Posts: 235
Joined: Tue Feb 20, 2007 4:56 pm

### Re: Newton Was Wrong ?

I'm not convinced yet that VV is any different than Symplectic Euler. I've got an Excel sheet here that gives the *exact* same results for both methods no matter what the inputs are (again, comparing the VV's midpoint velocity to the SE's final velocity). Maybe there's an error in the Excel sheet, but this is pretty dead simple so I'm not sure how.

Looking at it from a different perspective, rotating VV's loop so the force computation is at the top, and labeling the velocity and position changes:

Starting with v0 and r0 (starting velocity and position, respectively):
a = f / m
v.5 = v0 + 0.5 * a * dt // v.5 is really the final velocity in the original order
r1 = r0 + v.5 * dt + 0.5 * a * dt * dt
v1 = v.5 + 0.5 * a * dt // v1 is really the midpoint velocity in the original order

Substituting v.5 into the r1 and v1 equations:

r1 = r0 + ( v0 + 0.5 * a * dt ) * dt + 0.5 * a * dt * dt
v1 = ( v0 + 0.5 * a * dt ) + 0.5 * a * dt

Now simplify and rearrange:

r1 = r0 + ( v0 + a * dt ) * dt // look at the term multiplied by dt,
v1 = v0 + a * dt // it's the same as the right-side here

So ...

r1 = r0 + v1 * dt
v1 = v0 + a * dt

These are the equations for Symplectic Euler. It's the same. I fail to see how re-arranging it changes the results, and yet you're still claiming that it's an order of magnitude more accurate in your original post. It is not.

Yes, the energy measured at the final point wobbles with SE, but over the long term this is completely inconsequential - over a million timesteps, it doesn't gain or lose any more energy than VV. In fact you can simply adjust your energy measurement for SE to report the same thing as VV if you want (I believe by taking the average of v0 and v1).

Posts: 7
Joined: Fri Aug 28, 2009 8:54 am

### Re: Newton Was Wrong ?

I think it matters to put the force calculation where it is, since at v(t+dt/2) u calculate velocity with acceleration a1, and at v(t+dt/2+dt/2) u calculate the velocity with a new acceleration a2, and it changed because the position of the body changed (I simulate a spring making acceleration obay m*a=F=-k*x).

those are my functions :

inline void EulerMethod(float dt){
v += (a * dt); // dv/dt = a
r += (v * dt);// + (a * 0.5f * dt * dt); // dr = v*dt + 0.5*a*dt^2
}
inline void VerletVelocityMethod(float dt){
r += (v * dt) + (a * 0.5f * dt * dt); //.. r(t+dt)
v += a * 0.5f * dt; //.. v(t+dt/2)
a = pState(this,r,v,a)/mass; // pState points to function to get the force on object
v += a * 0.5f * dt; //.. v(t+dt)
}

VV show me the energy of the system being between 25 (integer!) and 24.99996
Euler shows energy between 24.91 .. 25.08.
so ye if make average of Euler energy it can be pretty precice, although BOTH methods are pretty precice just VV takes more time/memory since I need to call a function pState to calculate the force all over again (while hoping the force calculation is short enough).

h4tt3n
Posts: 23
Joined: Wed Mar 12, 2008 9:08 am

### Re: Newton Was Wrong ?

Please note that both symplectic Euler and velocity verlet requires just one force update per simulation loop! In your euler example, you also need to call the Pstate() function, it just doesn't show from the posted example. In other words, vv is only very slightly slower than symplectic euler.

Posts: 7
Joined: Fri Aug 28, 2009 8:54 am

### Re: Newton Was Wrong ?

*Aw my messege was deleted ... what is this 'save' option anyway ?*

To topic, I did calculate the force in the start of each loop I just didnt put it on the functions I wrote above.
About calculating the force twice in Verlet, if you may look here:
http://www.compsoc.man.ac.uk/~lucky/Dem ... erlet.html
at the buttom of the page you can see f(n) and f(n+1) which states its not the same force!

Plus at wiki they give algorithem which I followed:
http://en.wikipedia.org/wiki/Verlet_int ... ity_Verlet

So from what I see I do need to calculate the force twice which is also the 'bad' thing in Verlet.

h4tt3n
Posts: 23
Joined: Wed Mar 12, 2008 9:08 am

### Re: Newton Was Wrong ?

Indeed f(n) and f(n+1) are two different forces: f(n) is the force you calculated in the previous simulation loop and f(n+1) is the force you calculate during the current loop. Still, all you need is one force update per loop!

This also shows from the wikipedia article, where the new acceleration (which is just force / mass) is derived once per loop (in step 3 where it says "derive a").

Besides, I've done all sorts of experiments with this integration algorithm, including adding a second force update, but it doesn't improve accuracy or stability.

IMHO velocity verlet is the integration algorithm which gives you the best and cheapest allround result. Popular higher order algorithms like the Runge-Kutta 4th order requires four force updates per loop, and it is not even symplectic (ie. it does not preserve energy). Swapping the RK4 out with four VV itrations at timestep dt/4 gives you better force preservation over time and better stability for slightly less cpu power.

Cheers,
Mike

Posts: 7
Joined: Fri Aug 28, 2009 8:54 am

### Re: Newton Was Wrong ?

Kinda confusing to be honest ! I thought I got it right ..
I did some adjustments to my VV algorithem according to what you said, now energy goes between 25 to 25.00003 which is still better then Euler.

so I changed from :

Code: Select all

a = pState(this,r,v,a);
r += (v * dt) + (a * 0.5f * dt * dt); //.. r(t+dt)
v += (a*mass + pState(this,r,v,a)) * (0.5f * dt /mass); //.. v(t+2*dt/2)

to

Code: Select all

r += (v * dt) + (prev_a * 0.5f * dt * dt); //.. r(t+dt)
v += prev_a*0.5f*dt;
a = pState(this,r,v,prev_a) / mass;
prev_a = a;
v += a * 0.5f * dt; //.. v(t+dt)

so the trick is that I just calculate the force in a difference place in the algorithem (didnt thought about it before).

One problem is that when i do speed-check (using milliseconds clock()) I get about the same running time for both, I think its because the latter use more lines of codes and its harder to 'simplify' it like the first one which is just 2 lines and no division of vector by scalar.
EDIT: I made some adjustments, VV now is 300 ms (while previously was 360 ms) and euler stays at about 215 ms.

Thanks for clarifying it to me

h4tt3n
Posts: 23
Joined: Wed Mar 12, 2008 9:08 am

### Re: Newton Was Wrong ?

You're welcome. Try this:

Code: Select all

float halfdt = 0.5f * dt;
float invmass = 1.0f / mass;
...
do
...
r += v * dt + a * halfdt * dt;
v += a * halfdt;
a = pState(this, r, v, a) * invmass;
v += a * halfdt;
...
loop
This should be even faster (but since I'm not sure how your pState func works the syntax might not be entirely right).