I am solving the one dimensional wave equation:
$0=\Box\phi = -\partial_t^2\phi + \partial_r^2\phi ,$
using a Crank-Nicolson finite difference scheme, in the domain $r\in[0,R]$. First, I define $\xi\equiv\partial_r\phi$ and $\Pi\equiv\partial_t\phi$, so that the equations of motion can be written as follows:
$ 0 = \partial_t\xi - \partial_r\Pi, \\ 0 = \partial_t\Pi - \partial_r\xi, $
Specifically, the stencils I'm using for the interior grid points are:
$ (\partial_t u)^{n+1/2}_j = \frac{u^{n+1}_j - u^n_j}{\Delta t} + \mathcal{O}\left((\Delta t)^2\right) \\ (\partial_r u)^{n+1/2}_j = \frac{1}{4\Delta r}\left(u^{n+1}_{j+1} - u^{n+1}_{j-1} + u^n_{j+1} - u^n_{j-1}\right) + \mathcal{O}\left((\Delta t)^2, (\Delta r)^2\right) $
The Crank-Nicolson scheme is implicit but I am solving the implicit equations using Gauss-Seidel iteration to solve for $\xi^{n+1}_j$ and $\Pi^{n+1}_j$. I'm having no trouble evolving my initial data (a Gaussian profile for the $\phi$ field, moving the the left) in the interior grid points. I'm only running into trouble at the boundaries of my simulation.
I would like to implement outgoing wave boundary conditions at the boundaries of my simulation, i.e. at $r=0$
$ 0 = \partial_t\xi - \partial_r\xi, \\ 0 = \partial_t\Pi - \partial_r\Pi, $
and at $r=R$
$ 0 = \partial_t\xi + \partial_r\xi \\ 0 = \partial_t\Pi + \partial_r\Pi $
For some reason I am having a hard time numerically implementing these boundary conditions. I have tried several approaches (including first and second order upwind schemes) to implement the boundary conditions but when I run my code either I get what appears to be a numerical instability or the the wave solution 'bounces back' into my domain; i.e. it doesn't exit my domain contrary to the boundary condition!
I've been having difficulty finding a good reference/discussion on outgoing wave boundary conditions. Could someone provide some insight on how to implement outgoing wave boundary conditions using an implicit method? Thank you in advance for the help.