In [1]:

    

%pylab inline

def p(t):
    if t < 1.00 : return 4E5 * t
    if t < 3.00 : return 2E5 * (3-t)
    return 0.00

mass = 6E05
T_n  = 0.60
wn   = 2*pi/T_n
k    = mass*wn**2
zeta = 0.02
wd   = wn * sqrt(1.00-zeta**2)
damp = 2*zeta*mass*wn

h = 0.025
cz = cos(wd*h)*exp(-zeta*wn*h)
sz = sin(wd*h)*exp(-zeta*wn*h)
              
x_ = [] ; v_ = [] ; t_ = []

t = 0.00 ; X = 0.00 ; V = 0.00 ; P = p(t)

while t < 6.00:
    x_.append(X) ; v_.append(V) ; t_.append(t)
    # print "%6.3f   %+12.10f %+12.10f" % (t, X, V)
    t = t+h
    Ph = p(t)
    dx = P/k
    ddx = (Ph-P)/k
    B = X + 2*zeta*ddx/wn/h - dx
    A = (V + zeta*wn*B - ddx/h)/wd
    X = A*sz + B*cz + dx + ddx*(1-2*zeta/wn/h)
    V = (A*(wd*cz-zeta*wn*sz) -
         B*(wd*sz+zeta*wn*cz) + ddx/h)
    P = Ph
plot(t_, x_); grid()


    

    




Populating the interactive namespace from numpy and matplotlib






    

In [2]:

    

plot(t_, v_);grid()


    

In [ ]: