# The staggered finite-difference method (交错网格有限差分方法)

## One-dimensional problem

We consider 1d model problem for acoustic problem:$\frac{\partial u}{\partial t} – \frac{\partial v}{\partial x} =0, \quad \frac{\partial v}{\partial t} – \frac{\partial u}{\partial x} =0. \tag{1}$

### FD scheme

Let $$u_i^n$$ be the approximation of $$u(x_i, t_n)$$ for the uniform mesh with $$x_i=i h, t^n = n \tau$$. The finite difference scheme can be expressed as follows\begin{align} \frac{u_i^{n+\frac{1}{2}}-u_i^{n-\frac{1}{2}}}{\tau} – \frac{v_{i+\frac{1}{2}}^n – v_{i-\frac{1}{2}}^n}{h} =0 , \\ \frac{v_{i+\frac{1}{2}}^{n+1}-v_{i+\frac{1}{2}}^{n}}{\tau} – \frac{u_{i+1}^{n+\frac{1}{2}} – u_{i}^{n+\frac{1}{2}}}{h} =0 . \end{align} \tag{2} For other schemes, see, e.g., [1].

We can interpret (2) in another viewpoint. Making a integration of $$u_t-v_x=0$$ on the interval $$(x_{i-1/2}, x_{i+1/2})$$, we have$\int_{x_{i-1/2}}^{x_{i+1/2}} \frac{\partial u}{\partial t}(x,t_n) dx = \int_{x_{i-1/2}}^{x_{i+1/2}} \frac{\partial v}{\partial x}(x,t_n) dx =v(x_{i+1/2},t_n)-v(x_{i-1/2},t_n).$By using the mid-point quadrature formula, we have (note $$x_{i+1/2}-x_{i-1/2}=h$$)$h \frac{\partial u}{\partial t}(x_i,t_n) = v(x_{i+1/2},t_n)-v(x_{i-1/2},t_n).$Next, using the center difference formula for $$\frac{\partial u}{\partial t}$$ at $$t_n$$ leads to (note $$t_{n+1/2}-t_{n-1/2}=\tau$$)$h \frac{u(x_i,t_{n+1/2})-u(x_i,t_{n-1/2})}{\tau} = v(x_{i+1/2})-v(x_{i-1/2}), \tag{3}$which is the same as (2).

### Numerical results

Example 1 (Homogenous BC, [1]). Let $$\Omega=[-L,L]$$. Suppose the acoustics has the pressure$u(x,0)=\begin{cases} 1-\cos (x+\pi), & \text{if } |x|\leq \pi, \\ 0, & \text{if } \pi <|x|\leq L, \\ \end{cases}$The boundary condition is set to be zero, i.e., $$u(-L,t)=u(L,t)=0$$. The velocity satisfies $$v(x,0)=0$$ for all $$x\in \Omega$$. In numerical test, we set $$L=10$$.

Example 2 (Non-homogenous BC). Let $$\Omega=(0,b)$$. Suppose the acoustics has the boundary condition:$\begin{cases} u(0,t)=\sin t, \\ u(b,t)=0, \end{cases}\quad \begin{cases} v(0,t)=0, \\ v(b,t)=0. \end{cases} .$The initial conditions are set to be zero, i.e., $$u(x,0)=v(x,0)=0$$. The velocity satisfies $$v(x,0)=0$$ for all $$x\in \Omega$$.

## Two-dimensional problem

Given a bounded rectangular domain $$\Omega\subset \mathbb{R}^2$$. We consider 2d model problem for acoustic problem: Find $$u\in \mathbb{R},\mathbf{v}=[v,w]^T\in \mathbb{R}^2$$ s.t.$\frac{\partial u}{\partial t} – \nabla\cdot \mathbf{v} =0, \quad \frac{\partial \mathbf{v}}{\partial t} – \nabla u = 0,\quad \text{in }\Omega.$

### Staggered FD scheme 2D

Let $$u_{i,j}^n$$ be the approximation of $$u(x_i, y_i, t_n)$$ for the uniform mesh with $$x_i=i h, y_j=jh, t^n = n \tau$$. The staggered finite difference scheme can be expressed as follows$\frac{u_{i,j}^{n+\frac{1}{2}}-u_{i,j}^{n-\frac{1}{2}}}{\tau} – \frac{v_{i+\frac{1}{2},j}^n – v_{i-\frac{1}{2},j}^n}{h} – \frac{w_{i,j+\frac{1}{2}}^n – w_{i,j-\frac{1}{2}}^n}{h} = 0, \\ \frac{v_{i+\frac{1}{2},j}^{n+1}-v_{i+\frac{1}{2},j}^{n}}{\tau} – \frac{u_{i+1,j}^{n+\frac{1}{2}} – u_{i,j}^{n+\frac{1}{2}}}{h} =0,\\ \frac{w_{i,j+\frac{1}{2}}^{n+1}-w_{i,j+\frac{1}{2}}^{n}}{\tau} – \frac{u_{i,,j+1}^{n+\frac{1}{2}} – u_{i,j}^{n+\frac{1}{2}}}{h} =0.$For the first step $$n=0$$, we change the update procedure for $$u$$ to be$\frac{u_{i,j}^{n+\frac{1}{2}}-u_{i,j}^{n}}{\tau/2} – \frac{v_{i+\frac{1}{2},j}^n – v_{i-\frac{1}{2},j}^n}{h} – \frac{w_{i,j+\frac{1}{2}}^n – w_{i,j-\frac{1}{2}}^n}{h} = 0.$

After allocating memory for $$u_{(N_x+1)\times (N_y+1)}, v_{N_x\times (N_y-1)}, v_{(N_x-1)\times N_y}$$ and updating $$u^{1/2}$$, the main step is to follow the FD scheme. The kernel code is as the following (don’t forget to impose boundary condition after each updating):

for i in range(0,Nx):
for j in range(1,Ny):
v[i,j-1] += (u[i+1,j]-u[i,j])*dt/dx

for i in range(1,Nx):
for j in range(0,Ny):
w[i-1,j] += (u[i,j+1]-u[i,j])*dt/dx

for i in range(1,Nx):
for j in range(1,Ny):
u[i,j] += (v[i,j-1]-v[i-1,j-1])*dt/dx + (w[i-1,j]-w[i-1,j-1])*dt/dy


Example 3 (homogenous BC). Let $$\Omega=[-4,4]\times[-4,4]$$. Suppose the acoustics has the initial pressure$u(x,y,0)=\begin{cases} 1+\cos (x^2+y^2), & \text{if } x^2+y^2\leq \pi, \\ 0, & \text{if } \pi < x^2+y^2 \leq 4. \end{cases}$The boundary condition is set to be zero, i.e., $$u(x,y,t)=0,\ x\in \partial \Omega$$. The initial velocity satisfies $$v(x,0)=0$$ for all $$x\in \Omega$$.

## References

[1] Maria Yuliani Danggo and Sudi Mungkasi. A staggered grid finite difference method for solving the elastic wave equations. 2017. J. Phys.: Conf. Ser. 909 012047.

## Appendix: Python code for Example 1 and 2.

# Staggered finite difference method for acoustic wave.
# u_t - v_x = 0,   v_t - u_x = 0, in [-L,L]
# u(-L,t) = u(L,t) = 0.
#
# class list:
#   Wave1d - pde class
#   Mesh - for mesh
#   Wave1DSolver - Solver
#
# -- https://numanal.com 2022/2/2
#
import numpy as np
import matplotlib.pyplot as plt

class Wave1d(object):
# pde class
def __init__(self, a, b, T):
self.a = a
self.b = b
self.umax = 4
self.T = T

def u(self,x,t):
pass

def u_l(self,t):
# return 0.*t # ex1
return np.sin(np.pi*.5*t) # ex2
def u_r(self,t):
return 0.*t

def u0(self,x):
# ex1
# r = 0*x
# for i in range(r.shape[0]):
#     if np.abs(x[i]) <= np.pi:
#         r[i] = 1-np.cos(x[i]+np.pi)
# return r
# ex2
return 0*x

def v0(self,x):
return 0.*x

def f(self,x,t):
return 0.*x

# uniform spatial and time space
class Mesh():
def __init__(self,x,t):
self.x = x
self.t = t

# solver class
class Wave1DSolver():
def __init__(self, pde, mesh):
self.pde = pde
self.mesh = mesh
self.Nx = mesh.x.shape[0]-1
self.Nt = mesh.t.shape[0]-1
self.dx = mesh.x[1]-mesh.x[0]
self.dt = mesh.t[1]-mesh.t[0]

def solve(self):
dt, dx = self.dt, self.dx
Nx = self.Nx

# set initial values
x = self.mesh.x
u = self.pde.u0(x)             # current u
v = self.pde.v0(x[:-1]+dx*.5)  # current v (staggered mesh)

# compute u^{1/2}
for i in range(1,Nx):
u[i] = u[i] + (v[i]-v[i-1])*.5*dt/dx
# impose boundary condition
u[0] = self.pde.u_l(dt*.5)
u[-1] = self.pde.u_r(dt*.5)

# gif plot - 1/3
plt.figure(figsize=(8, 6), dpi=80)
plt.ion()

for j in range(1,self.Nt):
for i in range(0,Nx):
v[i] += (u[i+1]-u[i])*dt/dx
for i in range(1,Nx):
u[i] += (v[i]-v[i-1])*dt/dx

# u(a,t) and u(b,T)
u[0] = self.pde.u_l(dt*(j+.5))
u[-1] = self.pde.u_r(dt*(j+.5))

# gif plot - 2/3
if j%np.ceil(t.shape[0]/50) == 1:
plt.cla()
plt.title("Wave 1D")
plt.grid(True)
plt.plot(x,u, '-r', label="FDM u: t="+str(int(1e5*self.mesh.t[j+1])/1e5))
plt.plot(x[:-1]+.5*dx,v, '--b', label="FDM v: t="+str(int(1e5*self.mesh.t[j+1])/1e5))
plt.xlim(x[0], x[-1])
plt.ylim(-1.2*self.pde.umax, 1.2*self.pde.umax)
plt.pause(.1)

# gif plot - 3/3
plt.ioff()
plt.show()

return 0

if __name__ == "__main__":
pde = Wave1d(a=0, b=6, T=10.)  # pde
Nx = 200        # Nx
x = np.linspace(pde.a, pde.b, num=Nx, endpoint=True)
dt = 0.05*(pde.b-pde.a)/Nx; print(dt)
Nt = int(pde.T/dt+1); print(Nt)
t = np.linspace(0, pde.T, num=Nt, endpoint=True)
mesh = Mesh(x, t)

solver = Wave1DSolver(pde, mesh)
solver.solve()


0 0 投票数

2 评论

4 月 前