Exercise_13 : Waves

Problem 6.12 & 6.16

冉峰 弘毅班 2014301020064

Abstract

• This article is about using numerical method to solve wave function, especially for a wave spreading in a string. Further more, it will discuss the semirealistic case and realistic case,just like piano strings, with which whether we considering friction.
  

Background

• In physics, a wave is an oscillation accompanied by a transfer of energy that travels through a medium (space or mass). Frequency refers to the addition of time. Wave motion transfers energy from one point to another, which displace particles of the transmission medium–that is, with little or no associated mass transport. Waves consist, instead, of oscillations or vibrations (of a physical quantity), around almost fixed locations.
  

The Main Body

Ideal case

• The central equation of wave motion is:
  $$\begin{eqnarray*} \frac{\partial^2y}{\partial t^2}&=&c^2\frac{\partial^2y}{\partial x^2} \end{eqnarray*}$$
  
  From the figure, if we assume the string displacement is small, so that the angles $\theta_i$ are also small,we can write the equation of motion of segment $i$ as :
  $$\begin{eqnarray*} (\mu\Delta x)\frac{\partial^2y_i}{\partial t^2}&=Tsin\theta_{i+1}-Tsin\theta_{i} \end{eqnarray*}$$
  Then we use a finite difference approximation for the angles:
  $$\begin{eqnarray*} sin\theta_{i}\approx\frac{y_i-y_{i-1}}{\Delta x} \end{eqnarray*}$$
  We obtain:
  $$\begin{eqnarray*} \frac{\partial^2y_i}{\partial t^2}\approx\frac{T}{\mu}\frac{y_{i+1}-2y_i+y_{i-1}}{(\Delta x ^2)}\approx(\frac{T}{\rho})\frac{\partial^2y}{\partial x^2} \end{eqnarray*}$$
  By taking time and x into finite setps $x=i\Delta x$ and $t=n\Delta t$. In the numerical approach, we rewrite the equation into this form:
  $$\begin{eqnarray*} \frac{y(i,n+1)+y(i,n-1)-2y(i,n)}{(\Delta t)^2}\approx c^2[\frac{y(i,n+1)+y(i,n-1)-2y(i,n)}{(\Delta x)^2}] \end{eqnarray*}$$
  We can express $y(i,n+1)$ in terms of $y$ at previous time steps, with result:
  $$y(i,n+1)=2[1-r^2]y(i,n)-y(i,n-1)+r^2[y(i+1,n)+y(i-1,n)]$$

Gaussian pluck (semirealistic)

• That is, we have taken the initial string profile to be :
  $$y_0(x)=exp[-k(x-x_0)^2]$$
  

Realistic case

• In realistic case, if we consider the friction, the equation of vertical motion becomes:
  $$\begin{equation*} \frac{\partial ^2 y}{\partial t^2} = c^2(\frac{\partial ^2 y}{\partial x^2} - \epsilon L^2\frac{\partial ^4 y}{\partial x^4}) \end{equation*}$$
  Where $\epsilon$ is a dimisionless stiffness parameter.
  Using the same approach that we used to obtain that second partial derivation in C hapter 5. The result is:
  $$\frac{\partial^4y}{\partial x^4}=\frac{y(i+2,n)-4y(i+1,n)+6y(i,n)-4y(i-1,n)+y(i-2,n)}{(\Delta x)^4}$$
  Also, in the case of realistic string, the numerical equation becomes:
  $$y(i,n+1) = [2-2r^2 - 6\epsilon r^2M^2]y(i,n) - y(i,n-1) \\ +r^2[1 + 4\epsilon M^2][y(i+1,n)+y(i-1,n)] -\\ \epsilon r^2 M^2[y(i+2,n)-y(i-2,n)]$$
  Where $M=L/\Delta x$ is the number of spatial units along the string.

Code

import matplotlib.pyplot as plt
from matplotlib import animation
import numpy as np
from scipy.fftpack import rfft, irfft, fftfreq
class Guass_seidel:
    def __init__(self,L,dx,T,e):
        self.L = L
        self.dx = dx
        self.r = 0.25
        self.c = 300
        self.T = T
        self.dt = self.dx/(4*self.c)
        self.e = e 
        pass
    def initialization(self):
        self.string = []
        for i in range(int(self.T/self.dt)):
            self.string.append([])
            for j in range(int(self.L/self.dx)):
                self.string[i].append(0)
        for l in range(1,int(int(self.L/self.dx) - 1)):
            self.string[0][l] = np.e**(-0.1*(l - 30)**2) 
        self.string[0][1] = 0
        self.string[0][-2] = 0
        self.string[0][0] =  - self.string[0][2]
        self.string[0][-1] = - self.string[0][-3]
        self.string[1] = self.string[0][:]
        #plt.plot(self.string[1])
    def calculation(self):
        M = int(self.L/self.dx)
        for i in range(2,int(self.T/self.dt)):
            for j in range(2,M-2):
                self.string[i][j] = (2 - 2*self.r**2 - 6*self.e*self.r**2*M**2)*self.string[i - 1][j] - self.string[i - 2][j] + self.r**2*(1 + 4*self.e*M**2)*(self.string[i-1][j+1]+self.string[i-1][j-1]) - self.e*self.r**2*M**2*(self.string[i-1][j+2]+self.string[i-1][j-2])
                #self.string[i][j] = 2*(1-self.r**2)*self.string[i - 1][j] - self.string[i - 2][j] + self.r**2*(self.string[i - 1][j + 1] +self.string[i - 1][j - 1])
                self.string[i][0] = -self.string[i][2]
                self.string[i][-1] = -self.string[i][-3]
        #plt.plot(self.string[-1])
        return self.string
    def fft(self):
        t_plot = np.arange(0,self.T,self.dt)
        amplitude_record = []
        for i in range(int(self.T/self.dt)):
            amplitude_record.append(self.string[i][20])
        freq = fftfreq(len(amplitude_record), d=self.dt)
        freq = np.array(abs(freq))
        f_signal = rfft(amplitude_record)
        f_signal = np.array(f_signal**2)
        #plt.subplot(122)
        plt.title('Power spectra')
        plt.ylabel('Power (arbitrary units)')
        plt.xlabel('Frequency (Hz)')
        plt.xlim(2000,8000)
        plt.plot(freq,f_signal,label = 'epsilon = '+str(self.e))
        return 0
A = Guass_seidel(1,0.01,0.05,0)
A.initialization()
A.calculation()
A.fft()        
# ---- animation ---------------------------------
fig = plt.figure(figsize = (8,6))
ax = plt.axes(xlim=(0, 1),ylim = (-1,1))
line, = ax.plot([], [], 'k')
def init():
    line.set_data([], [])
    return line,    
def animate(i):
    x_plot = np.arange(0.01,0.99,0.01)
    y_plot = []
    for j in range(1, int(A.L/A.dx) - 1):
        y_plot.append(data_record[i][j])
    line.set_data(x_plot,y_plot)
    return line,
anim = animation.FuncAnimation(fig, animate, init_func=init, frames=100, interval=50, blit=True)
anim.save('chapter6_string_guass_e.gif', fps=20, writer='Feng_Chen')
plt.show()

Result

• First we can investigate the wave on a string, and chose the condition used in figure 6.2 that $x_0=0.3m$ and $k=1000m^{-2}$ and string runs from $x=0$ to 1m.
  
  Then we want to find the result that give two Gaussian pluck, one is at $x_0=0.3$, and the other at $x_0=0.7$.
  
  From the results above, we can draw the conclusion that when there are two Gaussian wave packets located at different places on the string, the wave packets may then propagate and collide but the wave packets are unaffected by the collisions.

• Also I want to find the result that frequency spectrum of waves on a Gaussian-excited string.
  Suppose the initial wavepacket is $$y_0(x)=exp[-1000\times(x-x_{excite})]$$
  And the total length of the string is 1 unit.
  When $x_{excite}=0.5$ and $x_{observe}=0.05$
  We have :
  
  We can polt the power spectrum figure:
  
  When the point deviation $x_{excite}=0.5$, we chose $x_{excite}=0.45$, to find the power spectrum figure.
  
  Because we have the possible frequencies as $f=mc/2L=150mHz$. This explains why the peaks in the spectral analysis in the above figures occur at regularly spaced frequencies. Each of the peaks correspond to one value of interger m. But some frequancies are missing and this can be traced to the operation of Fourier Transformation. So when $f$ is the times of 150, there will be a wave crest.

• When the intial wave packet was a real one.
  
  We can polt the power spectrum figure:
  
  when When the point deviation from the center of spring,
  We can polt the power spectrum figure:
  
  Because of the loss of symmetry, there be more wave crest.

Conclusion

• In realistic case, the friction of the string will cause a frequency shift especially for the high frequency. Also, the frequencies which excited in a guitar are always low frequencies.

Thanks For

• The book of Computational Physics, chapter 6
• Wu Yuqiao
• Wikipedia