Mass_and_spring_system.bas


'Two masses suspended by two springs in tandem are subject to a driving motion.
'The driving motion is sinusoidal, with an increasing frequency characterized by an increasing sweep rate.
'The system is described by mass and by spring length and stiffness.
'The output graph shows the driver motion on top, and the response of the two masses below.
'Notice the two resonance intervals where mass displacement is greatest.
'From the left, the masses are initially in phase, and later resonate out of phase.

WINDOW 0, -1200,1000, 100
'write function for the differential equations of motion.
func osct (mt, kt, kb, tn, bn, xt, xb, fr, zt) = -kt * (xt - tn + fr) / mt + kb * (xb - xt - bn) / mt - .1 * zt
func
oscb (mb, kb, nb, xt, xb, zb) = -kb * (xb - xt - bn) / mb - .1 * zb

'input initial conditions
PRINT "driven coupled masses with dampening program "
PRINT "vertical range is (top) 0 to (bottom) 120"
INPUT "top mass initial speed, m/s (0)"; zt
INPUT "top mass initial position, m (40)"; xt
INPUT "top spring natural length, m (40)"; tn
INPUT "top mass, kg (40)"; mt
INPUT "top spring constant, N/m (40)"; kt
INPUT "bottom mass initial speed, m/s (0)"; zb
INPUT "bottom mass initial position, m (80)"; xb
INPUT "bottom spring natural length, m (40)"; bn
INPUT "bottom mass, kg (30)"; mb
INPUT "bottom spring constant, N/m (30)"; kb
INPUT "final time (300)"; tf: n = 60 * tf: h = tf / n
INPUT "frequency sweeprate, Hz/sec (.004)"; hz

'time step loops to find the system response
t = 0: fact = hz * tf / n: CLS
FOR
q = 1 TO n
fr = 5 * SIN(fact * q * t)
kt1 = h * zt: lt1 = h * osct(mt, kt, kb, tn, bn, xt, xb, fr, zt)
kb1 = h * zb: lb1 = h * oscb(mb, kb, bn, xt, xb, zb)
kt2 = h * (zt + .5 * lt1)
kb2 = h * (zb + .5 * lb1)
lt2 = h * osct(mt, kt, kb, tn, bn, xt + .5 * kt1, xb + .5 * kb1, fr, zt)
lb2 = h * oscb(mb, kb, bn, xt + .5 * kt1, xb + .5 * kb1, zb)
kt3 = h * (zt + .5 * lt2)
kb3 = h * (zb + .5 * lb2)
lt3 = h * osct(mt, kt, kb, tn, bn, xt + .5 * kt2, xb + .5 * kb2, fr, zt)
lb3 = h * oscb(mb, kb, bn, xt + .5 * kt2, xb + .5 * kb2, zb)
kt4 = h * (zt + lt3)
kb4 = h * (zb + lb3)
xt = xt + (kt1 + 2 * kt2 + 2 * kt3 + kt4) / 6
xb = xb + (kb1 + 2 * kb2 + 2 * kb3 + kb4) / 6
zt = kt4 / h: zb = kb4 / h: t = t + h
xtp = -xt: xbp = -xb: tme = 100 * q / n

'plot the graph
PSET 10*tme, 10*fr
PSET 10*tme, 10*xtp, 9
PSET 10*tme, 10*xbp, 3
NEXT
q
end

Comments

Hi Dave,

Thanks for offering these programs - very interesting. Hope I pasted the correct version, there were a few duplicate submissions.

I've got a suggestion for the handing default input. This avoids needing to enter the default values:


func get_input(prompt, default)
input prompt + " (" + default + ")", result
if (result == 0) then
result = default
fi
get_input = result
end

kt = get_input("top spring constant, N/m", 40)
? kt