Strange oscillation when solving the advection equation by finite-difference with fully closed Neumann boundary conditions (reflection at boundaries)

Question

I am trying to solving the advection equation but have a strange oscillation appearing in the solution when the wave reflects from the boundaries. If anybody has seen this artefact before I would be interested to know the cause and how to avoid it!

This is an animated gif, open in separate window to view the animation (it will only play once or not at once it has been cached!)

Notice that the propagation seems highly stable until the wave begins to reflect from the first boundary. What do you think could be happening here? I have spend a few days double checking my code and cannot find any errors. It is strange because there seems to be two propagating solutions: one positive and one negative; after the reflection from the first boundary. The solutions seems to be travelling along adjacent mesh points.

The implementation details follow.

The advection equation,

$\frac{\partial u}{\partial t} = \boldsymbol{v}\frac{\partial u}{\partial x}$

where $\boldsymbol{v}$ is the propagation velocity.

The Crank-Nicolson is an unconditionally (pdf link) stable discretization for the advection equation provided $u(x)$ is slowly varying in space (only contains low frequencies components when Fourier transformed).

The discretization I have applied is,

$ \frac{\phi_{j}^{n+1} - \phi_{j}^{n}}{\Delta t} = \boldsymbol{v} \left[ \frac{1-\beta}{2\Delta x} \left( \phi_{j+1}^{n} - \phi_{j-1}^{n} \right) + \frac{\beta}{2\Delta x} \left( \phi_{j+1}^{n+1} - \phi_{j-1}^{n+1} \right) \right]$

Putting the unknowns on the right-hand side enables this to be written in the linear form,

$\beta r\phi_{j-1}^{n+1} + \phi_{j}^{n+1} -\beta r\phi_{j+1}^{n+1} = -(1-\beta)r\phi_{j-1}^{n} + \phi_{j}^{n} + (1-\beta)r\phi_{j+1}^{n}$

where $\beta=0.5$ (to take the time average evenly weighted between the present and future point) and $r=\boldsymbol{v}\frac{\Delta t}{2\Delta x}$.

These set of equation have the matrix form $A\cdot u^{n+1} = M\cdot u^n$, where,

$ \boldsymbol{A} = \left( \begin{matrix} 1 & -\beta r & & & 0 \\ \beta r & 1 & -\beta r & & \\ & \ddots & \ddots & \ddots & \\ & & \beta r & 1 & -\beta r \\ 0 & & & \beta r & 1 \\ \end{matrix} \right) $

$ \boldsymbol{M} = \left( \begin{matrix} 1 & (1 - \beta)r & & & 0 \\ -(1 - \beta)r & 1 & (1 - \beta)r & & \\ & \ddots & \ddots & \ddots & \\ & & -(1 - \beta)r & 1 & (1 - \beta)r \\ 0 & & &-(1 - \beta)r & 1 \\ \end{matrix} \right) $

The vectors $u^n$ and $u^{n+1}$ are the known and unknown of the quantity we want to solve for.

I then apply closed Neumann boundary conditions on the left and right boundaries. By closed boundaries I mean $\frac{\partial u}{\partial x} = 0$ on both interfaces. For closed boundaries it turns out that (I won't show my working here) we just need to solve the above matrix equation. As pointed out by @DavidKetcheson, the above matrix equations actually describe Dirichlet boundary conditions. For Neumann boundary conditions,

$ \boldsymbol{A} = \left( \begin{matrix} 1 & 0 & & & 0 \\ \beta r & 1 & -\beta r & & \\ & \ddots & \ddots & \ddots & \\ & & \beta r & 1 & -\beta r \\ 0 & & & 0 & 1 \\ \end{matrix} \right) $

Update

The behaviour seems fairly independent of the choice of constants I use, but these are the values for the plot you see above:

• $\boldsymbol{v}$=2
• dx=0.2
• dt=0.005
• $\sigma$=2 (Gaussian hwhm)
• $\beta$=0.5

Update II

A simulation with non-zero diffusion coefficient, $D=1$ (see comments below), the oscillation goes away, but the wave no longer reflects!? I don't understand why?

Answer 1

The equation you're solving does not permit right-going solutions, so there is no such thing as a reflecting boundary condition for this equation. If you consider the characteristics, you'll realize that you can only impose a boundary condition at the right boundary. You are trying to impose a homogeneous Dirichlet boundary condition at the left boundary, which is mathematically invalid.

To reiterate: the method of characteristics says that the solution must be constant along any line of the form $x-\nu t = C$ for any constant $C$. Thus the solution along the left boundary is determined by the solution at earlier times inside your problem domain; you cannot impose a solution there.

Unlike the equation, your numerical scheme does admit right-going solutions. The right-going modes are referred to as parasitic modes and involve very high frequencies. Notice that the right-going wave is a sawtooth wave packet, associated with the highest frequencies that can be represented on your grid. That wave is purely a numerical artifact, created by your discretization.

For emphasis: you have not written down the full initial-boundary value problem that you are trying to solve. If you do, it will be clear that it is not a mathematically well-posed problem.

I'm glad you posted this here, though, as it's a beautiful illustration of what can happen when you discretize a problem that's not well-posed, and of the phenomenon of parasitic modes. A big +1 for your question from me.

Answer 2

I learned a lot from the answers above. I want to include this answer because I believe it offers different insights to the problem.

Let us consider the equation \begin{eqnarray} u_{xx} + \frac{1}{c} u_{tt} = 0. \end{eqnarray} with constant wavespeed $c$.

Without initial and boundary conditions this equation has a solution of the form $u(x,t)=f(x-ct)$. (a pulse moving to the right)

If we impose an initial condition $u(x,t_0)=p(x)$, then the solution of the equation in the interval $x \in (-\infty, \infty)$ is $p[x-c(t-t_0)]$. There is no boundary and that is the solution.

Now, let us assume that we define a bounded domain $x \in [a,b]$, after all computers have limited memory. Then we need to specify the values at $a$ and $b$, otherwise computationally we are stuck .

One way to define the boundary conditions is to use Dirichlet on the left boundary and a condition that is consistent with the solution being propagated. That is, we can define (assume initial time $t_0$) \begin{eqnarray} u(a,t_0)=0 \quad , \quad u(b, t_0) = p[b-c(t-t_0)]. \end{eqnarray}

This produces a pulse running to the right until it disappears on the right edge.

push here for animation on Dirichlet on left boundary

I still get some noise which I cannot understand (anybody could help here please?)

The other option is to impose periodic boundary conditions. That is instead of imposing the Dirichlet boundary condition on the left we can impose the wave packet which is consistent with the boundary on the left. That is:

\begin{eqnarray} u(a,t_0)=p[a-c(t-t_0) ] \quad , \quad u(b, t_0) = p[b-c(t-t_0)]. \end{eqnarray}

However $a-c(t-t_0) <a$ for $t>t_0$ and $c>0$, and since we need to put data inside the interval $[a,b]$ we add one "period" of length $b-a$ and then we find $a-c(t-t_0)+b-a=a+b(t-t_0)$, so that the condition on the left would be $u(a,t_0)=p[b-c(t-t_0]$ (same!) and we would have a pulse exiting on the right and entering on the left.

This link shows what I would call periodic boundary conditions.

I made the animations in python and the scheme is a Crank-Nicholson scheme as indicated in the question here.

I am still struggling with the noise pattern after the wave hits the first (right) boundary.