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"}]

enter image description here

  • 3
    $\begingroup$ You're just solving a differential equation, right? Why not just use NDSolve for this? It also has special methods for preserving energy (link) $\endgroup$ Oct 3, 2023 at 15:06
  • $\begingroup$ 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. $\endgroup$
    – JEM
    Oct 4, 2023 at 2:08
  • $\begingroup$ 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. $\endgroup$ Oct 4, 2023 at 7:42

1 Answer 1


In the discussion about the momentum principle in Matter and Interactions by Chabay and Sherwood they approximate the new velocity by dividing the new momentum by the mass and then getting the new position as rf = ri + v dt. In their discussion they mention that the force is assumed to be constant over the interval dt and by taking dt smaller the resulting approximation would be better. Indeed this is the case. Here is a graph with step size decreased by 2 orders of magnitude and the number of iterations increased by 2 OOM.

enter image description here


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.