stupid me or stupid Euler method?

[HomesickAlien]

Hi all,
One more question about my Euler method problem. The usual Euler method is not energy conserving (it diverges with time), so it has to be modified. If you solve the simple harmonic oscillator with spring constant k and mass m, position and velocity should have this form:

x(t + Δt) = x(t) + Δt * v(t + Δt) + O(Δt)²
v(t + Δt) = v(t) - k * Δt * x(t) / m + O(Δt)²

So I used the following code fragment to calculate my x and v:

iX[0] = aSpring.iPosition;
iV[0] = aSpring.iVelocity;
	for (int i = 1; i <= aSteps; i++)
	{
 iX[i] = aSpring.iPosition + aSpring.iVelocity * aDt;
 iV[i]  = aSpring.iVelocity - aSpring.iK / aSpring.iM * aSpring.iPosition *aDt;
 aSpring.iPosition = iX[i];
 aSpring.iVelocity = iV[i];
	}

The problem is that I don't get a circle in the (x,v) diagram, I get a spiral which will have increased amplitude and velocity with increasing time. I cannot find any major errors in the code. Somebody has an idea?

[ne555]

When talking about Euler you need to be more precise.
According to your formula, you need the future velocity to calculate the position. But in the code you're using the current one.
Also, integer division returns an integer

[HomesickAlien]

I don't think I have to write the derivation of the formulae here. It is sufficient to know that in order to solve the divergence problem one replaces v(t) by v(t+Δt) in the first formula. I think I have covered that by setting aSpring.iVelocity = iV[i] in the code. If I put it like this it looks exactly like in the formula, i.e. the current velocity is replaced by the velocity in the future.

[ne555]

The problem is that I don't get a circle in the (x,v) diagram, I get a spiral which will have increased amplitude and velocity with increasing time

x(t + Δt) = x(t) + Δt * v(t) + O(Δt^2)
v(t + Δt) = v(t) - k * Δt * x(t) / m + O(Δt^2)
You are discarding a term.

[HomesickAlien]

The error terms O(Δt)² are not included in the implementation.

[ne555]

Yep. And those errors are making diverge the solution. (it is not stable)

[HomesickAlien]

The error for constant Δt is also constant. The divergence comes from the replacement of the the first order derivates by their forward finite differences. All I can say that "everywhere" it is stated that if you replace the t by t+Δt in the first equation you should get a stable solution which obeys energy conservation (of course with the mentioned error squared).

[ne555]

Sure, use backward Euler instead of forward Euler.
But again, your code is forward Euler. When calculating X[K] you are using V[K-1], it should be V[K]

The error for constant Δt is also constant

You can't say that the error is constant. It depends on the step, but also on a derivative of the function evaluated in some point of the interval.

[HomesickAlien]

The problem is how to use v[i] for the first equation while using v[i-1] for the second.

[ne555]

Like this

for (int K = 1; K <= aSteps; ++K)
{
  V[K]  = V[K-1] - X[K-1]*C*Dt;
  X[K] = X[K-1] + V[K]*Dt; //as we need the calculated value for V, this need to be after.
}

[HomesickAlien]

Oh, yes that's cool, exchanging the lines works and I am getting a perfect circle, thanks!!!