🔒 Locked

integration

Started by juxie Feb 23, 2008 at 3:32 PM 6 replies 1.4k views

Original Post

Feb 23, 2008 at 3:32 PM

Hi, I am a bit confused now. Should RK4 or Euler actually calculates the new positions first and then the new velocity or the other way around? Thanks.

chronos

Feb 23, 2008 at 8:19 PM

You integrate velocity to obtain position.

taby

Feb 23, 2008 at 9:17 PM

Good question. Always add acceleration to obtain the new velocity first, then add the new velocity to obtain the new position. If you do it the other way around, things will get wonky.

Here are source snippets of both RK4 and Euler doing the same thing:

// Runge-Kutta 4 integration, GRinline void proceed_RK4(const double &dt){	// Performance is increased on some architectures when division is avoided	static const double one_sixth = 1.0/6.0;	vector_3 k1_velocity = mercury.velocity;	vector_3 k1_acceleration = acceleration(mercury.position, k1_velocity);	vector_3 k2_velocity = mercury.velocity + k1_acceleration*dt*0.5;	vector_3 k2_acceleration = acceleration(mercury.position + k1_velocity*dt*0.5, k2_velocity);	vector_3 k3_velocity = mercury.velocity + k2_acceleration*dt*0.5;	vector_3 k3_acceleration = acceleration(mercury.position + k2_velocity*dt*0.5, k3_velocity);	vector_3 k4_velocity = mercury.velocity + k3_acceleration*dt;	vector_3 k4_acceleration = acceleration(mercury.position + k3_velocity*dt, k4_velocity);	mercury.velocity += (k1_acceleration + k2_acceleration*2.0 + k3_acceleration*2.0 + k4_acceleration)*dt*one_sixth;	mercury.position += (k1_velocity + k2_velocity*2.0 + k3_velocity*2.0 + k4_velocity)*dt*one_sixth;}// Euler integration, GRinline void proceed_E(const double &dt){	mercury.velocity += acceleration(mercury.position, mercury.velocity)*dt;	mercury.position += mercury.velocity*dt;}

chronos

Feb 24, 2008 at 1:46 AM

Ah... I misunderstood the question.

Tabby,

What is your source for the claim that one should obtain the new velocity before the new position? According to the book "Numerical Methods" by Fausett, the way to solve higher order ODE's is as follows:

If y'''' = g(x,y,y') is your second-order ODE, then

v = y'

u = y

v' = g(x,u,v)

u' = v = f(x,u,v).

By Euler's method:

a = initial value of x;b = final value of x;n = number of steps;h = (b - a) / n; // step sizeu[0] = u0;       // initial condition for u at av[0] = v0;       // initial condition for v at afor i = 0 to n-1    x = a + h * i;    u[i+1] = u + h * f(x, u, v);    v[i+1] = v + h * g(x, u, v);next ireturn u, v; // arrays of approximate values

So, both the new position and the new velocity are obtained at the same time.

[Edited by - chronos on February 24, 2008 2:46:46 AM]

juxie

Feb 24, 2008 at 2:32 AM

Quote:
Original post by chronos
Ah... I misunderstood the question.

Tabby,

What is your source for the claim that one should obtain the new velocity before the new position? According to the book "Numerical Methods" by Fausett, the way to solve higher order ODE's is as follows:

If y'''' = g(x,y,y') is your second-order ODE, then

v = y'

u = y

v' = g(x,u,v)

u' = v = f(x,u,v).

By Euler's method:

*** Source Snippet Removed ***

So, both the new position and the new velocity are obtained at the same time.

Should I then do Euler as:

proceedEuler(float dt){    nextPos = initPos + initVel * dt;    nextVel = initVel + acc * dt;}

MrRowl

Feb 24, 2008 at 5:55 AM

Quote:
Original post by juxie
Should I then do Euler as:
*** Source Snippet Removed ***

Yes - this is the standard Euler formulation.

There's a slight variation, called either Euler–Cromer or Symplectic Euler, where you use your predicted velocity as an estimate that is used in the position update:

nextVel = initVel + acc * dt;
nextPos = initPos + nextVel * dt;

This is generally considered to be better, and is a trivial change to make.

Higher order methods don't allow you to seperate out the steps into a simple "do velocity" and "do position" update - the two things are tied together as one operation (calculate the new state form the old state, where state means position and velocity).

Higher order integration methods differ in their details - but also in one important respect: If you have a system of N bodies, some of which are connected by forces that depend on the relative position between bodies, you can perform an Euler integration on each one independantly, using forces calculated at the start of the integration. This means it's OK to write a Body::integrate(dt) function, for example.

When implementing RK4 (for example), the system needs to be integrated as a whole, since the inter-object forces change "during" the integration. This means you can't implement RK4 for individual bodies - you need to write a solver that integrates the whole system at once. If you expect you need to use RK4 (or similar) eventually, its worth bearing this in mind so that you don't write a lot of code that assumes each body can be integrated independantly.

juxie

Feb 24, 2008 at 3:44 PM

Quote:
Original post by MrRowl
Quote:
Original post by juxie
Should I then do Euler as:
*** Source Snippet Removed ***

Yes - this is the standard Euler formulation.

There's a slight variation, called either Euler–Cromer or Symplectic Euler, where you use your predicted velocity as an estimate that is used in the position update:

nextVel = initVel + acc * dt;
nextPos = initPos + nextVel * dt;

This is generally considered to be better, and is a trivial change to make.

Higher order methods don't allow you to seperate out the steps into a simple "do velocity" and "do position" update - the two things are tied together as one operation (calculate the new state form the old state, where state means position and velocity).

Higher order integration methods differ in their details - but also in one important respect: If you have a system of N bodies, some of which are connected by forces that depend on the relative position between bodies, you can perform an Euler integration on each one independantly, using forces calculated at the start of the integration. This means it's OK to write a Body::integrate(dt) function, for example.

When implementing RK4 (for example), the system needs to be integrated as a whole, since the inter-object forces change "during" the integration. This means you can't implement RK4 for individual bodies - you need to write a solver that integrates the whole system at once. If you expect you need to use RK4 (or similar) eventually, its worth bearing this in mind so that you don't write a lot of code that assumes each body can be integrated independantly.

Thanks it's clearer now.
The reason I am asking this is also because my current Runge-Kutta method for cloth simulation has been acting weird with collision detection.
The particles will stick to object it is colliding with and don't response correctly to change in velocity due to impact.
Encountered this kind of problem before?
How can it be solved or rather how can I implement Runge-Kutta with correct collision detection?

Thank you.

juxie

Mar 2, 2008 at 8:04 AM

Quote:
Original post by taby
Good question. Always add acceleration to obtain the new velocity first, then add the new velocity to obtain the new position. If you do it the other way around, things will get wonky.

Here are source snippets of both RK4 and Euler doing the same thing:

*** Source Snippet Removed ***

Hi Taby,

Just to ask.
I was told that K1 in RK4 is supposed to be the same as Euler.
Which means if I replace k2_velocity, k3_velocity, k4_velocity with k1_velocity like below, I should get back Euler.

mercury.velocity += (k1_acceleration + k1_acceleration*2.0 + k1_acceleration*2.0 + k1_acceleration)*dt*one_sixth;	mercury.position += (k1_velocity + k1_velocity*2.0 + k1_velocity*2.0 + k1_velocity)*dt*one_sixth;

From your code, it seems that K1 and Euler are not the same, aren't they?

[Edited by - juxie on March 4, 2008 1:04:20 AM]

Topic Locked

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