0
$\begingroup$

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$

enter image description here

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:

enter image description here

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()
$\endgroup$
2
$\begingroup$

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:

enter image description here

$\endgroup$
  • $\begingroup$ Thank you very much for taking the time to look into my stupid code. This will help tremendously. $\endgroup$ – user32882 Jul 28 '17 at 11:31

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.