Neumann boundary conditions for the upwind scheme applied to the advection equation (Python)

I'm trying to solve the linear advection equation

$$u_{t} = cu_{x}, \\ x \in [x_{0}, x_{e}], \quad t \in (0, T], \quad c \in \mathbb{R} \\ u(x,0) = f(x)$$

Note that for $c > 0$, the solution is a left travelling wave so requires a boundary condition on the left, while $c < 0$ is a right travelling wave so requires a boundary condition on the right. Now, implementing Dirichlet ($u(\cdot, t) = 0$) and periodic ($u(x_{0},t) = u(x_{e},t)$) boundary conditions has been quite straight forward. However, I've been trying to implement Neumann conditions $u_{x}(\cdot,t) = \alpha$) and I'm not exactly sure how to or if it is possible?

If $c > 0$, then the discretisation in matrix form is given as

$$\begin{pmatrix} u_{1}^{k+1} \\ u_{2}^{k+1} \\ u_{3}^{k+1} \\ \vdots \\ u_{N}^{k+1} \end{pmatrix} = \begin{pmatrix} 1+\gamma \\ -\gamma & 1+\gamma \\ & -\gamma & 1+\gamma \\ & & \ddots & \ddots \\ & & & -\gamma & 1+\gamma \end{pmatrix} \begin{pmatrix} u_{1}^{k} \\ u_{2}^{k} \\ u_{3}^{k} \\ \vdots \\ u_{N}^{k} \end{pmatrix}$$

If we are to use the ghost point approach, then applying a central difference we get

\begin{align} u_{x}(x_{N},t) &= \frac{u_{N+1} - u_{N-1}}{2 \Delta x} \\ &= \alpha \\ \implies u_{N+1} &= u_{N-1} + 2 \alpha \Delta x, \quad \forall k \end{align}

The problem is, nowhere in our spatial discretisation up to $n = N$ do we have our function indexed at $N+1$, so we can't use that approach. If we use a backward difference instead, we get

\begin{align} u_{x}(x_{N},t) &= \frac{u_{N} - u_{N-1}}{\Delta x} \\ &= \alpha \\ \implies u_{N} &= u_{N-1} + \alpha \Delta x, \quad \forall k \end{align}

and hence our final ODE becomes

$$u_{N}^{k+1} = u_{N}^{k} + \alpha \Delta x$$

but that didn't seem to work either. So my question is, how do we implement Neumann boundary conditions using the upwind method?

Here is my code

def Upwind(N,x0,xe,t0,te,m,c,alpha,f,boundary='Dirichlet'):

# ----------------------
### Timestep

t = np.linspace(t0, te, N)
dt = t[1]-t[0]

# ----------------------
### Direction conditions

if c == 0:

print('u(x,t) = f(x)')

else:

if c > 0:

x = np.linspace(xe, x0, N)

else:

x = np.linspace(x0, xe, N)

dx = x[1]-x[0]
r = c*np.ones(N)*dt/dx

# ----------------------
### Stability conditions

if abs(r[0]) > 1+1e-15:

print('r = {}'.format(r[0]))
print('This solution is unstable, abs(r) must be less \
than or equal to one for stability.')

# ----------------------
### Derivative matrix

DL = spdiags([np.ones(N)], (0), N, N).tocsr() # NOT NECESSARY, MORE FOR
# CONTINUITY WITH OTHER SCHEMES
# THAT HAVE A LHS
DR = spdiags([-r, 1+r], (-1, 0), N, N).tocsr()

# ----------------------
### Boundary conditions

if boundary == 'Dirichlet': # AGAIN, NOT NECESSARY, MORE FOR
# CONTINUITY. HERE u(.,t) = 0.

BL = spdiags([0], (0), N, N)
BR = spdiags([0], (0), N, N)

elif boundary == 'Neumann': # NEEDS TO BE FIXED

BL = spdiags([0], (0), N, N)
BR = spdiags([0], (0), N, N)

else:

BL = spdiags([0], (0), N, N)
BR = spdiags([-r], (N-1), N, N)

# ----------------------
### Initial condition

u = f(x)

# ----------------------
### Loop

for k in range(m+1):

u = spsolve(DL+BL, (DR+BR)*u)

if k % 25 == 0:

plt.plot(x,u)


and to test

Upwind(200,-40,40,0,40,100,-2,1,lambda x: np.exp(-x**2),boundary='Periodic')

• A few years ago I was working on a more complicated version of this problem (advection-diffusion-reaction) and arrived at a general solution for that problem using the finite volume method. My notes might be useful danieljfarrell.github.io/FVM/… . I hope you don't mind a comment about your Python code, it looks like MATLAB. Take a look at PEP8 python.org/dev/peps/pep-0008 and the associated tools. – boyfarrell Jul 27 '17 at 21:37
• The PDE is not well-posed with the Neumann boundary condition, so even if you do manage to implement it, don't expect the numerical solution to approximate anything sensible. – ekkilop Aug 16 '17 at 21:34
• @ekkilop Thanks, I didn't actually know that. Do you have a reference for the well posedness? – mattos Aug 18 '17 at 2:37

The problem is, nowhere in our spatial discretisation up to n=N do we have our function indexed at N+1

You just add an extra element (i.e. at N+1).

Then, every loop set it to the value of the last element times $\alpha\Delta x$. When the forward difference is calculated, it uses a value at $i+1$ to work out $i$. So although there's an extra element, it isn't part of the domain and doesn't appear in the result: it is purely a condition on the boundary.

Alternatively, in your code, you can handle the forward difference for the last element as a special case: instead of looking up the next element (i.e. past the last element), use the last element times $\alpha\Delta x$.

BTW: I think Barba introduces boundary conditions and other issues more simply in her 12 Steps to Navier Stokes (which don't use a matrix solver). Might be worth a look; though it is more work.

The time-derivative scheme you are using is referred to as Forward Euler as is one of the simplest schemes out there, which I will assume for the rest of my answer. Even in the case that you would like to extend this, it is quite straightforward.

Coming to your implementation, we basically need to use $$u_{j}^{N+1} = u_{j}^{N} + c \Delta t u_{x}$$

where you replace $u_{x}$ with the numerical derivative at the interior points and with the boundary condition for the end points.

Given that you've stuck to a matrix notation, which complicates the implementation, the only change you'll need to make is

RHS = (DR + BR)*u
if boundary == 'Neumann':
RHS[0] += c*dt*(alpha-u[0]/dx)
u = spsolve(DL+BL, RHS)


This replaces the left boundary derivative with the boundary condition.

The Wikipedia article for the first-order upwind scheme is clear and would recommend going through it (https://en.wikipedia.org/wiki/Upwind_scheme).

You seem to be using explicit Euler in time, so your scheme would be $$u^{k+1}_j = u^k_j + \frac{c \Delta t}{\Delta x} (u^k_j-u^k_{j-1})$$ (if so, I don't really understand why you would make a sparse matrix). With the second approach that would give $$u^{k+1}_N = u^k_N + \alpha c \Delta t$$