Verlet Integration Problem

Dec 5, 2008 at 5:14 AM

Hello, Im working on Verlet Integration. I got a problem during Integration. Thats the source:



void Particle::Integrate(real time) {

	real dt = 0.01;
	Vector3 temp = position;
	position += (position - prev_position) + acceleration * dt * dt;
	prev_position = temp;

}

I started testing my Verlet. I got Results:



/* Testing */

//Specs: 

position = Vector3(0, 50, 0)
acceleration = Vector3(0, -10, 0)
velocity = Vector3(0, 5, 0)
Mass = 10
Integration Timestep = 0.01

//Results of Posiion:
Seconds | Position (y-coorinate)|
--------+-----------------------|
0(start | 50                    |
of App) |                       |
--------+-----------------------|
1       | 99.999001             |
--------+-----------------------|
2       | 149.997001            |
--------+-----------------------|
3       | 199.994019            |
--------+-----------------------|
4       | 249.990021            |
--------+-----------------------|

If the Acceleration is -10 why the position needs to increase i didn't understand. Please Help! Thanks, Kasya

Dec 5, 2008 at 5:41 AM

Do you ever initialize prev_position?

You set your initial velocity to (0 5 0), so i guess prev_position should be position - velocity * dt ?

Hope it helps.
//Emil

Mod

198

Dec 5, 2008 at 10:06 AM

It is correct for the position to increase even when acceleration is negative. The reason is that the object (for your case) is initially moving with a positive velocity and velocity has a larger impact on the change of position than acceleration, since velocity is O(dt) in the equations of motion while acceleration is O(dt²) and for small time steps, dt > dt². The negative acceleration reduces the velocity every time step, but it will take several time steps before velocity goes negative. After velocity goes negative, then the position will start to decrease. So, just run your code for more time steps and you'll see position decrease.

Now, it might seem confusing to talk about velocity when your basic Verlet integration does not include velocity. But, your first two positions, position(t = 0.0) and position(t = dt) are based on that positive velocity. The Verlet integration that you have implemented implicitly represents velocity, since velocity is baked into the first two values of position.

(NOTE: I originally documented a correction to the integration, a correction that cannonicus later pointed out was incorrect. So, I have now removed the incorrect correction.)

[Edited by - grhodes_at_work on December 5, 2008 10:06:14 AM]

Author

Dec 5, 2008 at 10:07 AM

Hello,

I just changed the prev_position = temp - velocity * dt;

void Particle::Integrate(real time) {	real dt = 0.01;	Vector3 temp = position;	position += ((position*2.0) - prev_position) + acceleration * dt * dt;	prev_position = temp - velocity * dt;}

But i shocked here because the Results Are Absolutely same with: prev_position = temp - velocity * dt

What can i do?

[EDIT] I Added position*2.0

[grhodes_at_work note: the 2.0 multiplier is incorrect, as cannonicus points out below. I apologize for having asked you to make that change...it was my mistake to point this out.]

I tested the code without velocity too. But the results are same.

Thanks,
Kasya

[Edited by - grhodes_at_work on December 5, 2008 10:07:14 AM]

Mod

198

Dec 5, 2008 at 10:11 AM

Quote:
Original post by cannonicus
Do you ever initialize prev_position?

You set your initial velocity to (0 5 0), so i guess prev_position should be position - velocity * dt ?

Hope it helps.
//Emil

Kasya is using the basic Verlet integration, which represents velocity implicitly, e.g., velocity is not needed in the integration step. The derivation of the scheme cancels out the velocity terms.

Author

Dec 5, 2008 at 10:33 AM

I don't understand that Verlet Integration. I read all resources in the Internet but found nothing. I want to know if i use Verlet will there be any change if there is velocity?

I read the Wikipedia:

x(dt) = x(0) + v(0)*dt + a(0) * dt * dt / 2 + O(dt*dt*dt)

x(0) -> prev_position
v(0) -> velocity
a(0) -> acceleration
dt -> Time Step

Am i right about that? And what is O(dt^3)?

Mod

198

Dec 5, 2008 at 10:39 AM

Quote:
Original post by Kasya
Hello,

I just changed the prev_position = temp - velocity * dt;

*** Source Snippet Removed ***

But i shocked here because the Results Are Absolutely same with: prev_position = temp - velocity * dt

What can i do?

Thanks,
Kasya

Kasya, your original code, prev_position = temp, was correct. Change it back. Check my reply, as you have an error in your position integration.

I explained in my other post why it is possible for position to increase for a while until velocity (implicitly) becomes negative.

But, in your case, there is another problem. cannonicus identified it, but you did not address it correctly. You did not initialize prev_position to a meaningful value. You appear to have initialized prev_position to (0,0,0). So, when you begin your integration for the first time, you have:

prev_position = (0,0,0)
position = (0,50,0)

BEFORE you call Integrate. So, for your first integration, you get:

position += ((0,50,0) - (0,0,0)) + (0,-10,0) * .001;

Or, position += (0,50,0) + (0,-.1,0) = (0,99.9,0)

The issue here is...since you did not compute a prev_position that reflects your desired initial velocity of (0,5,0), you get a much higher initial velocity. Actually, for the time step you use, dt = 0.1, you have an initial velocity of (0,5000,0)!!! Rather large! It'll take many, many, many time steps before your velocity goes negative and the position starts decreasing!

So, here is the solution. Before you begin the simulation loop, compute prev_position based on your desired starting position and velocity, using a simple two point difference:

desired_velocity = (position - prev_position) / dt;

So,

prev_position = position - (desired_velocity * dt);

Or,

prev_position = (0,50,0) - 0.01*(0,5,0);

Or,

prev_position = (0,49.95,0);

So, then, if you START your simulation with prev_position having been set to (0,49.95,0) and position = (0,50,0), you will be working with values that properly represent your initial velocity of (0,5,0). Your simulation (with my correction as shown in the other post) should go fine. It will still take a few steps before position will start to reduce, but it will happen.

Hope this helps.

Dec 5, 2008 at 10:40 AM

Damn, you beat me to it ;)

Quote:
Original post by grhodes_at_work
Quote:
Original post by cannonicus
Do you ever initialize prev_position?

You set your initial velocity to (0 5 0), so i guess prev_position should be position - velocity * dt ?

Hope it helps.
//Emil

Kasya is using the basic Verlet integration, which represents velocity implicitly, e.g., velocity is not needed in the integration step. The derivation of the scheme cancels out the velocity terms.

Yes, but the previous position of the particle is. What i meant is that he initializes position to (0 50 0) but doesnt (as far as i see) initialize prev_position (leaving it to 0). Thus the particle will have initial velocity of (0; 50/dt; 0) and not (0 5 0), which gives the exact numbers posted in the first post.

Kasya:

I didnt mean you should do prev_position = position - velocity*dt every integration tick, only before you start your simulation.

So, basically what im saying:

/* Testing *///Specs: position = Vector3(0, 50, 0)acceleration = Vector3(0, -10, 0)velocity = Vector3(0, 5, 0)Mass = 10Integration Timestep = 0.01prev_position = position - velocity*Integration timestep /////////////////

You shouldnt have 2*position when you use += operator. You were right in your first post.

Also, about your first post: you do have an initial velocity of (0 5 0), so its natural for the particle to move upwards for some time before reaching its peak height. Try velocity (0 0 0) to debug. It shouldnt give you any positive motion in y-direction.

Quote:

I tested the code without velocity too. But the results are same.

This is just more proof that prev_position is you problem.

//Emil

Mod

198

Dec 5, 2008 at 10:46 AM

Quote:
Original post by Kasya
I don't understand that Verlet Integration. I read all resources in the Internet but found nothing. I want to know if i use Verlet will there be any change if there is velocity?

I read the Wikipedia:

x(dt) = x(0) + v(0)*dt + a(0) * dt * dt / 2 + O(dt*dt*dt)

x(0) -> prev_position
v(0) -> velocity
a(0) -> acceleration
dt -> Time Step

Am i right about that? And what is O(dt^3)?

Kasya,

Please go back to your original code. Your original code was mostly fine. The big problem was that you had not computed prev_position and so your simulation had a VERY large initial velocity represented in the positions. I think this has been covered adequately here by now!

The Wikipedia entry you are looking at shows a varient called "Velocity Verlet." It is a two-step version of Verlet integration. It does have the benefit of not having to pre-compute prev_position before starting the simulation.

If you compute a proper value for prev_position before beginning (see my 2nd reply to you), then your original code will work.

The O(dt³) is an annotation that documents the truncation error represented by the equation. It is not part of the equation that you need to code. All it is saying is that the truncation error present in the formula is the same order of magnitude, or proportional to, dt³. the "O" notation, often called "big-O" notation, does indicate "order of magnitude." Another way of stating O(dt³) is that the quation is "3rd order accurate in time."

Mod

198

Dec 5, 2008 at 10:49 AM

Quote:
Original post by cannonicus
Damn, you beat me to it ;)

Actually, didn't you beat me to it in your first reply? Kasya just didn't understand where to set prev_position to position - velocity * dt, and I clarified that! And, then, you clarified that also! Ha Ha! I'm glad to see you here!

Mod

198

Dec 5, 2008 at 10:51 AM

Quote:
Original post by cannonicus
You shouldnt have 2*position when you use += operator. You were right in your first post.

Crap, you're right! Oops...wow, I made a mistake! Thank you for correcting it. I'll make a note in my other posts.

Author

Dec 5, 2008 at 11:07 AM

I DID THAT! :) Thank you very very much! It goes with very very small rate. Does it need to be like that? Or not:

sec pos0   50.0000001   49.9990012   49.9970023   49.9940034   49.9900055   49.9850086   49.979012

Dec 5, 2008 at 11:39 AM

Well, its only 0.01 sec between each line, so its not surprising the value is changing so slowly. (Assuming you now have velocity = (0 0 0), else something is still rotten)

Congrats by the way and have fun with your verlet integrator :)

//Emil

Author

Dec 5, 2008 at 1:34 PM

Now im going to write the calculation of timestep and Time Corrected Verlet Integration. Thank you very very much! :)

Good Luck!
Kasya

Author

Dec 6, 2008 at 5:04 AM

Hello,
I have finished the my Verlet Integration System:

My Verlet System

Please look at it and play with acceleration and tell me if that is correct please because i think there is a bug here.

The bug im thinking of:
When i change acceleration does it need to do before simulation thing:
prev_position = position - velocity*Integration timestep, or not?

Note: Timestep is 0.01 here

Thanks,
Kasya

h4tt3n

Dec 6, 2008 at 6:57 AM

Hi Kasya,

It's a bit hard to comment on the program when you don't know what it's supposed to do. If the downwards movement is caused by a continuous gravitational pull, then something is probably wrong, since the object doesn't seem to accelerate, just move downwards at a constant speed. If on the other hand you just apply one "knock" at the object on program startup, then the constant velocity is right.

Please add a few details, or maybe even the source code.

Cheers,
Mike

Author

Dec 6, 2008 at 7:16 AM

Thats the Integration

void Particle::startSimulation(real dt) {	prev_position = position - velocity * dt;}void Particle::Integrate(real time) {	real dt = 0.01;	Vector3 temp = position;	position += (position - prev_position) + acceleration * dt * dt;	prev_position = temp;}

And Thats the Test Code:

void InitOpenGL() {//OpenGL Functions like glEnable ...	particle.SetMass(10);	particle.SetPosition(Vector3(0,3,0));	particle.SetVelocity(Vector3(0,2,0));		Physics::ApplyForce(&particle, "Gravity"); //Just add -10 to y acceleration	particle.startSimulation(0.01);}void RenderOpenGL() {//Render Functions like LoadIdentity...	glPointSize(particle.getMass()/0.4); //For making object real :)	glBegin(GL_POINTS);	glColor3f(1,0,1);	glVertex3f(2,particle.getPosition().y,0);	glEnd();	if(simulate) {		particle.Integrate(0);	}      //The Texts Function (Mouse Position, Information, Simulation Button ...)}

Author

Dec 6, 2008 at 7:24 AM

I just Added the Position Difference between current and last position. It is decreasing. If it is accelerating downward it needs to decrease right?

*EDIT* Just Forgotten to add Absolute Value for difference. Its always increasing now. LoL.

Now the problem here is when i make acceleration positive, there comes a bug i guess so. It goes downward a little bit and then goes upward. Why? (Change the acceleration from -10 to 10 and see what happens)

Thanks,
Kasya

[Edited by - Kasya on December 6, 2008 10:24:35 AM]

Author

Dec 6, 2008 at 1:39 PM

So What you think?