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

k_ = k + 3*damp/h + 6*mass/h/h

cv = 3*damp + 6*mass/h
ca = damp*h/2 + 3*mass

x_ = [] ; v_ = [] ; t_ = []

t = 0.00 ; X = 0.00 ; V = 0.00 ; P = p(t)
A = (P - V*damp - X*k)/mass

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)
    dp_ = (Ph-P) + cv*V + ca*A
    dx  = dp_/k_
    dv  = 3*dx/h - 3*V - A*h/2
    X = X+dx ; V = V+dv
    P = Ph ; A = (P - damp*V - k*X)/mass
plot(t_, x_)


    

    




Populating the interactive namespace from numpy and matplotlib






    
Out[1]:





[<matplotlib.lines.Line2D at 0x7f7d57ba2750>]






    

In [1]:

    




    

In [ ]: