<SYSTEM_TASK:>
Given the following text description, write Python code to implement the functionality described below step by step
<END_TASK>
<USER_TASK:>
Description:
 Constant Acceleration
Define the Dynamical System
Step1:  Define the Loading
Our load is
$$p(t) = 1 \text{N} \times \begin{cases}\sin^2\frac{\omega_n}2t&0 \le t \le 5\0&\text{otherwise}\end{cases}$$
Step2:  Numerical Constants
We want to compute, for each step
$$\Delta x = \frac{\Delta p + a^\star a_0 + v^\star v_0}{k^\star}$$
where we see the constants
$$a^\star = 2m, \quad v^\star = 2c + \frac{4m}{h},$$
$$ k^\star = k + \frac{2c}{h} + \frac{4m}{h^2}.$$
Step3:  Vectorize the time and the load
We want to compute the response up to 8 seconds
Step4:  Integration
Prepare containers
write initial conditions
loop on load and load increments
save initial conditions (x0, v0, a0),
compute the _corrected_load increment,
compute dx and dy,
compute the vaues of x0 and y0 for the beginning of the next step
compute a0 at the beginning of the next step.
save the results for the last step (they are usually saved at the beginning of the loop, but here we need to special case, why?) and vectorize the whole of the results
Step5:  Results
The response
Step6:  Comparison
Using black magic it is possible to divine the analytical expression of the response during the forced phase,
$$x(t) = \frac1{2k}\left(\left(\frac{1-2\zeta^2}{2\zeta\sqrt{1-\zeta^2}}\sin\omega_Dt-\cos\omega_Dt\right)\exp(-\zeta\omega_nt)+1-\frac{1}{2\zeta}\sin\omega_nt\right), \qquad 0\le t \le 5.$$
and hence plot a comparison within the exact response and the (downsampled) numerical approximation, in the range of validity of the exact response.
Step7:  Eventually we plot the difference between exact and approximate response,
<ASSISTANT_TASK:>
Python Code:
m = 1.00
k = 4*pi*pi
wn = 2*pi
T = 1.0
z = 0.02
wd = wn*sqrt(1-z*z)
c = 2*z*wn*m
Explanation: Constant Acceleration
Define the Dynamical System
End of explanation
NSTEPS = 200 # steps per second
h = 1.0 / NSTEPS
def load(t):
    return np.where(t<0, 0, np.where(t<5, sin(0.5*wn*t)**2, 0))
t = np.linspace(-1, 6, 7*NSTEPS+1)
plt.plot(t, load(t))
plt.ylim((-0.05, 1.05));
Explanation: Define the Loading
Our load is
$$p(t) = 1 \text{N} \times \begin{cases}\sin^2\frac{\omega_n}2t&0 \le t \le 5\0&\text{otherwise}\end{cases}$$
End of explanation
kstar = k + 2*c/h + 4*m/h/h
astar = 2*m
vstar = 2*c + 4*m/h
Explanation: Numerical Constants
We want to compute, for each step
$$\Delta x = \frac{\Delta p + a^\star a_0 + v^\star v_0}{k^\star}$$
where we see the constants
$$a^\star = 2m, \quad v^\star = 2c + \frac{4m}{h},$$
$$ k^\star = k + \frac{2c}{h} + \frac{4m}{h^2}.$$
End of explanation
t  = np.linspace(0, 8+h, NSTEPS*8+1)
P = load(t) 
Explanation: Vectorize the time and the load
We want to compute the response up to 8 seconds
End of explanation
x, v, a = [], [], []
x0, v0 = 0.0, 0.0
a0 =(P[0]-k*x0-c*v0)/m
for p0, p1 in zip(P[:-1], P[+1:]):
    x.append(x0), v.append(v0), a.append(a0)
    dx = ((p1-p0) + astar*a0 + vstar*v0)/kstar
    dv = 2*(dx/h-v0)
    x0, v0 = x0+dx, v0+dv
    a0 = (p1 - k*x0 - c*v0)/m
x.append(x0), v.append(v0), a.append(a0)
x, v = np.array(x), np.array(v)
Explanation: Integration
Prepare containers
write initial conditions
loop on load and load increments
save initial conditions (x0, v0, a0),
compute the _corrected_load increment,
compute dx and dy,
compute the vaues of x0 and y0 for the beginning of the next step
compute a0 at the beginning of the next step.
save the results for the last step (they are usually saved at the beginning of the loop, but here we need to special case, why?) and vectorize the whole of the results
End of explanation
plt.plot(t, x);
Explanation: Results
The response
End of explanation
xe = (((1 - 2*z**2)*sin(wd*t) / (2*z*sqrt(1-z*z)) - cos(wd*t)) *exp(-z*wn*t) + 1. - sin(wn*t)/(2*z))/2/k
plt.plot(t[:1001], xe[:1001], lw=2)
plt.plot(t[:1001:10], x[:1001:10], 'ko');
Explanation: Comparison
Using black magic it is possible to divine the analytical expression of the response during the forced phase,
$$x(t) = \frac1{2k}\left(\left(\frac{1-2\zeta^2}{2\zeta\sqrt{1-\zeta^2}}\sin\omega_Dt-\cos\omega_Dt\right)\exp(-\zeta\omega_nt)+1-\frac{1}{2\zeta}\sin\omega_nt\right), \qquad 0\le t \le 5.$$
and hence plot a comparison within the exact response and the (downsampled) numerical approximation, in the range of validity of the exact response.
End of explanation
plt.plot(t[:1001], x[:1001]-xe[:1001]);
Explanation: Eventually we plot the difference between exact and approximate response,
End of explanation
<END_TASK>