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 '17 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 '17 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 '17 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 '17 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 '17 at 5:42

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Browse other questions tagged or ask your own question.