Instability at the boundary of a finite difference simulation of a hyperbolic PDE

Question

I want to simulate the hyperbolic partial differential equation $$W_{tt} + V W_{tx} + k_E V W_x + k W_t = 0,$$ but I am having trouble finding a discrete analog of this equation which is numerically stable.

I introduced central differences in both space and time. Time and space are discretized via $t_n = n \Delta t$ and $x_i = i \Delta x$ for $j=1,\dots, N_{iter}$ and $i = 1,\dots,N$. The discrete version of $W$ is written $$W(x_i,t_n) = W_i^n.$$ Defining the numerical derivatives as $$\partial_t W_i^{n+1} = \frac{W_i^{n+1}-W_i^{n-1}}{2\Delta t},$$ $$\partial_t^2 W_i^n = \frac{W_i^{n+1}-2 W_i^{n-1}+W_i^{n-1}}{\Delta t^2}$$ and similarly for the space component gives the difference equations $$ -a W_{i-1}^{n+1} + (1+c)W_i^{n+1} + a W_{i+1}^{n+1} = b W_{i-1}^n + 2 W_i^n -b W_{i+1}^n \\ - a W_{i-1}^{n-1} + (c-1) W_i^{n-1} + a W_{i+1}^n.$$ Here the parameters $a,b,c$ involve the coefficients in the differential equation, the timestep $\Delta t$, and the grid spacing $\Delta x$.

Assuming the absorbing boundary condition $W_1^n = W_N^n = 0$, these can be written in a matrix form as $$ A W^{n+1}= B W^n + C W^{n-1},$$ where each of the matrices $A,B,C$ are tri-diagonal and involve constant coefficients composed of $a,b,c$. The vectors are $W^n = [W_2^n, W_3^n, \dots W_{N-1}^n]^T$.

I formulated the matrices, then set up the following type of loop to simulate the PDE:

W0=W1=P # The first two iterations are both to the same initial condition
Ainv = ... # the inverse of A is calculated
for j in range(Niter): # time loop
    W = Ainv@(B@W1+C@W0) # solution of the matrix equation for W
    W0 = W1 # redefine old W
    W1 = W # redefine new W 

The result seems in a sense half correct. Every second point does what I expect... while the remaining points seem to diverge until they overflow. Here is an image of what I mean. This image demonstrates 5 seconds of evolution from an initially triangular profile. On the right one can see the splitting between two different evolutions, one which makes sense (the lower one), and another which doesn't (the upper one).

Am I missing anything here regarding finite difference simulation of hyperbolic equations? Is there a reason central differences would not work here? I suppose the issue is really that the right-most boundary is not properly behaving as an absorbing boundary. Is this a common issue? I am relatively certain that my code is bug-free (and I am happy to share). Thanks for any help!

Answer 1

I've got a two comments that hope can be of help.

1. In my experience, trying to impose boundary conditions analytically and come up with the corresponding matrices for the interior nodes as you have done, is more bug-prone than imposing boundary conditions numerically. Given that you have a problem with boundary conditions, I wouldn't discard easily that your code is indeed bug-free. At first sight, and judging from what you've written, it looks like it's not.

2. A second thing I can add is that, usually, when solving hyperbolic problems involving second time-derivatives, we have initial data on $W$ and $W_t$. You can use this data to make a first iteration (for example with an explicit Euler method, or a better one if you need a higher-order approximation), to come up with an estimate for W1. But you have used

W0=W1=P

This looks odd to me.