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.

  • 1
    $\begingroup$ What spatial discretization are you using? $\endgroup$ May 3, 2017 at 3:33
  • 1
    $\begingroup$ First of all, you are using a first-order absorbing conditions. As such, you cannot expect it to have zero reflections. Also, your condition is on the derivatives, and will not absorb the constant components of $\phi$ $\endgroup$
    – DanielRch
    May 3, 2017 at 15:51
  • $\begingroup$ I agree with the comment from @DanielRch; your BC definitions are wrong. Also, this is a two-way wave equation so your Gaussian profile should split and travel both to the left and right. Regarding the correct BC, take a look at this existing post: scicomp.stackexchange.com/questions/26210/… $\endgroup$ May 3, 2017 at 21:22
  • $\begingroup$ Thank you for the comments. I set the wave to be moving to the left by setting the initial $\xi$ and $\Pi$ variables to satisfy a left moving advection equation, $\partial_t\phi-\partial_r\phi=0$, so I use $\xi=\partial_r\phi$ and $\pi = \xi$ for the initial condition. In terms of the boundary condition, @BillGreene I think what you're saying is the correct boundary condition would be, for example $\partial_r\phi + \partial_r\phi=0$ for the boundary on the right. I would like to write down boundary conditions for the variables $\xi$ and $\Pi$ though. $\endgroup$ May 4, 2017 at 1:39
  • $\begingroup$ @Justin, do you notice that your problem, written out in $\xi$ and $\Pi$ will not give a unique solution for $\phi$? An arbitrary constant can be added to $\phi$. There may be a larger problem than the absorbing conditions here. Or are you just solving for the "potentials" $\xi$ and $\Pi$ in your model and ignoring $\phi$ altogether? $\endgroup$
    – DanielRch
    May 4, 2017 at 5:42


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

Browse other questions tagged or ask your own question.