Zane

05-29-2009, 09:17 PM

I am programming in python and know some C++ as well so I can interpret something similar to those two. I have written a programm that models a mass on a spring in harmonic motion with a damping force. It displays a graph in phase space and on a velocity and position vs. time graph. I would like to add in an extra force that will only activate if the amplitude of the velocity and position curve drops below some minimum. Essentially like giving the spring a little push when it's velocity drops too low. Any help on how to do so would be greatly appreciated. The program is below in python programming language which is similar to C++, comments follow #

from scipy import *

from pylab import *

k=0.05 #Define constants (k=spring constant,b=damping force,m=mass on spring)

b=0.1

m=5

npts= 100000

v=zeros(npts) #Define arrays for velocity, position, and time

x=zeros(npts)

t=zeros(npts)

v[0]=30.00 #Initialize arrays

x[0]=5.00

t[0]=0.00

dt=0.01

for i in arange(npts-1):

v[i+1]=v[i]+dt*(-k*x[i]-b*v[i]+sin(0.01*t[i]))/m #Difference equation for velocity

x[i+1]=x[i]+ v[i]*dt #Difference equation for position

t[i+1]=t[i]+dt #time euation

figure(1) #the phase plot

plot(v,x)

grid(True) #shows grid on the graph

ylabel('velocity')

xlabel('position')

figure(2) #time plots of v and x

plot(t,v,'-',label='velocity')

plot(t,x,'--',label='position')

legend() #shows the labels on graph

xlabel('Time')

show()

from scipy import *

from pylab import *

k=0.05 #Define constants (k=spring constant,b=damping force,m=mass on spring)

b=0.1

m=5

npts= 100000

v=zeros(npts) #Define arrays for velocity, position, and time

x=zeros(npts)

t=zeros(npts)

v[0]=30.00 #Initialize arrays

x[0]=5.00

t[0]=0.00

dt=0.01

for i in arange(npts-1):

v[i+1]=v[i]+dt*(-k*x[i]-b*v[i]+sin(0.01*t[i]))/m #Difference equation for velocity

x[i+1]=x[i]+ v[i]*dt #Difference equation for position

t[i+1]=t[i]+dt #time euation

figure(1) #the phase plot

plot(v,x)

grid(True) #shows grid on the graph

ylabel('velocity')

xlabel('position')

figure(2) #time plots of v and x

plot(t,v,'-',label='velocity')

plot(t,x,'--',label='position')

legend() #shows the labels on graph

xlabel('Time')

show()