NRL-Plasma-Physics-Division / block-on-spring

TurboPy app to compute the motion of a block on a spring
0 stars 5 forks source link

Error In ForwardEuler Push Method #23

Closed RushilSidhu closed 3 years ago

RushilSidhu commented 3 years ago

The FowardEuler class in the spring.py file has an error. The push function copies the current momentum, then updates the momentum to a new one, then uses the copied momentum in the calculation for the position. This seems intentional but over a long period of time causes inconsistency with the calculations. This is most easily seen by going from t=0 to t=100 with 1000 time steps and graphing position vs. momentum. This should be a circle, but in this case, it becomes a spiral, which is obviously not correct. The simplest solution is swapping the order of calculation of momentum and position in the ForwardEuler class.

arichar6 commented 3 years ago

This is not an error in the implementation in spring.py. The spiraling behavior that you describe is one of the consequences of using the Euler Method. It is known that this method has low-order numerical errors which lead to the energy not being conserved (and hence the spiral in position-momentum phase space).

The solution that you propose is a variant of the Leapfrog Method. This is a much better numerical algorithm for modeling the motion of a system such as this block on a spring, which should be conserving energy. A slightly different variant of the Leapfrog Method is implemented in spring.py in the Leapfrog class. You can select it as an alternative ComputeTool for your simulation, and show that the results are much better than those obtained with the Euler Method.

The mathematics involved in analyzing the errors in these various update algorithms is actually really fascinating. It turns out that the Leapfrog Method has nice properties because it is a "symplectic" integrator. That means that it explicitly conserves areas in phase space. More information is in the wiki article on symplectic integration, and I'd be happy to chat with you about it sometime.