Estimate velocity on a spring by iterative approach

Question

The problem:

Consider a system with a mass and a spring as shown in the picture below. The stiffness of the spring and the mass of the object are known. Therefore, if the spring is stretched the force the spring exerts can be calculated from Hooke`s law and the instantaneous acceleration can be estimated from Newton´s laws of motion. Integrating the acceleration twice yields the distance the spring would move and subtracting that from the initial length results in a new position to calculate the acceleration and start the loop again. Therefore as the acceleration decreases linearly the speed levels off at a certain value (top right) Everything after that point, spring compressing & decelerating is neglected for this case.

My question is how would to go about coding that up in python. So far I have written some pseudocode.

instantaneous_acceleration = lambda x: 5*x/10    # a = kx/m
delta_time = 0.01     #10 milliseconds
a[0] = instantaneous_acceleration(12)     #initial acceleration when stretched to 12 m
v[0] = 0        #initial velocity 0 m/s
s[0] = 12       #initial length 12 m
i = 1
while a[i] > 12:
    v[i] = a[i-1]*delta_time + v[i-1]      #calculate the next velocity
    s[i] = v[i]*delta_time + s[i-1]        #calculate the next position
    a[i] = instantaneous_acceleration (s[i])          #use the position to derive the new accleration
    i = i + 1

Any help or tips are greatly appreciated.

Answer 1

You have a sign error in the force, for a spring or any other oscillation it should always be opposite to the excitation direction. Correcting this gives instantly an oscillation. However, your loop condition will now never be satisfied, so you have to also adapt that.

You can immediately increase the order of your method by elevating it from the current symplectic Euler method to Leapfrog-Verlet. You only have to change the interpretation of v[i] to be the velocity at t[i]-dt/2. Then the first update uses the acceleration in the middle at t[i-1] to compute the velocity at t[i-1]+dt/2=t[i]-dt/2 from the velocity at t[i-1]-dt/2 using a midpoint formula. Then in the next line the position update is a similar midpoint formula using the velocity at the middle time between the position times. All you have to change in the code to get this advantage is to set the initial velocity to the one at time t[0]-dt/2 using the Taylor expansion at t[0].

instantaneous_acceleration = lambda x: -5*x/10    # a = kx/m
delta_time = 0.01     #10 milliseconds
s0, v0 = 12, 0 #initial length 12 m, initial velocity 0 m/s
N=1000
s = np.zeros(N+1); v = s.copy(); a = s.copy()
a[0] = instantaneous_acceleration(s0)     #initial acceleration when stretched to 12 m
v[0] = v0-a[0]*delta_time/2        
s[0] = s0       
for i in range(N):
    v[i+1] = a[i]*delta_time + v[i]               #calculate the next velocity
    s[i+1] = v[i+1]*delta_time + s[i]             #calculate the next position
    a[i+1] = instantaneous_acceleration (s[i+1])  #use the position to derive the new acceleration
#produce plots of all these functions
t=np.arange(0,N+1)*delta_time;
fig, ax = plt.subplots(3,1,figsize=(5,3*1.5))
for g, y in zip(ax,(s,v,a)):
    g.plot(t,y); g.grid();
plt.tight_layout(); plt.show();

This is obviously and correctly an oscillation. The exact solution is 12*cos(sqrt(0.5)*t), using it and its derivatives to compute the errors in the numerical solution (remember the leap-frogging of the velocities) gives via

w=0.5**0.5; dt=delta_time;
fig, ax = plt.subplots(3,1,figsize=(5,3*1.5))
for g, y in zip(ax,(s-12*np.cos(w*t),v+12*w*np.sin(w*(t-dt/2)),a+12*w**2*np.cos(w*t))): 
    g.plot(t,y); g.grid();
plt.tight_layout(); plt.show();

the plot below, showing errors in the expected size delta_time**2.