The order of calculating velocity and position alters the solution?

[TachyonLord]

So I tried solving the differential equation for a spring - mass system using Euler's Algorithm in Python. The equation being
d²x/dt²= -4π²x

(The equation was obtained by Dimensional Analysis)
here x and t are both dimensionless equivalents of position and time.

Formula used for Energy Equivalent : 0.5*v*v + 0.5*4.*np.pi*np.pi*x*x (by dimensional analysis)

So in the loop, where we successively keep adding to the position and velocity I had 2 cases :

1. x = x + v*h;
   v = v + a*h

   here, h is the small time interval (h≈dt)
   When this is executed, I get this :
   
   
   
   Which is kinda fine, considering the fact that my end time is 10.0 and I have kept n = 10,000 (h = end_time/n)
   
2. But when the order is changed, i.e.
   

   v = v + a*h
   x = x + v*h

   here, h is the small time interval (h≈dt)
   When this is executed, I get this :
   
   
   
   As we clearly see, the energy is clearly Oscillating here whereas in the first case, it was increasing (within the limits of Euler's Algorithm, Euler's is a poor algorithm but its fine for small stuff like this)

I can't really seem to wrap my head around why the nature of the way energy is varying differs, in fact the second case has less of a margin(as in like it varies rapidly but stays within a small range)

As far as I understand, the order is just an indication of what changes first : velocity/displacement/acceleration. I just don't get it. Help ?

 

[DrClaude]

You didn't give how you calculate the acceleration, so it is hard to say exactly. Since v is used to calculate x, updating v first means that you are not calculating x the same way (or at the same time), which may be giving a higher-order term for the leading error.

 

[TachyonLord]

You didn't give how you calculate the acceleration, so it is hard to say exactly. .

acceleration = d²x/dt²= -4π²x

Since v is used to calculate x, updating v first means that you are not calculating x the same way (or at the same time), which may be giving a higher-order term for the leading error.

But then, I have initialized the varibales to their respective initial values, so the values at (t=0) are given. The successive values are calculated using Euler's Algorithm.
At the time instant t = h
we should be able to simultaneously calculate x and v. A s far as v is concerned, even v is explicitly a function of position since acceleration is a function of position.
So i still don't really get it. If you want, I can post the code, but its like kinda long.

 

[DrClaude]

When you discretize time, the order in which you do things is important. The actual equations you get from the Euler method are
$$
x(t_{i+1}) = x(t_i) + v(t_i) h \\
v(t_{i+1}) = v(t_i) + a(t_i) h \\
a(t_{i+1}) = -4 \pi x(t_{i+1})
$$

In the second case, you have instead ##x(t_{i+1}) = x(t_i) + v(t_{I+1}) h##, which definitely not the same. Look up for stance the leap-frog algorithm.

 

[TachyonLord]

When you discretize time, the order in which you do things is important. The actual equations you get from the Euler method are
$$
x(t_{i+1}) = x(t_i) + v(t_i) h \\
v(t_{i+1}) = v(t_i) + a(t_i) h \\
a(t_{i+1}) = -4 \pi x(t_{i+1})
$$

In the second case, you have instead ##x(t_{i+1}) = x(t_i) + v(t_{I+1}) h##, which definitely not the same. Look up for stance the leap-frog algorithm.

Thank you so much ! :) Although can you tell me as to why the energy is kinda oscillating, also that the error is less the wrong case and tbh seems more accurate for the same values.

As you can clearly see, the wrong method seems to be more accurate, I mean the energy for case 1 rises by 8 units whereas the enrgy in case 2 (wrong algorithm) just alternates between a range of 0.12 units (more close to actual energy conservation). Also, in the second case, the amplitude seems to remain constant(as it should) but the other case has a slope ( obviously because the energy is increasing). I honestly don't get why this is happening.(Also let's just ignore the fact that I used MS-Paint for editing and drawing those lines)

 

[TachyonLord]

I honestly need some clarification to this, can someone please respond ?

 

[FactChecker]

Since your calculations are done in a loop with small steps in time, h, the calculations should be identical except for the initial value. In one case, the position has already been incremented before the velocity is first incremented and in the other case, it has not. When there is a difference, it usually indicates that your time step, h, is not small enough to make the difference insignificant. (There are cases where the solutions will not settle down no matter how small h is.)

PS. That being said, such a large difference usually means that there is a mistake in the calculations. You should double-check them.

 

[TachyonLord]

Since your calculations are done in a loop with small steps in time, h, the calculations should be identical except for the initial value. In one case, the position has already been incremented before the velocity is first incremented and in the other case, it has not. When there is a difference, it usually indicates that your time step, h, is not small enough to make the difference insignificant. (There are cases where the solutions will not settle down no matter how small h is.)

PS. That being said, such a large difference usually means that there is a mistake in the calculations. You should double-check them.

Thank you so much ! If I decrease the time step, it obviously gets better, but what I don't seem to understand is why the energy oscillates in one case whereas it doesn't in the other. I can send the code if you want to.

 

[DrClaude]

Make a phase-space plot of the trajectory (plot p vs x for all time steps).

 

[TachyonLord]

Make a phase-space plot of the trajectory (plot p vs x for all time steps).

Could you provide me with some resources or anything because I honestly don't know what that is. Thanks a lot though.

 

[DrClaude]

Could you provide me with some resources or anything because I honestly don't know what that is. Thanks a lot though.

I said what a phase-space plot was in the parenthesis. Each time step correspond to one point on a graph where the x-axis is x and the y-axis is p (or v). For harmonic motion, this has a very particular shape. If you make such plots for your results and post them, I will help you understand what is going on with the numerical integration.

 

[FactChecker]

In your post #5 are you using the leap-frog method that @DrClaude recommended?

 

[TachyonLord]

I said what a phase-space plot was in the parenthesis. Each time step correspond to one point on a graph where the x-axis is x and the y-axis is p (or v). For harmonic motion, this has a very particular shape. If you make such plots for your results and post them, I will help you understand what is going on with the numerical integration.

Its coming out to an ellipse of sort, Although in the proper Euler's algorithm, it's coming out to be really smudgy(cos the energy isn't really conserved at the bigger interval) . That's all I understand from this.

In your post #5 are you using the leap-frog method that @DrClaude recommended?

Yes. With bit of an alteration in Case 2.

 

[Chestermiller]

This is probably a stupid question, but shouldn't the energy be constant? And, have you compared your results with the analytic solution?

 

[TachyonLord]

This is probably a stupid question, but shouldn't the energy be constant? And, have you compared your results with the analytic solution?

The energy isn't exactly constant because the Euler's algorithm is a poor algorithm and so due to the error in rounding off Taylor's expansion, there is obviously going to be an error of the order of "h²". I haven't yet compared it though. I will and I'll inform here.

 

[A.T.]

As far as I understand, the order is just an indication of what changes first : velocity/displacement. I just don't get it. Help ?

What matters is the order of all 3 updates: velocity/displacement/acceleration not just velocity/displacement. Just go though your steps, one by one, to understand why exactly you get either behavior.

 

[TachyonLord]

What matters is the order of all 3 updates: velocity/displacement/acceleration not just velocity/displacement.

I meant that. I am sorry, I'll edit my post. So ehm what changes first ? Acceleration (-4π²x) is explicitly a function of position.

 

[A.T.]

So ehm what changes first ?

In reality they all change simultaneously. In the Euler algorithm you get different inconsistencies for different update orders. For some update orders the errors are direction dependent, so for a whole oscillation period they cancel out. For some update orders the errors are direction independent, so they accumulate.

 

[TachyonLord]

In reality they all change simultaneously. In the Euler algorithm you get different inconsistencies for different update orders. For some update orders the errors are direction dependent, so for a whole oscillation period they cancel out. For some update orders the errors are direction independent, so they accumulate.

But I still don't get it as to why one caused the energy to oscillate, while the other just caused it to increase, also the fact that the supposedly "wrong" algorithm was more accurate as compared to the proper algorithm. I also tried plotting the global error(|Analytical Value - Computed Value|) vs time for both displacement and velocity for both the cases. And well, the results are weird. Here, h(time step) = 0.0001
The error is clearly very small but I can clearly say that the second algorithm seems to have negligbly low error, but then they vary werdly.

 

[Chestermiller]

I would never have solved this using forward Euler. Here is a scheme that guarantees conservation of energy and is second order in time:
$$v(t+\Delta t)-v(t)=-2\pi^2(x(t+\Delta t)+x(t))\Delta t$$
$$x(t+\Delta t)-x(t)=\frac{(v(t+\Delta t)+v(t))}{2}\Delta t$$Just solve these two simultaneous linear equations algebraically for ##v(t+\Delta t)## and ##x(t+\Delta t)## in terms of ##v(t)## and ##x(t)##, and apply the resulting difference scheme (which will, of course, be independent of the order of application).

 

[TachyonLord]

I would never have solved this using forward Euler. Here is a scheme that guarantees conservation of energy and is second order in time:
$$v(t+\Delta t)-v(t)=-2\pi^2(x(t+\Delta t)+x(t))\Delta t$$
$$x(t+\Delta t)-x(t)=\frac{(v(t+\Delta t)+v(t))}{2}\Delta t$$Just solve these two simultaneous linear equations algebraically for ##v(t+\Delta t)## and ##x(t+\Delta t)## in terms of ##v(t)## and ##x(t)##, and apply the resulting difference scheme (which will, of course, be independent of the order of application).

Thank you so much for your suggestion. I appreciate it a lot ! :) By the way, what is the maximum error for this algorithm ?
But I still kinda want to know about the reason for energy having these little disturbances. Thanks again :)

 

[FactChecker]

The order of calculation matters because there is a difference between the prior frame values and the current frame values. It makes a difference which one you use in each calculation. For instance, the velocity is not an integral of either the prior frame acceleration or the current frame acceleration. It is the integral of the continuously changing acceleration from one frame to the next. What the "leap-frog" method does is to take a (weighted?) average of frame values to get a better estimate.

 

[TachyonLord]

The order of calculation matters because there is a difference between the prior frame values and the current frame values. It makes a difference which one you use in each calculation. .

Yes ! But why does one give "an energy" which is oscillating while the other just gives a normal increasing function ?

 

[FactChecker]

The study of stability, convergence, divergence, and oscillatory behavior is, in general, rather deep. There can be a very fine line between stable behavior and exponential instability. In this example, I will have to defer to others who have more expertise.

 

[TachyonLord]

The study of stability, convergence, divergence, and oscillatory behavior is, in general, rather deep. There can be a very fine line between stable behavior and exponential instability. In this example, I will have to defer to others who have more expertise.

Thank you so so much for your help !

 

[DrClaude]

View attachment 235975Its coming out to an ellipse of sort, Although in the proper Euler's algorithm, it's coming out to be really smudgy(cos the energy isn't really conserved at the bigger interval) . That's all I understand from this.

These are known as trajectories in phase-space. The energy is related to the position in phase space. The exact solution gives an ellipse.

If you start by the result on the lower left, you will notice that the trajectory is not closed. Instead of coming back to the initial point, after one orbit the position obtained is slightly greater than it was before. The trajectory is spiralling out, and consequently the energy keeps increasing.

The result for case 2 is better, and at the scale of the plot, the trajectory appears to be closed (this might not be exactly true). Nevertheless, the trajectory is not perfect, which leads to an oscillation of the energy. You will notice that the energy oscillates twice during one orbit, corresponding to x = 1 → x = 0 → x = -1 (first oscillation), then x = -1 → x = 0 → x = 1 (second oscillation). By symmetry, it doesn't matter if you are on the bottom part of the orbit or the top part.

An algorithm that conserves energy will preserve the area of the trajectory in phase space, but the shape of the trajectory in phase space can still be deformed (compared to the exact solution).

 

[A.T.]

But I still don't get it as to why one caused the energy to oscillate, while the other just caused it to increase

I just told you:
direction dependent error -> cancels over a period -> energy oscillates
direction independent error -> accumulates over a period -> energy accumulates

To understand why the errors are direction dependent/independent in your specific cases, compare how your values change during steps for different directions. Note that "direction dependent" could also be a combination, like acceleration and velocity are opposite vs. in the same direction.

 

[TachyonLord]

An algorithm that conserves energy will preserve the area of the trajectory in phase space, but the shape of the trajectory in phase space can still be deformed (compared to the exact solution).

I just told you:
direction dependent error -> cancels over a period -> energy oscillates
direction independent error -> accumulates over a period -> energy accumulates

Thank you so so much to the both of you. I think I'm finally getting it(Kind of) Thank you for your time

 

[FactChecker]

To study this, you might want to plot everything: position, velocity, acceleration, and energy of the two cases side-by side. There should be some situations where the energy of one case is decreasing and the other is increasing. Then you can look for the cause of the difference.

 

[TachyonLord]

To study this, you might want to plot everything: position, velocity, acceleration, and energy of the two cases side-by side. There should be some situations where the energy of one case is decreasing and the other is increasing. Then you can look for the cause of the difference.

Thank you so much ! :) I'll definitely try that.

 

[Chestermiller]

Thank you so much for your suggestion. I appreciate it a lot ! :) By the way, what is the maximum error for this algorithm ?
But I still kinda want to know about the reason for energy having these little disturbances. Thanks again :)

For this scheme, the error is going to be on the order of ##(\Delta t)^2##, compared to ##\Delta t## for the forward Euler scheme.

The explicit difference equations are going to be $$v(t+\Delta t)=\frac{1-(\pi \Delta t)^2}{1+(\pi \Delta t)^2}v(t)-\frac{4\pi^2\Delta t}{1+(\pi \Delta t)^2}x(t)$$
$$x(t+\Delta t)=\frac{\Delta t }{1+(\pi \Delta t)^2}v(t)+\frac{1-(\pi \Delta t)^2}{1+(\pi \Delta t)^2}x(t)$$

Try it. You'll like it.

 

[Chestermiller]

Have you tried this difference scheme yet? You will find that it exactly conserves energy with no discretization error (and only tiny roundoff error). You will also find that this scheme is much more accurate than forward Euler, for equal values of the time step.