# Open boundary condition for 1d wave equation with variable wave speed using finite differences

I have implemented a finite difference solver for the 1d wave equation with variable wave speed:

$$u_{tt} = c(x)u_{xx}, \hspace{10mm}c(x) = \dfrac{6 -x^2}{2} \hspace{5mm}$$ on $$-2 \leq x \leq 2, t > 0$$ with initial conditions: $$u(x,0) = \cos(\dfrac{\pi x}{0.2}), \hspace{5mm} u_t(x,0) = 0$$ for $$-0.1 \leq x \leq 0.1$$.

I use second order finite differences in order to solve for each grid point in time: $$\frac{u_j^{n+1} - 2u_j^n + u_j^{n-1}}{k^2} = c(x_j)\frac{u_{j+1}^n - 2u_j^n + u_{j-1}^n}{h^2}$$

Rearranging for $$u_j^{n+1}$$: $$u_{j}^{n+1} = 2u_j^n - u_j^{n-1} + c(x_j)q^2(u_{j+1}^n - 2u_j^n + u_{j-1}^n)$$

Where $$q = k/h$$. The boundary condition to be enforced is the open condition, so $$qu_{x} = u_{t}$$ at $$x = -2$$, and $$qu_{x} = -u_{t}$$ at $$x = 2$$. Using the first order finite differences, we can obtain an equation for the ghost points at the boundaries:

$$q(u_{1}^{n} - u_{-1}^n) = u_0^{t+1} - u_0^{t-1}$$ at $$x = -2$$.

Next, we can substitute this into the second order equation ($$c(x_j)$$ disappears since it is 1 at $$\pm 2$$) in order to end up with our iterative formula for the boundary at $$x=-2$$:

$$u_{0}^{n+1} = \dfrac{2u_0^n + (q - 1) u_0^{n-1} + 2r^2(u_{1}^n - u_0^n)}{q + 1}$$

The same can be done for the boundary at $$x = 2$$.

My python implementation is shown below.

def one_dimensional_wave_solver(x_steps, t_steps):
# discrete points in space
x = np.linspace(-2, 2, x_steps + 1)
# distance of step
dx = x[1] - x[0]
# discrete points in time
t = np.linspace(0, 10, t_steps + 1)
# time steps
dt = t[1] - t[0]

# define cur_u and previous u
cur_u = np.zeros(x_steps + 1)
prev_u = np.zeros((x_steps + 1, t_steps + 2))

# Define mesh constant
c = 1
mesh = (c * dt) / dx

# Initial conditions
for i in range(0, x_steps + 1):
if (-0.2 <= x[i] and x[i] <= 0.2):
prev_u[i, 0] = np.cos((np.pi * x[i]) / (0.2))
else:
prev_u[i, 0] = 0

# Iterate through time steps
for t in range(0, t_steps + 1):
# For the first time step
if t == 0:
for i in range(1, x_steps):
cur_u[i] = prev_u[i, t] + (cx(x[i]) * (mesh ** 2) / 2) * (prev_u[i - 1, t] - 2 * prev_u[i, t] + prev_u[i + 1, t])
# Boundary conditions
cur_u[0] = prev_u[0, t] +  mesh * (prev_u[1, t] - prev_u[0, t])
cur_u[-1] = prev_u[-1, t] - mesh * (prev_u[-1, t] - prev_u[-2, t])

else:
# equation for computing u
for i in range(1, x_steps):
cur_u[i] = 2 * prev_u[i, t] - prev_u[i, t - 1] + cx(x[i]) * (mesh ** 2) * (prev_u[i - 1, t] - 2 * prev_u[i, t] + prev_u[i + 1, t])
# Boundary conditions
cur_u[0] = (2 * prev_u[0, t] + (mesh - 1) * prev_u[0, t - 1] + (2 * mesh ** 2) * (prev_u[1, t] - prev_u[0, t])) / (mesh + 1)
cur_u[-1] = (2 * prev_u[-1, t] + (mesh - 1) * prev_u[-1, t - 1] + 2 * (mesh ** 2) * (prev_u[-2, t] - prev_u[-1, t])) / (mesh + 1)

# Swap variables
prev_u[:, t + 1] = cur_u.copy()

# Return final u value
return prev_u[:, -1], prev_u


Using this allows me to get an almost fully transparent boundary. But there is still some slight reflection. Is there an issue with the first order equations for the boundaries?

• That doesn't look like a "slight" reflection. It's clearly visible. Without having examined the question in detail, I vote for error. (When I've done similar things years ago, I had really slight reflections of the order $10^{-6}$ or something like that). Commented Mar 26, 2022 at 14:09
• 1. The initial conditions you give do not match the implemented ones, also you do not state the initial conditions outside the interval $[-0.1, 0.1]$ or $[-0.2, 0.2]$. 2. Why do you have only a plain $c$ in your wave equation, and not $u_{tt} = c{^\color{red} 2}(x) u_{xx}$? 3. In what sense is the B.C. you give the open condition? How about the simple outflow $u_{n+1} = u_n, u_{-1} = u_{0}$? 4. The image you provide looks quite underresolved, did you try smaller time-steps? 5. Related to 4: I am not sure that a CFL numbe rof 1 works here, since the max. of $c$ can be 3. Commented Mar 31, 2022 at 9:43