<SYSTEM_TASK:>
Given the following text description, write Python code to implement the functionality described below step by step
<END_TASK>
<USER_TASK:>
Description:
 Piecewise Exact Integration
The Dynamical System
We want to study a damped SDOF system, so characterized
Step1:  The excitation is given by a force such that the static displacement is 5 mm, modulated by a sine in resonance with the dynamic sistem, i.e., $\omega=\omega_n$.
Step2:  For such a system, we know exactly the response.  The particular integral is 
$$\xi(t)=-\frac{\cos\omega t}{2\zeta}$$
and imposing rest initial conditions it is
$$x(t) = \frac{\Delta_{st}}{2\zeta} ((\frac{\zeta}{\sqrt{1=\zeta^2}}\sin\omega_Dt + \cos\omega_Dt)\exp(-\zeta\omega t) - \cos\omega t),\qquad \omega=\omega_n.
$$
Step3:  Numerical integration
Now we prepare for the numerical integration, first the constants that represent the homogeneous response at end of step
Step4:  then the actual step computations, where in terms of the initial state vector and the load variation the final state is derived.
Step5:  With those pieces in place, we can define a function that, for a given number of steps per period computes the response on the interval $0 \le t \le 2.0$.
Step6:  Let's compute the responses for different numbers of steps, and store them away too...
Step7:  Eventually we can plot the numerical responses along with the exact response
Step8:  But... there are only two numerical curves and I've plotted three of them.
Let's plot the difference between the exact response and the response computed at 16 samples per period...
<ASSISTANT_TASK:>
Python Code:
T=1.0                  # Natural period of the oscillator
w=2*pi                 # circular frequency of the oscillator
m=1000.0               # oscillator's mass, in kg
k=m*w*w                # oscillator stifness, in N/m
z=0.05                 # damping ratio over critical
c=2*z*m*w              # damping
wd=w*sqrt(1-z*z)       # damped circular frequency
ratio=sqrt(1-z*z)      # ratio damped/undamped frequencies
Explanation: Piecewise Exact Integration
The Dynamical System
We want to study a damped SDOF system, so characterized
End of explanation
D=0.005                # static displacement, 5mm
P=D*k                  # force amplitude
Explanation: The excitation is given by a force such that the static displacement is 5 mm, modulated by a sine in resonance with the dynamic sistem, i.e., $\omega=\omega_n$.
End of explanation
def exact(t):
    return D*((z*sin(wd*t)/ratio+cos(wd*t))*exp(-z*w*t)-cos(w*t))/(2*z)
t = np.linspace(0.0, 2.0, 1001)
plt.plot(t, exact(t));
Explanation: For such a system, we know exactly the response.  The particular integral is 
$$\xi(t)=-\frac{\cos\omega t}{2\zeta}$$
and imposing rest initial conditions it is
$$x(t) = \frac{\Delta_{st}}{2\zeta} ((\frac{\zeta}{\sqrt{1=\zeta^2}}\sin\omega_Dt + \cos\omega_Dt)\exp(-\zeta\omega t) - \cos\omega t),\qquad \omega=\omega_n.
$$
End of explanation
def initstep(h):
    cdh=cos(wd*h)*exp(-z*w*h)
    sdh=sin(wd*h)*exp(-z*w*h)
    return cdh, sdh
Explanation: Numerical integration
Now we prepare for the numerical integration, first the constants that represent the homogeneous response at end of step
End of explanation
def step(x0,v0,p0,p1,h,cdh,sdh):
    dst=p0/k
    ddst=(p1-p0)/k
    B = x0 - dst + ((2*z)/w)*(ddst/h)
    A = (v0 + z*w*B - ddst/h)/wd
    x1 = A*sdh + B*cdh + dst + ddst - ddst/h * 2*z/w
    v1 = A*(wd*cdh-z*w*sdh) - B*(z*w*cdh+wd*sdh) + ddst/h
    return x1, v1
Explanation: then the actual step computations, where in terms of the initial state vector and the load variation the final state is derived.
End of explanation
def resp(nstep):
    T = np.linspace(0.0, 2.0, 2*nstep + 1)
    X = np.zeros(2*nstep + 1)
    
    h=1./float(nstep)
    cdh, sdh = initstep(h)
    x1=0. ; v1=0. ; p1=0
    
    for i, t in enumerate(T):
        X[i] = x1
        x0=x1 ; v0=v1 ; p0=p1 ; p1=P*sin(w*(t+h))
        x1,v1=step(x0,v0,p0,p1,h, cdh, sdh)
    return T, X
Explanation: With those pieces in place, we can define a function that, for a given number of steps per period computes the response on the interval $0 \le t \le 2.0$.
End of explanation
t_x = {n:resp(n) for n in (4, 8, 16)}
Explanation: Let's compute the responses for different numbers of steps, and store them away too...
End of explanation
plt.plot(t, exact(t), label='Exact', lw=3)
plt.plot(*t_x[4], *t_x[8], *t_x[16])
plt.grid()
plt.legend(loc=3)
plt.xlabel('Time t/s')
plt.ylabel('Displacement x/m');
Explanation: Eventually we can plot the numerical responses along with the exact response
End of explanation
t16, x16 = t_x[16]
plt.plot(t16, exact(t16)-x16)
Explanation: But... there are only two numerical curves and I've plotted three of them.
Let's plot the difference between the exact response and the response computed at 16 samples per period...
End of explanation
<END_TASK>