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

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.

$\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?

• what did you take for $\mathbf{v}$? – chris Mar 4 '13 at 14:22
• $\boldsymbol{v=2}$ in those simulations. I will update with the simulation setting. Good idea. – boyfarrell Mar 4 '13 at 14:23
• Then I would expect the initial condition to be advected to the right and disappear through the right boundary. All that comes to mind is that the central scheme can give spurious oscillations unless applied to the advection-diffusion equation with cell Peclet number below 2. Perhaps try the upwind scheme? – chris Mar 4 '13 at 14:36
• Do you think there could be sign error with the equation. Actually, my end goal is to apply this with the advection-diffusion equation. I am currently testing different limiting cases. In the above example the diffusion coefficient was set to zero. I have included a new animation above. I don't understand why the peak doesn't reflect when the diffusion coefficient is nonzero? It does exactly as you mentioned (apart from the direction). – boyfarrell Mar 4 '13 at 14:49
• I was thinking of $\partial_t u + v \partial_x u =0$, so the sign is ok. The second plot look ok to me. Why would you expect anything to reflect? That could only happen if $v$ changes sign somehow. Try with the upwind scheme for advection instead of the central scheme, then you should see something similar for $D=0$. – chris Mar 4 '13 at 15:19

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.

• thank you for the discussion and corrections. I hadn't appreciated that the matrix will not have the same properties as the differential equation. More comments to follow... – boyfarrell Mar 4 '13 at 23:59
• Yes, I see now how these are actually Dirichlet boundary conditions I have made a note above for the correction. This is the first time I have tried (in ernest) to really understand the process of solving these equations, I keep on hitting bumps on the road. I'm happy to post my progress! – boyfarrell Mar 5 '13 at 0:42
• @David Ketcheson : I am running into the same problem and posted my problem in the following link scicomp.stackexchange.com/questions/30329/… Can you please explain me why do you say that the problem is not "mathematically well-posed" ? Thanks. – Herman Jaramillo Oct 15 '18 at 13:39
• @HermanJaramillo You are trying to impose a solution value at the left boundary, where it is already determined by the PDE. Any textbook that includes a discussion of the advection will also indicate what the valid boundary conditions are and why. I added a second paragraph with additional explanation; hope that helps. – David Ketcheson Oct 15 '18 at 15:12
• @HermanJaramillo: not "mathematically well-posed" basically means you have two equations for one function value at the boundary, the boundary condition as well as the PDE itself. In general, those two equations contradict each other. More generally, one can consider this as an optimization problem in which both objectives are to be satisfied as good as possible. – davidhigh Oct 17 '18 at 9:09

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) 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.

• I could not see the animation in my mobile phone, but I believe your noise pattern is caused by missing numerical accuracy. The absorption only works because you impose the exact solution on the boundary. W.r.t. this exact solution, the numerical solution arriving at the right boundary slightly differs in frequency and phase. This again leads to reflection and thus interference. – davidhigh Oct 17 '18 at 8:38
• @davidhigh : Thanks for your information. I will check this. I am sorry the animation did not work in your phone. It did not work in my phone (Samsung) either. This might be some missing software in the phones. It should work in a computer. Thanks again. – Herman Jaramillo Oct 17 '18 at 20:53
• You're welcome. The animation itself is not so important, I just wanted to see how large the errors are. Btw, by imposing the exact solution on the boundary, you avoid the "ill-posedness" by design. That is, you still have two equations for one value at the boundary, but by using the exact result you force them to be consistent. But this works only approximately, as the numerical result is not completely exact. – davidhigh Oct 18 '18 at 8:56
• And another comment: the general analytical solution to the wave equation is a superposition of a pulse moving left and one moving right. In your case, you only consider the right-moving pulse, so you already applied initial conditions -- in contrast to what you state in the text. – davidhigh Oct 18 '18 at 9:20
• @davidhigh : I have been thinking a bit on your insights about the noise after the pulse hits the boundary. I believe you are right and there is a difference between the exact analytic solution and the numerical propagated pulse. In the boundary that little difference generates the small noise pattern seen there. The C.N. advection system shown in this discussion is dispersive and I believe that while dispersion is not noticed before the pulse hits the boundary, it might triggers but the small perturbation (difference between analytical and numerical solutions) at the boundary. Thanks again. – Herman Jaramillo Oct 20 '18 at 23:58