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?