Sign in to follow this  

Verlet Integration Problem

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

If you intended to correct an error in the post then please contact us.

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 this post


Link to post
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 this post


Link to post
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?

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

Share this post


Link to post
Share on other sites
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.

Share this post


Link to post
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?

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

Share this post


Link to post
Share on other sites
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.

Share this post


Link to post
Share on other sites
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 = 10
Integration Timestep = 0.01

prev_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 this post


Link to post
Share on other sites
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(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 this post


Link to post
Share on other sites
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!

Share this post


Link to post
Share on other sites
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.

Share this post


Link to post
Share on other sites
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 pos
0 50.000000
1 49.999001
2 49.997002
3 49.994003
4 49.990005
5 49.985008
6 49.979012



?

Share this post


Link to post
Share on other sites
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

Share this post


Link to post
Share on other sites
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

Share this post


Link to post
Share on other sites
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

Share this post


Link to post
Share on other sites
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 ...)

}


Share this post


Link to post
Share on other sites
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]

Share this post


Link to post
Share on other sites
Quote:

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


No, dont do that. You only need to do that before you start the simulation. When the simulation is running prev_position will always be the position of the particle in the last frame, you can apply any acceleration each simulation tick.

Honestly i cant see anything wrong with the simulation i downloaded. Seemed to work allright. But if you want to be sure you can always calculate the particles total energy each frame. As long as you dont change the acceleration the total energy should be constant.

Total energy = Kinetic energy + potential energy

Kinetic energy = 1/2*m*v^2
potential energy = -m*y-acceleration*y-coordinate

Quote:

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)



I couldnt reproduce this.

Share this post


Link to post
Share on other sites
i don't understand what you wrote in Potensial Energy:

potential energy = -m*y-acceleration*y-coordinate

Isn't that energy = m * acceleration * height

And for my code its:

potensial energy = particle.getMass() * particle.getAccel().y * particle.getPos().y

Share this post


Link to post
Share on other sites
Does my Integration System need to go downward a little bit when i change acceleration from -10 to 10 (or something else positive)? Velocity is 2 m/s. (In demo change the acceleration from -10 to 10 and you'll see it goes downward a little bit and then upward (look at the position here))

Thanks,
Kasya

Share this post


Link to post
Share on other sites
No, i really cant see the particle moving downwards at all with acc = +10 and vel = 2. Are you sure youre using the same version of the program as i am?

Quote:

And for my code its:

potensial energy = particle.getMass() * particle.getAccel().y * particle.getPos().y


No, you need the minus sign. Think about it. You want the potential energy to increase as you move to higher positions, right. So when y-acc is -10 and y-coord goes from 0 to 2 you will need the minus sign to get a higher energy for y=2.

Your probably thinking about the textbook example potential energy Ep = mgh but in this case they have introduced the y-direction as pointing downwards, while your y axis is pointing upwards, hence the minus sign.

//Emil

Share this post


Link to post
Share on other sites
That was my Mistake. When you put 10 it goes downward as needed. But when you put less, It doesn't go upward straight. It goes downward (Look at the position and you'll see it goes downward a little bit) then upward

That's what i saw:
I put acceleration +4.0 During Simulation (in any position). What i saw is something i don't understand it starts decreasing. Then starts increasing as needed. Its like when it starts it goes upward then goes downward. Maybe it needs to be like that.

Just please try that once and see what happens with +2, +3 and +4. And tell me do you see results like mine.

Thanks,
Kasya

Note: The Demos are same

Share this post


Link to post
Share on other sites
No, just no. I cant reproduce your problem. Maybe ive missunderstood something. I tried putting the acceleration to +4, and the particle moved upwards as expected with what looked like 2 m/s initial velocity and some acceleration. The position never dropped below the initial value.

Share this post


Link to post
Share on other sites

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

If you intended to correct an error in the post then please contact us.

Create an account or sign in to comment

You need to be a member in order to leave a comment

Create an account

Sign up for a new account in our community. It's easy!

Register a new account

Sign in

Already have an account? Sign in here.

Sign In Now

Sign in to follow this