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

Question

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.

TimeStep=0.1s

TimeStep=0.005s

Answer 1

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).

Answer 2

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[1]-t[0]

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.