# Verlet Integration Problem

This topic is 3654 days old which is more than the 365 day threshold we allow for new replies. Please post a new topic.

## Recommended Posts

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

##### Share on other sites
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

##### Share on other sites
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(dt2) and for small time steps, dt > dt2. 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]

##### Share on other sites
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?

[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]

##### Share on other sites
Quote:
 Original post by cannonicusDo 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.

##### Share on other sites
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?

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)?

##### Share on other sites
Quote:
 Original post by KasyaHello,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 * dtWhat 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.

##### Share on other sites
Damn, you beat me to it ;)

Quote:
Original post by grhodes_at_work
Quote:
 Original post by cannonicusDo 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

##### Share on other sites
Quote:
 Original post by KasyaI 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_positionv(0) -> velocitya(0) -> accelerationdt -> Time StepAm 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(dt3) 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, dt3. the "O" notation, often called "big-O" notation, does indicate "order of magnitude." Another way of stating O(dt3) is that the quation is "3rd order accurate in time."

##### Share on other sites
Quote:
 Original post by cannonicusDamn, 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!

1. 1
2. 2
3. 3
Rutin
15
4. 4
khawk
13
5. 5

• 9
• 9
• 11
• 11
• 23
• ### Forum Statistics

• Total Topics
633665
• Total Posts
3013239
×