# Tracking сode performance with the conservation of energy

I am working with the Momentum Principle as defined by Chabay and Sherweood in their text Matter and Interactions (pf = pi+Fnet dt). The following code has been created to automate the process of updating the position and velocities of a mass on a spring in 3D. The mass is placed at the end of the spring and then the spring is pulled to its initial position r[0] and released. Since the forces involved are conservative, the total energy should be constant. My plot of the Energy is non-constant. I'm hoping someone would critique this program. Thanks

Clear[m, k, L0, r, v, p, dt, Fg, En]

m = .5;
k = 20.;
L0 = {0., .50, 0.};
r[0] = {.2, .7, -.1};
s[0] = Norm[r[0]] - Norm[L0];
Fg = m {0., -9.8, 0.};
p[0] = {0., 0., 0.};
v[0] = p[0]/m;
dt = .01;
t[0] = 0;
En[0] = .5 m v[0] . v[0] +
Norm[Fg] r[0][[2]] + .5 k (Norm[r[0] - Norm[L0]])^2

data = Table[
{s = Norm[r[n]] - Norm[L0],
p[n + 1] = p[n] + (Fg - (k s r[n])/Norm[r[n]]) dt,
v[n + 1] = p[n + 1]/m,
r[n + 1] = r[n] + v[n + 1] dt,
t[n + 1] = t[n] + dt,
En[n + 1]  = .5 m v[n + 1] . v[n + 1] +
Norm[Fg] r[n + 1][[2]] + .5 k s^2},
{n, 0, 2000}];

ListPlot[data[[All, 6]], AxesLabel -> {"Time","Energy"}]


• You're just solving a differential equation, right? Why not just use NDSolve for this? It also has special methods for preserving energy (link) Oct 3, 2023 at 15:06
• My students are learning about the momentum principle and how to iteratively update the position of a system. NDSolve will come eventually. This is a check on my own understanding as we haven't talked about energy yet.
– JEM
Oct 4, 2023 at 2:08
• In that case: what you're seeing here is that the energy isn't constant because of the integration approximation you choose for this problem. If you're going to invent your own NDSolve, you'll run into all the problems that the field of numeric DE solving had to overcome in the past. Which is why, generally, it's a good idea to just use the right tool for the job. Oct 4, 2023 at 7:42