<SYSTEM_TASK:>
Given the following text description, write Python code to implement the functionality described below step by step
<END_TASK>
<USER_TASK:>
Description:
 Exact Integration for an EP SDOF System
We want to compute the response, using the constant
acceleration algorithm plus MNR, of an Elasto Plastic (EP)
system... but how we can confirm or reject our results?
It turns out that computing the exact response of an EP system with 
a single degree of freedom is relatively simple.
Here we discuss a program that computes the analytical solution of
our problem.
The main building blocks of the program will be two functions
that compute, for the elastic phase and for the plastic phase,
the analytical functions that give the displacement and the
velocity as functions of time.
Elastic response
We are definining a function that, for a linear dynamic system,
returns not the displacement or the velocity at a given time,
but rather a couple of functions of time that we can use afterwards
to compute displacements and velecities at any time of interest.
The response depends on the parameters of the dynamic system $m,c,k,$
on the initial conditions $x_0, v_0,$ and on the characteristics of the external load.
Here the external load is limited to a linear combination of a cosine modulated,
a sine modulated (both with the same frequency $\omega$) and a constant force,
<center>$P(t) = c_C \cos\omega t + c_S \sin\omega t + F,$</center>
but that's all that is needed for the present problem.
The particular integral being
<center>$\xi(t) = S \cos\omega t + R \sin\omega t + D,$</center>
substituting in the equation of motion and equating all the corresponding
terms gives the undetermined coefficients in $\xi(t)$, then evaluation of the general
integral and its time derivative for $t=0$ permits to find the constants in the homogeneous
part of the integral.
The final step is to define the displacement and the velocity function, according to 
the constants we have determined, and to return these two function to the caller
Step1:  Plastic response
In this case the equation of motion is
<center>$m\ddot x + c \dot x = P(t),$</center>
the homogeneous response is 
<center>$x(t)=A\exp(-\frac{c}{m}t)+B,$</center>
and the particular integral, for a load described
as in the previous case, is (slightly different...)
<center>$\xi(t) = S \cos\omega t + R \sin\omega t + Dt.$</center>
Having computed $R, S,$ and $D$ from substituting $\xi$ 
in the equation of motion, $A$ and $B$ by imposing 
the initial conditions,we can define the displacement and velocity 
functions and, finally, return these two functions to the caller.
Step2:  An utility function
We need to find when
the spring yields
the velocity is zero
to individuate the three ranges of different behaviour
elastic
plastic
elastic, with permanent deformation.
We can use the simple and robust algorithm of bisection 
to find the roots for 
<center>$x_{el}(t)=x_y \text{ and } \dot{x}_{ep}(t)=0$.</center>
Step3:  The system parameters
Step4:  Derived quantities
The damping coefficient $c$ and the first yielding displacement, $x_y$.
Step5:  Load definition
Our load is a half-sine impulse
<center>$p(t)=\begin{cases}p_0\sin(\frac{\pi t}{t_1})&0\leq t\leq t_1,\\
    0&\text{otherwise.}\end{cases}$</center>
In our exercise
Step6:  The actual computations
Elastic, initial conditions, get system functions
Step7:  Yielding time is
The time of yielding is found solving the equation $x_\text{next}(t) = x_y$
Step8:  Forced response in elastic range is
Step9:  Preparing for EP response
First, the system state at $t_y$ is the initial condition for the EP response
Step10:  now, the load must be expressed in function of a restarted time,
<center>
    $\tau=t-t_y\;\rightarrow\;t=\tau+t_y\;\rightarrow\;\sin(\omega t)=\sin(\omega\tau+\omega t_y)$
 </center> <center>
    $\rightarrow\;\sin(\omega t)=\sin(\omega\tau)\cos(\omega t_y)+\cos(\omega\tau)\sin(\omega t_y)$
</center>
Step11:  Now we generate the displacement and velocity functions for the yielded phase,
please note that the yielded spring still exerts a constant force $f_y$ on the mass,
and that this fact must be (and it is) taken into account.
Step12:  At this point I must confess that I have already peeked the numerical solution,
hence I know that the velocity at $t=t_1$ is still greater than 0 and I know that the current
solution is valid in the interval $t_y\leq t\leq t_1$.
Step13:  In the next phase,  still it is $\dot x> 0$ so that
the spring is still yielding,
but now $p(t)=0$, so we must compute two new state functions,
starting as usual from the initial conditions (note that the yielding
force is still applied)
Step14:  Elastic unloading
The only point worth commenting is the constant force
that we apply to our system.
The force-displacement relationship for an EP spring is 
<center>$f_\text{E} = k(x-x_\text{pl})= k x - k (x_\text{max}-x_y)$</center>
taking the negative, constant part of the last expression
into the right member of the equation of equilibrium we have a constant
term, as follows
Step15:  now we are ready to plot the whole response
Step17:  Numerical solution
first, we need the load function
Step18:  and also a function that, given the displacement,
the velocity and the total plastic deformation, returns
the stiffness and the new p.d.; this function is defined
in terms of the initial stiffness and the yielding load
Step19:  Problem data
Step20:  Initialize the algorithm
compute the functions that return the load and the tangent sstiffness + plastic deformation
compute the damping constant
for a given time step, compute all the relevant algorithmic constants, with $\gamma=\frac12$ and
   $\beta=\frac14$
Step21:  System state initialization
and a bit more, in species we create two empty vectors to hold the computation results
Step22:  Iteration
We iterate over time and, if there is a state change, over the single time step to equilibrate
the unbalanced loadings
Step23:  Plotting our results
we plot red crosses for the numericaly computed response and a continuous line for the results
of the analytical integration of the equation of motion.
<ASSISTANT_TASK:>
Python Code:
def resp_elas(m,c,k, cC,cS,w, F, x0,v0):
    wn2 = k/m ; wn = sqrt(wn2) ; beta = w/wn
    z = c/(2*m*wn)
    wd = wn*sqrt(1-z*z)
    # xi(t) = R sin(w t) + S cos(w t) + D
    det = (1.-beta**2)**2+(2*beta*z)**2
    R = ((1-beta**2)*cS + (2*beta*z)*cC)/det/k
    S = ((1-beta**2)*cC - (2*beta*z)*cS)/det/k
    D = F/k
    A = x0-S-D
    B = (v0+z*wn*A-w*R)/wd
    
    def x(t):
        return exp(-z*wn*t)*(A*cos(wd*t)+B*sin(wd*t))+R*sin(w*t)+S*cos(w*t)+D
    
    def v(t):
        return (-z*wn*exp(-z*wn*t)*(A*cos(wd*t)+B*sin(wd*t))
                  +wd*exp(-z*wn*t)*(B*cos(wd*t)-A*sin(wd*t))
                  +w*(R*cos(w*t)-S*sin(w*t)))
    return x,v   
Explanation: Exact Integration for an EP SDOF System
We want to compute the response, using the constant
acceleration algorithm plus MNR, of an Elasto Plastic (EP)
system... but how we can confirm or reject our results?
It turns out that computing the exact response of an EP system with 
a single degree of freedom is relatively simple.
Here we discuss a program that computes the analytical solution of
our problem.
The main building blocks of the program will be two functions
that compute, for the elastic phase and for the plastic phase,
the analytical functions that give the displacement and the
velocity as functions of time.
Elastic response
We are definining a function that, for a linear dynamic system,
returns not the displacement or the velocity at a given time,
but rather a couple of functions of time that we can use afterwards
to compute displacements and velecities at any time of interest.
The response depends on the parameters of the dynamic system $m,c,k,$
on the initial conditions $x_0, v_0,$ and on the characteristics of the external load.
Here the external load is limited to a linear combination of a cosine modulated,
a sine modulated (both with the same frequency $\omega$) and a constant force,
<center>$P(t) = c_C \cos\omega t + c_S \sin\omega t + F,$</center>
but that's all that is needed for the present problem.
The particular integral being
<center>$\xi(t) = S \cos\omega t + R \sin\omega t + D,$</center>
substituting in the equation of motion and equating all the corresponding
terms gives the undetermined coefficients in $\xi(t)$, then evaluation of the general
integral and its time derivative for $t=0$ permits to find the constants in the homogeneous
part of the integral.
The final step is to define the displacement and the velocity function, according to 
the constants we have determined, and to return these two function to the caller
End of explanation
def resp_yield(m,c, cC,cS,w, F, x0,v0):
    # csi(t) = R sin(w t) + S cos(w t) + Q t
    Q = F/c
    det = w**2*(c**2+w**2*m**2)
    R = (+w*c*cC-w*w*m*cS)/det
    S = (-w*c*cS-w*w*m*cC)/det
    # x(t) = A exp(-c t/m) + B + R sin(w t) + S cos(w t) + Q t
    # v(t) = - c A/m exp(-c t/m) + w R cos(w t) - w S sin(w t) + Q
    #
    # v(0) = -c A / m + w R + Q = v0
    A = m*(w*R + Q - v0)/c
    # x(0) = A + B + S = x0
    B = x0 - A - S
    
    def x(t):
        return A*exp(-c*t/m)+B+R*sin(w*t)+S*cos(w*t)+Q*t
    def v(t):
        return -c*A*exp(-c*t/m)/m+w*R*cos(w*t)-w*S*sin(w*t)+Q
    
    return x,v
Explanation: Plastic response
In this case the equation of motion is
<center>$m\ddot x + c \dot x = P(t),$</center>
the homogeneous response is 
<center>$x(t)=A\exp(-\frac{c}{m}t)+B,$</center>
and the particular integral, for a load described
as in the previous case, is (slightly different...)
<center>$\xi(t) = S \cos\omega t + R \sin\omega t + Dt.$</center>
Having computed $R, S,$ and $D$ from substituting $\xi$ 
in the equation of motion, $A$ and $B$ by imposing 
the initial conditions,we can define the displacement and velocity 
functions and, finally, return these two functions to the caller.
End of explanation
def bisect(f,val,x0,x1):
    h  = (x0+x1)/2.0
    fh = f(h)-val
    if abs(fh)<1e-8 : return h
    f0 = f(x0)-val
    if f0*fh > 0 : 
        return bisect(f, val, h, x1)
    else:
        return bisect(f, val, x0, h)
Explanation: An utility function
We need to find when
the spring yields
the velocity is zero
to individuate the three ranges of different behaviour
elastic
plastic
elastic, with permanent deformation.
We can use the simple and robust algorithm of bisection 
to find the roots for 
<center>$x_{el}(t)=x_y \text{ and } \dot{x}_{ep}(t)=0$.</center>
End of explanation
mass = 1000.   # kg
k    = 40000.  # N/m
zeta = 0.03    # damping ratio
fy   = 2500.   # N
print('Limit displacement Uy =', fy*1000/k, 'mm')
Explanation: The system parameters
End of explanation
damp = 2*zeta*sqrt(k*mass)
xy = fy/k    # m
Explanation: Derived quantities
The damping coefficient $c$ and the first yielding displacement, $x_y$.
End of explanation
t1 = 0.3     # s
w  = pi/t1   # rad/s
Po = 6000.   # N
Explanation: Load definition
Our load is a half-sine impulse
<center>$p(t)=\begin{cases}p_0\sin(\frac{\pi t}{t_1})&0\leq t\leq t_1,\\
    0&\text{otherwise.}\end{cases}$</center>
In our exercise
End of explanation
x0=0.0     # m
v0=0.0     # m/s
x_next, v_next = resp_elas(mass,damp,k, 0.0,Po,w, 0.0, x0,v0)
Explanation: The actual computations
Elastic, initial conditions, get system functions
End of explanation
t_yield = bisect(x_next, xy, 0.0, t1)
print(t_yield, x_next(t_yield)*k)
Explanation: Yielding time is
The time of yielding is found solving the equation $x_\text{next}(t) = x_y$
End of explanation
t_el = linspace( 0.0, t_yield, 201)
x_el = vectorize(x_next)(t_el)
v_el = vectorize(v_next)(t_el)
# ------------------------------
figure(0)
plot(t_el,x_el,
     (0,0.25),(xy,xy),'--b',
     (t_yield,t_yield),(0,0.0699),'--b')
title("$x_{el}(t)$")
xlabel("Time, s")
ylabel("Displacement, m")
# ------------------------------
figure(1)
plot(t_el,v_el)
title("$\dot x_{el}(t)$")
xlabel("Time, s")
ylabel("Velocity, m/s")
Explanation: Forced response in elastic range is
End of explanation
x0=x_next(t_yield)
v0=v_next(t_yield)
print(x0, v0)
Explanation: Preparing for EP response
First, the system state at $t_y$ is the initial condition for the EP response
End of explanation
cS = Po*cos(w*t_yield)
cC = Po*sin(w*t_yield)
print(Po*sin(w*0.55), cS*sin(w*(0.55-t_yield))+cC*cos(w*(0.55-t_yield)))
Explanation: now, the load must be expressed in function of a restarted time,
<center>
    $\tau=t-t_y\;\rightarrow\;t=\tau+t_y\;\rightarrow\;\sin(\omega t)=\sin(\omega\tau+\omega t_y)$
 </center> <center>
    $\rightarrow\;\sin(\omega t)=\sin(\omega\tau)\cos(\omega t_y)+\cos(\omega\tau)\sin(\omega t_y)$
</center>
End of explanation
x_next, v_next = resp_yield(mass, damp,    cC,cS,w, -fy, x0,v0)
Explanation: Now we generate the displacement and velocity functions for the yielded phase,
please note that the yielded spring still exerts a constant force $f_y$ on the mass,
and that this fact must be (and it is) taken into account.
End of explanation
t_y1 = linspace(t_yield, t1, 101)
x_y1 = vectorize(x_next)(t_y1-t_yield)
v_y1 = vectorize(v_next)(t_y1-t_yield)
figure(3)
plot(t_el,x_el, t_y1,x_y1,
     (0,0.25),(xy,xy),'--b',
     (t_yield,t_yield),(0,0.0699),'--b')
xlabel("Time, s")
ylabel("Displacement, m")
# ------------------------------
figure(4)
plot(t_el, v_el,  t_y1, v_y1)
xlabel("Time, s")
ylabel("Velocity, m/s")
Explanation: At this point I must confess that I have already peeked the numerical solution,
hence I know that the velocity at $t=t_1$ is still greater than 0 and I know that the current
solution is valid in the interval $t_y\leq t\leq t_1$.
End of explanation
x0 = x_next(t1-t_yield)
v0 = v_next(t1-t_yield)
print(x0, v0)
x_next, v_next = resp_yield(mass, damp, 0, 0, w, -fy, x0, v0)
t2 = t1 + bisect( v_next, 0.0, 0, 0.3)
print(t2)
t_y2 = linspace( t1, t2, 101)
x_y2 = vectorize(x_next)(t_y2-t1)
v_y2 = vectorize(v_next)(t_y2-t1)
print(x_next(t2-t1))
figure(5)
plot(t_el,x_el, t_y1,x_y1, t_y2, x_y2,
     (0,0.25),(xy,xy),'--b',
     (t_yield,t_yield),(0,0.0699),'--b')
xlabel("Time, s")
ylabel("Displacement, m")
# ------------------------------
figure(6)
plot(t_el, v_el,  t_y1, v_y1, t_y2, v_y2)
xlabel("Time, s")
ylabel("Velocity, m/s")
Explanation: In the next phase,  still it is $\dot x> 0$ so that
the spring is still yielding,
but now $p(t)=0$, so we must compute two new state functions,
starting as usual from the initial conditions (note that the yielding
force is still applied)
End of explanation
x0 = x_next(t2-t1) ; v0 = 0.0
x_next, v_next = resp_elas(mass,damp,k, 0.0,0.0,w, k*x0-fy, x0,v0)
t_e2 = linspace(t2,4.0,201)
x_e2 = vectorize(x_next)(t_e2-t2)
v_e2 = vectorize(v_next)(t_e2-t2)
Explanation: Elastic unloading
The only point worth commenting is the constant force
that we apply to our system.
The force-displacement relationship for an EP spring is 
<center>$f_\text{E} = k(x-x_\text{pl})= k x - k (x_\text{max}-x_y)$</center>
taking the negative, constant part of the last expression
into the right member of the equation of equilibrium we have a constant
term, as follows
End of explanation
# ------------------------------
figure(7) ;
plot(t_el, x_el, '-b',
     t_y1, x_y1, '-r',
     t_y2, x_y2, '-r',
     t_e2, x_e2, '-b',
     (0.6, 4.0), (x0-xy, x0-xy), '--y')
title("In blue: elastic phases.\n"+
      "In red: yielding phases.\n"+
      "Dashed: permanent plastic deformation.")
xlabel("Time, s")
ylabel("Displacement, m")
Explanation: now we are ready to plot the whole response
End of explanation
def make_p(p0,t1):
  make_p(p0,t1) returns a 1/2 sine impulse load function, p(t)
  def p(t):
    ""
    if t<t1:
      return p0*sin(t*pi/t1)
    else:
      return 0.0
  return p
Explanation: Numerical solution
first, we need the load function
End of explanation
def make_kt(k,fy):
  "make_kt(k,fy) returns a function kt(u,v,up) returning kt, up"
  def kt(u,v,up):
    f=k*(u-up)
    if (-fy)<f<fy:               return k,up
    if fy<=f    and v>0:  up=u-uy;return 0,up
    if fy<=f    and v<=0: up=u-uy;return k,up
    if f<=(-fy) and v<0:  up=u+uy;return 0,up
    else:                 up=u+uy;return k,up
  return kt
Explanation: and also a function that, given the displacement,
the velocity and the total plastic deformation, returns
the stiffness and the new p.d.; this function is defined
in terms of the initial stiffness and the yielding load
End of explanation
# Exercise from lesson 04
# 
mass = 1000.00       # kilograms
k  =  40000.00       # Newtons per metre
zeta  =   0.03       # zeta is the damping ratio
fy =   2500.00       # yelding force, Newtons
t1 =      0.30       # half-sine impulse duration, seconds
p0 =   6000.00       # half-sine impulse peak value, Newtons
uy =      fy/k       # yelding displacement, metres
Explanation: Problem data
End of explanation
# using the above constants, define the loading function
p=make_p(p0,t1)
# the following function, given the final displacement, the final
# velocity and the initial plastic deformation returns a) the tangent
# stiffness b) the final plastic deformation
kt=make_kt(k,fy) 
# we need the damping coefficient "c", to compute its value from the
# damping ratio we must first compute the undamped natural frequency 
wn=sqrt(k/mass)       # natural frequency of the undamped system
damp=2*mass*wn*zeta         # the damping coefficient
# the time step
h=0.005
# required duration for the response
t_end = 4.0
# the number of time steps to arrive at t_end
nsteps=int((t_end+h/100)/h)+1
# the maximum number of iterations in the Newton-Raphson procedure
maxiters = 30
# using the constant acceleration algorithm
# below we define the relevant algorithmic constants
gamma=0.5
beta=1./4.
gb=gamma/beta
a=mass/(beta*h)+damp*gb
b=0.5*mass/beta+h*damp*(0.5*gb-1.0)
Explanation: Initialize the algorithm
compute the functions that return the load and the tangent sstiffness + plastic deformation
compute the damping constant
for a given time step, compute all the relevant algorithmic constants, with $\gamma=\frac12$ and
   $\beta=\frac14$
End of explanation
t0=0.0
u0=0.0
up=0.0
v0=0.0
p0=p(t0)
(k0, up)=kt(u0,v0,up)
a0=(p0-damp*v0-k0*(u0-up))/mass
time = []; disp = []
Explanation: System state initialization
and a bit more, in species we create two empty vectors to hold the computation results
End of explanation
for i in range(nsteps):
    time.append(t0); disp.append(u0)
    # advance time, next external load value, etc
    t1 = t0 + h
    p1 = p(t1)
    Dp = p1 - p0
    Dp_= Dp + a*v0 + b*a0
    k_ = k0 + gb*damp/h + mass/(beta*h*h)
    # we prepare the machinery for the modified Newton-Raphson
    # algorithm.  if we have no state change in the time step, then the
    # N-R algorithm is equivalent to the standard procedure
    u_init=u0; v_init=v0 # initial state
    f_spring=k*(u0-up)   # the force in the spring
    DR=Dp_               # the unbalanced force, initially equal to the
                         # external load increment
    for j in range(maxiters):
        Du=DR/k_           # the disp increment according to the initial stiffness
        u_next = u_init + Du
        v_next = v_init + gb*Du/h - gb*v_init + h*(1.0-0.5*gb)*a0
                           # we are interested in the total plastic elongation
        oops,up=kt(u_next,v_next,up)
                           # because we need the spring force at the end
                           # of the time step
        f_spring_next=k*(u_next-up)
                           # so that we can compute the fraction of the
                           # incremental force that's equilibrated at the
                           # end of the time step
        df=f_spring_next-f_spring+(k_-k0)*Du
                           # and finally the incremental forces unbalanced
                           # at the end of the time step
        DR=DR-df
                           # finish updating the system state
        u_init=u_next; v_init=v_next; f_spring=f_spring_next
                           # if the unbalanced load is small enough (the
                           # criteria used in practical programs are
                           # energy based) exit the loop - note that we
                           # have no plasticization/unloading DR==0 at the
                           # end of the first iteration
        if abs(DR)<fy*1E-6: break
    # now the load increment is balanced by the spring force and
    # increments in inertial and damping forces, we need to compute the
    # full state at the end of the time step, and to change all
    # denominations to reflect the fact that we are starting a new time step.
    Du=u_init-u0
    Dv=gamma*Du/(beta*h)-gamma*v0/beta+h*(1.0-0.5*gamma/beta)*a0
    u1=u0+Du ; v1=v0+Dv
    k1,up=kt(u1,v1,up)
    a1=(p(t1)-damp*v1-k*(u1-up))/mass
    t0=t1; v0=v1; u0=u1 ; a0=a1 ; k0=k1 ; p0=p1
Explanation: Iteration
We iterate over time and, if there is a state change, over the single time step to equilibrate
the unbalanced loadings
End of explanation
figure(8)
plot(time[::4],disp[::4],'xr')
plot(t_el, x_el, '-b',
     t_y1, x_y1, '-r',
     t_y2, x_y2, '-r',
     t_e2, x_e2, '-b',
     (0.6, 4.0), (x0-xy, x0-xy), '--y')
title("Continuous line: exact response.\n"+
      "Red crosses: constant acceleration + MNR.\n")
xlabel("Time, s")
ylabel("Displacement, m");
Explanation: Plotting our results
we plot red crosses for the numericaly computed response and a continuous line for the results
of the analytical integration of the equation of motion.
End of explanation
<END_TASK>