[Scientific / technical calculation by Python] Numerical solution of one-dimensional and two-dimensional wave equations by FTCS method (explicit method), hyperbolic partial differential equations

Introduction

The vibration transmitted through the substance / space itself (medium) that spreads in space is generally called a wave or wave [1]. Elastic waves and fluid motions correspond to mechanical waves, and electromagnetic waves correspond to electromagnetic waves. The equation satisfied by the physical quantity $ u (x, y, z, t) $ oscillating in the medium is the wave equation.

** In this paper, Python is used for 2D wave equation (hyperbolic partial differential equation) ** $\frac{1}{c^2}\frac{\partial^2 u}{\partial t^2} = \frac{\partial^2 u}{\partial x^2} + \frac{\partial^2 u}{\partial y^2} $

Solving ** by the FTCS method (Forward in Time and Centered Difference method) [Addendum], which is an explicit difference scheme ** (Addition: ** Also included in the case of one dimension. See content (2) **)

Note that in deriving this equation in mechanics, the assumption of minute vibration and the uniformity and isotropic nature of the medium that transmits the wave are assumed (phase velocity c is a constant).

Contents

(1) Two-dimensional wave equation

Solve the two-dimensional wave equation by the FTCS method, which is an explicit difference scheme. The calculation conditions are as follows.

• Computation area: Square with Lx = 1, Ly = 1 In 41-point increments in both the x direction (Nx) and the y direction (Ny) (Nx = Ny = 41)

• The maximum value (t_max) of the calculation time is 0.5.

• Time step (delta_t): The value that takes numerical stability into consideration is halved.

• Wave phase velocity: c = 1

• Boundary condition: All four sides are fixed ends (does not vibrate and does not move)

• Initial condition 1: Gaussian distribution function is installed at the center of the square at time 0. This "pulls" the wave (see Figure 1).

• Initial condition 2: $ \ frac {\ partial u} {\ partial t} = 0 \ (t = 0) $ (initial velocity is 0)

In addition, the calculation requires the information of the wave $ u $ in the t = 1 step. In this simulation, the initial condition 2 was evaluated using the result obtained by center difference [4].

Figure 1. As an initial condition, Gaussian is installed in the center of the square.

(2) One-dimensional wave equation

$\frac{1}{c^2}\frac{\partial^2 u}{\partial t^2} = \frac{\partial^2 u}{\partial x^2} $ Solve with the FTCS method.

The boundary condition is $ u = 0 $ (fixed end) at the end point. The initial condition is $ \ frac {\ partial u} {\ partial t} = 0 \ (t = 0) $ (initial velocity is 0).

The initial conditions for the waveform are as shown in the figure below.

Figure 2. Initial waveform

Calculation code (1) Two-dimensional wave equation

#%matplotlib nbagg
"""
2D Wave Equation
Uniform medium:Phase velocity is constant in time and space
FTCS method
"""

import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
from matplotlib.animation import ArtistAnimation #Import methods for creating animations


c=1
print ("c=",c)
Lx = 1 #length[m]
Ly = 1
t_max = 0.5

Nx = 41
Ny =41

delta_x = Lx/Nx
delta_y = Ly/Ny
delta_t=(1/(1/delta_x+1/delta_y))/c*0.5 #Time step width. Take it small, though.

print("delta_t=",delta_t)

Nt = int(t_max/delta_t)
 

print("Nt=",Nt)

stx=delta_x/delta_t
sty=delta_y/delta_t



u =  np.empty((Nx,Ny,Nt),dtype=float)

#boundary condition
u[0,:,:], u[-1,:,:] , u[:,0,:] , u[:,-1,:] = 0, 0, 0, 0

#Initial conditions
for i in range (Nx) :
    for ii in range(Ny):
        u[i,ii,0]  = (np.exp(-((i*delta_x-Lx/2)**2+(ii*delta_y-Ly/2)**2)/0.01))  #Initial conditions. Installation of Gaussian distribution function.

#Hado u[i,ii,1]Calculation to make. It is obtained by centering the initial condition 2.
for i in range(1,Nx-1): 
    for ii in range(1,Ny-1):
        u[i,ii,1]= u[i,ii,0]+0.5*((delta_t/delta_x)**2)*(c**2)*(u[i+1,ii, 0]+u[i-1, ii,0]-2*u[i,ii,0])\
+0.5*((delta_t/delta_y)**2)*(c**2)*(u[i,ii+1, 0]+u[i,ii-1, 0]-2*u[i,ii,0])


# #Loop for time evolution of waves
for j in range(1,Nt-1):
    for i in range(1,Nx-1):
        for ii in range(1,Ny-1):
            u[i,ii, j+1] = 2*u[i,ii,j]-u[i,ii,j-1]+((c/stx)**2)* (u[i+1,ii,j]+u[i-1,ii,j]\
-2*u[i,ii,j])+((c/sty)**2)* (u[i,ii+1,j]+u[i,ii-1,j]-2*u[i,ii,j]) 
       

The following code is used to visualize the obtained wave function u [x, y, t]. The wireframe display [3] is plotted and the result is animated.

#Program for visualization

%matplotlib nbagg 
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
import matplotlib.animation as animation #Import methods for creating animations

fig = plt.figure(figsize = (8, 6))
ax = fig.gca(projection='3d')

x=list(range(Nx))
y=list(range(Ny))
X, Y = np.meshgrid(x,y)

def update(i,u,fig_title):
    if i != 0:
        ax.clear()
    ax.view_init(elev=60., azim=60.) #Angle setting
    #ax.plot_surface(X,Y,u[X,Y,i],rstride=1,cstride=1,cmap='Greys',linewidth=0.3) #Surface plot
    ax.plot_wireframe(X,Y,u[X,Y,i],rstride=1,cstride=1,color='blue',linewidth=0.3) #Wireframe display
    
    ax.set_title(fig_title + 'i=' + str(i))
    ax.set_zlim(-0.4,0.4) #Specifying the drawing range of z
    ax.set_xlabel('X') #label
    ax.set_ylabel('Y')
    ax.set_zlabel('U')


    
    return ax


ani = animation.FuncAnimation(fig,update,fargs = (u,'Wave motion: time step='), interval =1, frames = Nt,blit=False)
fig.show()

ani.save("output.gif", writer="imagemagick")

Result (1)

Shows the animated result. It can be seen that the wave generated at the center position propagates around.

Calculation code (2): One-dimensional wave equation

%matplotlib nbagg
"""
1D　Wave Equation
FTCS method

"""

import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
from matplotlib.animation import ArtistAnimation #Import methods for creating animations


T=40 #tension[N]
rho=0.01 #density[Kg/m]
c = np.sqrt(T/rho) #Phase velocity

print ("c=",c)
Lx = 1 #length[m]
t_max = 0.1

Nx = 100

delta_x = Lx/Nx

delta_t=delta_x/c   #Maximum time step that satisfies the stability condition
#print("delta_t=",delta_t)

Nt = int(t_max/delta_t)
st=delta_x/delta_t 
print("c-st=",c-st) #Checking stability conditions(Unstable if negative)

u =  np.zeros([Nx,Nt])

#boundary condition
u[0,:] = 0
u[-1,:] = 0

#Initial conditions
for i in range (Nx) :
    if i <= 0.8*Nx :
        u[i,0]=1.25*i/Nx
    
    else :  
        u[i,0] = 5.0*(1-i/Nx)

for i in range(1,Nx-1): # u[:,1]settings of
    u[i,1]  = u[i,0]+0.5*((delta_t/delta_x)**2)*(c**2)*(u[i+1,0]+u[i-1,0]-2*u[i,0])

for j in range(Nt-1):
    for i in range(1,Nx-1):
        u[i,j+1] = 2*u[i,j]-u[i,j-1]+(c**2/st**2)* (u[i+1,j]+u[i-1,j]-2*u[i,j])
        

x=list(range(Nx))
y=list(range(Nt))

X, Y = np.meshgrid(x,y)

def functz(u):
    z=u[X,Y]
    return z

Z = functz(u)
fig = plt.figure()
ax = Axes3D(fig)
ax.plot_wireframe(X,Y,Z, color='r')
ax.set_xlabel('x')
ax.set_ylabel('t')
ax.set_zlabel('U')

plt.show()

Result (2) One-dimensional wave equation

Figure 3. Solving the one-dimensional wave equation

Next, the result is illustrated by animation.

%matplotlib nbagg 
from matplotlib.animation import ArtistAnimation #Import methods for creating animations

fig = plt.figure()

anim = [] #A list for storing the data of the para-para diagram drawn for animation

for i in range(Nt):
    U=list(u[:,i])
    x=list(range(Nx))
    if i % int(Nt*0.01) ==0: 
        im=plt.plot(x,U, '-', color='red',markersize=10, linewidth = 2, aa=True)
        anim.append(im)

anim = ArtistAnimation(fig, anim) #Animation creation
plt.xlabel('x')
plt.ylabel('U')

fig.show() 
anim.save("t.gif", writer='imagemagick')   #Animation.Save as gif and create a gif animation file.

Animation of the solution of the one-dimensional wave equation

Addendum

(1) FTCS method

The ** FTCS method (Forward in Time and Centered Difference method) [2] ** in PDEs is a simple scheme that takes the forward difference for the time derivative and the central difference for the spatial derivative. ** Easy to implement, but need to be timed finely to avoid numerical instability **

(2) The simulation dealt with here does not consider the damping effect. Moreover, the time / space dependence of the phase velocity is not considered. References [4] provide a rudimentary explanation of these.

References

[1] Yosuke Nagaoka, ["Vibration and Waves"](https://www.amazon.co.jp/%E6%8C%AF%E5%8B%95%E3%81%A8%E6%B3%A2 -% E9% 95% B7% E5% B2% A1-% E6% B4% 8B% E4% BB% 8B / dp / 4785320451) [2] William H. Press, ["Numerical Recipe in Sea Japanese Version-Recipe for Numerical Calculation in C Language"](https://www.amazon.co.jp/%E3%83%8B% E3% 83% A5% E3% 83% BC% E3% 83% A1% E3% 83% AA% E3% 82% AB% E3% 83% AB% E3% 83% AC% E3% 82% B7% E3% 83% 94-% E3% 82% A4% E3% 83% B3-% E3% 82% B7% E3% 83% BC-% E6% 97% A5% E6% 9C% AC% E8% AA% 9E% E7 % 89% 88-C% E8% A8% 80% E8% AA% 9E% E3% 81% AB% E3% 82% 88% E3% 82% 8B% E6% 95% B0% E5% 80% A4% E8 % A8% 88% E7% AE% 97% E3% 81% AE% E3% 83% AC% E3% 82% B7% E3% 83% 94-Press-William-H / dp / 4874085601 / ref = sr_1_1? S = books & ie = UTF8 & qid = 15039997535 & sr = 1-1 & keywords =% E3% 83% 8B% E3% 83% A5% E3% 83% BC% E3% 83% A1% E3% 83% AA% E3% 82% AB% E3% 83% AB% E3% 83% AC% E3% 82% B7% E3% 83% 94) [3] Qiita article: [Science and technology calculation by Python] Drawing and visualization of 2D (color) contour lines, matplotlib [4] Rubin H. Landau, ["Computational Physics Application" (translation)](https://www.amazon.co.jp/%E8%A8%88%E7%AE%97%E7%89% A9% E7% 90% 86% E5% AD% A6-% E5% BF% 9C% E7% 94% A8% E7% B7% A8-Landau-Rubin-H / dp / 4254130872 / ref = sr_1_2? S = books & ie = UTF8 & qid = 15039997611 & sr = 1-2 & keywords = Landau + computational + physics)