I am trying to plot a reference solution for the 1D wave equation using python.
The above link states the following: For a rod fixed at the right end and free at the left end and subjected to a hammer blow at the left end at $t = 0$ s where the hammer blow induces a constant velocity $V_0$ to occur from $x<a<L$, the displacement is given as follows:
$$A_n = \frac{2V_0 a \sin((n-\frac 1 2)\pi a /L)}{v \pi (n-\frac 1 2)^2 \pi a/L}$$
and
$$\psi(x,t) = \sum_{n=1,\infty}\cos((n-\frac 1 2)\pi x/L)\sin((n-\frac 1 2)\pi t/\tau) $$
where $\tau = L/c$, $\psi$ is the displacement, $v = \sqrt{E/\rho}$ and $a$ is the distance where the initial velocity $V_0$ occurs.
It appears the author has used $c$ in $\tau = L/c$ instead of using $v$ for the wave speed? This is confusing to me but my intuition says that $c = v = \sqrt{E/\rho}$
For normalized rod displacement $\hat \psi(x,t) = (c / V_0 a) \psi(x,t)$ calculated from the preceding equations using the first 100 normal modes (i.e $n = 1,100$) and choosing $a/L = 0.1$, this is what we get for $t/\tau = 1.28$
Now coding in python:
import matplotlib.pyplot as plt
from numpy import *
plt.ion()
n = 100
L = 25
a = 2.5
E = 100
rho = 1
V0 = 3
c = (E/rho)**0.5
x = linspace(0, L, 100)
u = []
ti = 1.28*L/c
for xi in x:
SUM = 0
for i in range(1,n+1,1):
Ai = V0*a*2*sin((i-0.5)*a*pi/L)/(c*pi*(i-0.5)**2*pi*a/L)
SUM = SUM+Ai*cos((i-0.5)*pi*xi/L)*sin((n-0.5)*pi*ti/(L/c))
u.append(SUM)
u = [x * (c/V0*a) for x in u]
plt.plot(x/L,u)
plt.show()
Yields the following plot:
Which really has nothing to do with the correct plot above.
I have checked and rechecked the equations. I am not sure where I am going wrong and hope it is not an error from what the author has provided? The most important thing is the final result so I hope I am implementing it correctly? Some help please.
EDIT
Based on the answer below, here is the working code.
import matplotlib.pyplot as plt
from numpy import *
plt.ion()
n = 100
L = 25
a = 2.5
E = 100
rho = 1
V0 = 3
c = (E/rho)**0.5
x = linspace(0, L, 100)
u = []
ti = 1.28*L/c
for xi in x:
SUM = 0
for i in range(1,n+1,1):
Ai = V0*a*2*sin((i-0.5)*a*pi/L)/(c*pi*(i-0.5)**2*pi*a/L)
SUM = SUM+Ai*cos((i-0.5)*pi*xi/L)*sin((i-0.5)*pi*ti/(L/c))
u.append(SUM)
u = [x * (c/V0/a) for x in u]
plt.plot(x/L,u)
plt.show()
