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 codedef Upwind(N,x0,xe,t0,te,m,c,alpha,f,boundary='Dirichlet'):# ----------------------### Timestept = np.linspace(t0, te, N)dt = t[1]-t[0]# ----------------------### Direction conditionsif 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 testUpwind(200,-40,40,0,40,100,-2,1,lambda x: np.exp(-x**2),boundary='Periodic')
Neumann boundary conditions for the upwind scheme applied to the advection equation (Python)
python;boundary conditions;advection