How many boundary conditions do I need to use this finite difference spatial approximation to solve the wave equation?

Question

In the book, on page $233$, I found that the wave equation:

$$ u_{tt}=u_{xx} \tag 1$$

can be discretized using the following finite difference approximation:

$${ U_{i+1,j}-2U_{i,j}+U_{i-1,j} \over \Delta t^2 } = { U_{i+1,j+1} - 2U_{i+1,j} + U_{i+1,j-1} + U_{i-1,j+1} - 2U_{i-1,j} + U_{i-1,j-1}\over 2\Delta x^2} \tag 2$$

where $i$ represents the time step index and $j$ represents the spatial step index. This is the averaged spatial approximation method and it is unconditionally stable apparently. However, in the book, there is no example that shows how to use this approximation. When I reformulate equation $(2)$ I get the following expression:

$$ -kU_{i+1,j+1} + (1 + 2k)U_{i+1,j} -kU_{i+1,j-1} = 2U_{i,j}+kU_{i-1,j+1}-(1 + 2k)U_{i-1,j} + kU_{i-1,j-1} \tag 3$$

$$k = {\Delta t^2 \over 2 \Delta x^2} \tag 4$$

I know the boundary condition at $x=0$ and $x=x_{max}$. However, it seems to me that I need one more boundary condition in order to use this approximation. To calculate $U_{i+1,2}$, I need to know $U_{i+1,1}$ and $U_{i+1,0}$. The problem is, I only know $U_{i+1,0}$. My question is, do I truly need one more boundary condition to use this approximation? If not, how can I use this approximation correctly?

Answer 1

Wave equation PDE is: $u_{tt}=u_{xx}$

Example problem:

Initial condition: $u(0,x)=\sin(\pi x)$

Boundary conditions: $u(t,0)=0, u(t,1)=0$

Neumann condition: $u_{t}(0,x)=0$

General solution is: $u\! \left(t,x\right)=\sin\! \left(\pi x\right)\cdot \cos\! \left(\pi t\right)$

Solution by finite difference method with Python (using finite difference scheme (2)):

import numpy as np
import matplotlib.pyplot as plt

# Define helper functions
def index(p_row, p_col):
    return p_row * nx + p_col

def f(x):
    return np.sin(np.pi * x)
    # return 2 * x if x < 0.5 else 2 - 2 * x

# General solution for reference
def u_ref(t, x):
    return np.cos(np.pi * t) * np.sin(np.pi * x)

# Setup
t_min = 0  # Left endpoint of t interval
t_max = 5  # Right endpoint of t interval
x_min = 0  # Left endpoint of x interval
x_max = 1  # Right endpoint of x interval

nt = 81  # Number of grid points in t
nx = 61  # Number of grid points in x
ht = (t_max - t_min) / (nt - 1)  # Spacing in t direction
hx = (x_max - x_min) / (nx - 1)  # Spacing in x direction
n = nx * nt  # Dimension of system

fac_t = 1 / ht ** 2
fac_x = 1 / (2 * hx ** 2)

A = np.zeros((n, n), dtype='float')
b = np.zeros(n, dtype='float')

# Index of equation
# HINT: i represents the time step, j represents the spatial step index.
row = 0

# Set finite difference approximation
for i in range(1, nt - 1):
    for j in range(1, nx - 1):
        # LHS u_tt
        A[row, index(i + 1, j)] += 1 * fac_t
        A[row, index(i, j)] += -2 * fac_t
        A[row, index(i - 1, j)] += 1 * fac_t

        # RHS u_xx
        A[row, index(i + 1, j + 1)] += -1 * fac_x
        A[row, index(i + 1, j)] += 2 * fac_x
        A[row, index(i + 1, j - 1)] += -1 * fac_x

        A[row, index(i - 1, j + 1)] += -1 * fac_x
        A[row, index(i - 1, j)] += 2 * fac_x
        A[row, index(i - 1, j - 1)] += -1 * fac_x

        b[row] = 0
        row += 1

# Dirichlet boundary conditions
for j in range(nx):
    # u(0,x) = f(x)
    A[row, index(0, j)] = 1
    b[row] = f(x_min + j * hx)
    row += 1

for i in range(1, nt):
    # u(t,0) = 0
    A[row, index(i, 0)] = 1
    b[row] = 0
    row += 1

    # u(t,1) = 0
    A[row, index(i, nx - 1)] = 1
    b[row] = 0
    row += 1

# Neumann boundary condition: u_t(0,x) = 0
for j in range(1, nx - 1):
    if row >= n:
        break
    A[row, index(1, j)] = +1 / ht
    A[row, index(0, j)] = -1 / ht
    b[row] = 0
    row += 1

plt.spy(A)
plt.show()

# Check rank of matrix A
rank = np.linalg.matrix_rank(A)
print('rank=', rank, 'row=', row)

if rank < row:
    raise Exception('Insufficient rank {}: Expecting rank {}'.format(rank, n))

# Solve the linear system
u = np.linalg.solve(A, b)
np.savetxt('uuu.csv', u, delimiter=',')

# Evaluate quality of solution
error = np.matmul(A, u) - b
norm = np.linalg.norm(error)
print('error norm=', norm)

if norm > 1E-9:    
    raise Exception('Error norm: {}'.format(norm))

# Prepare visualization
t = np.linspace(t_min, t_max, nt)  # Grid points in t-direction
x = np.linspace(x_min, x_max, nx)  # Grid points in x-direction
tg, xg = np.meshgrid(t, x)
ug = u.reshape(nt, nx).transpose()
uu = u_ref(tg, xg)

# Plot the solution
fig = plt.figure(figsize=(8, 8))
ax = plt.axes(projection='3d')
surf = ax.plot_surface(tg, xg, ug - uu)
ax.set_xlabel('t')
ax.set_ylabel('x')
ax.set_zlabel('u')
ax.set_title('Deviation from analytical solution nt={}, nx={}'.format(nt, nx))
# ax.azim = 45
plt.show()

Unfortunately the deviation increases over time with this finite difference scheme.