# Mass connected to spring on a frictionless surface vibrating out of control

I am not good at Physics so please bear with me. I am trying to understand harmonics and in order to do that I thought its best to ensure my understanding of the equations of motion is correct. To do this I tried to model a mass vibrating on a frictionless surface with restoring force from a spring when acted on by a force F. The parameters are as follows:

M = 5 kg (Mass)
F = 10 kN (Initial force)
k = -1kN/m (Spring constant)


I expected the plot of displacement against time to be uniform i.e without dissipation (or no damping). However, the result shows that the displacement are increasing. Have I made a mistake, if not, could someone explain the phenomena.

The plot is shown below Algorithm is as follows:
F = 10kN
m = 5kg
k = -1kN/m
initial acceleration, a0 = 10/5 = 2 m/s

Assuming displacement is positive to the right
At t = 0,
u = 0
v = 0
s = 0
dt = 1
a0 = 2
ar = 0
anet= a0-ar
While (t<100),
ds = u*dt + 0.5anet*dt^2;
s = s + ds;
v = u + anet*dt;
t = t + dt
acceleration by restoring force, ar= (k*s/m);
updating net acceleration, anet = anet + ar; //K is already negative
Next;


Update: I plotted the same graph at different time steps.

• What formula are you using? – Kyle Kanos Aug 8 '17 at 2:40
• I am using the equations of motion, to summarize I use, s = ut + 0.5at^2, followed by v = u+at, and update of restoring force (k*s)/m, where k and m are spring constant and mass, respectively, for different time steps. I will add the precise notations and sequence in the original question after some time. – Stephen Jacob Aug 8 '17 at 3:50
• The equations you are using for s and v are for constant acceleration, which is not the case here. Instead, you could alculate s(t) =s(t-1)+delta_s – Wolphram jonny Aug 8 '17 at 5:37
• compsci.stackexchange.com – peterh - Reinstate Monica Aug 8 '17 at 6:08
• I'm voting to close this question as off-topic because this is debugging code. – Jon Custer Aug 8 '17 at 12:58

The problem is that you are calculating displacement on the assumption of a constant acceleration in this step:

From t=1 to t=n,
ds = ut + 0.5anet*t^2;


Because you are letting t vary from 1 to n. If you instead used a time step dt (much less than the period of the motion) you could then compute the updated acceleration, velocity and position. Example Python code:

# integrating equation of motion - SHO
import numpy as np
import matplotlib.pyplot as plt

m = 1.
k = 1.
N = 601  # steps
X0 = 1.  # initial position

t = np.linspace(0,50,N)
dt = t-t

vx = 0.  # initial velocity
x  = X0  # initial position

# simplistic integration: calculate a, v, x for each step
X1 = np.zeros(N)
for ii in range(N):
a = -k*x
vx = vx + a * dt
x += vx * dt
X1[ii] = x

plt.figure()
plt.plot(t, X1)
plt.title('simple harmonic motion')
plt.xlabel('time')
plt.ylabel('displacement')
plt.show()


Plotting the result: There are fancier algorithms to do this properly, and I strongly recommend you learn about them (from Kyle's answer as a good start) - but your problem is literally that instead of using a time step, you are using the total time since the start of the simulation.

## Update

You might find it helpful to read this question with associated answers and links for a better understanding of integration algorithms and how they affect stability and appropriate choice of time steps.

• I used the same algorithm as I used earlier and plotted it for shorter timesteps from 1 to 0.005 (from @KyleKano's answer) and it appears to behave as expected. So may I conclude that it depends on the selection of the timestep? (Please see updated question) – Stephen Jacob Aug 9 '17 at 1:39
• Note that the original code had t, not dt. That changes everything. – Floris Aug 9 '17 at 2:54
• You are right @Floris, I should have paid more attention to the clarity of the algorithm. I apologize for that. I thought it had more to do with my understanding of Physics, than the value of time-step. I would still like to pose the same question that I posed to Kyle, which is if their is a way to determine an appropriate time-step? And additionally, you mentioned dt being much smaller than the period, how would i go about determining that? – Stephen Jacob Aug 9 '17 at 4:35
• The correct size of the time step depends on the algorithm. With a simple algorithm like yours, your step needs to be small enough that curvature in the function is small over the step size. When you use "higher order" algorithms (for example 4th order Runge-Kutta, look it up) then you can take larger steps because you worry about the size step where the fourth derivative changes significantly. Incidentally the period is given by $2\pi\sqrt{\frac{m}{k}}$ - when m=k=1 (as in my example) period is about 6.2 seconds and 100 steps per period is plenty small. – Floris Aug 9 '17 at 11:12

When working with a situation like this (i.e. time-varying force), it's best to start with the differential equations \begin{align} \frac{{\rm d}x}{{\rm d}t}&=v \\ \frac{{\rm d}v}{{\rm d}t}&=F/m \end{align} rather than trying to use positions and velocities from kinematics.

The above can be re-witten using a time-stepping technique called leapfrog integration, in which you update positions and velocities at offset intervals: \begin{align} x^{n+1} &= x^n + v^{n+\frac12}\Delta t \\ a^{n+1} &= F(x^{n+1}) \tag{1}\\ v^{n+\frac32}&=v^{n+\frac12} + a^{n+1}\Delta t \end{align} Here, the fractional $n$ can be thought of the "cell wall" value (e.g., $x_{i+\frac12}=\frac12(x_i+x_{i+1})$) but in time instead of space.

Or you can use a modification called verlocity verlet. This turns (1) into a multi-step process: \begin{align} a_1 &= F\left(x^n_i\right)/m \\ x^{n+1} &= x_i^n + \left(v_i^n + \frac{1}{2}\cdot a_1\cdot\Delta t\right)\cdot\Delta t \tag{2}\\ a_2 & =F\left(x^{n+1}\right)/m \\ v^{n+1} &= v_i^n + \frac{1}{2}\left(a_1+a_2\right)\cdot\Delta t \end{align}

Both integration methods are reasonably simple to implement, requiring definitions for forces and some vectors. The latter of the two is more often the algorithm-of-choice due to it's symplectic nature.

The implementation for (2) is basically

while t < t_end
a1 = Force(x) / m
x += (v + 0.5 * a1 * dt) * dt
a2 = Force (x) / m
v += 0.5 * (a1 + a2) * dt
t += dt
output t, x, v


where dt is your time-step (ought to be small, 0.005 should suffice for your values).

• Some text overlap with two other of my answers: physics.stackexchange.com/a/240475/25301 and physics.stackexchange.com/a/116688/25301 – Kyle Kanos Aug 8 '17 at 15:05
• I used the same algorithm as I used earlier and plotted it for shorter timesteps from 1 to 0.005 and it appears to behave as expected. So may I conclude that it depends on the selection of the timestep? (Please see updated question) – Stephen Jacob Aug 9 '17 at 1:38
• @StephenJacob: yes, the timestep matters for some algorithms – Kyle Kanos Aug 9 '17 at 1:39
• So is there a definite way to determine in an appropriate selection of time-step? – Stephen Jacob Aug 9 '17 at 1:41
• @StephenJacob: For hydrodynamics there is, but not 100% sure for this. I can't imagine it's trial-and-error, but I don't know what the actual method is. – Kyle Kanos Aug 9 '17 at 1:45