# Trying to plot 1D wave equation for benchmarking

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()


There are two problems with your code:

1) You are mixing $i$ and $n$ in the following line:

SUM = SUM+Ai*cos((i-0.5)*pi*xi/L)*sin((n-0.5)*pi*ti/(L/c))


It should be:

SUM = SUM+Ai*cos((i-0.5)*pi*xi/L)*sin((i-0.5)*pi*ti/(L/c))


2) You have to normalize the rod displacement. The coefficients are then:

    Ai = 2*sin((i-0.5)*a*pi/L)/(pi*(i-0.5)**2*pi*a/L)


The result of these changes is: • Thank you very much for taking the time to look into my stupid code. This will help tremendously. Jul 28 '17 at 11:31