Imposing periodic boundary condition for linear advection equation - Node problem

Question

I've spent the whole day trying to figure out what is the correct way to impose (and implement) periodic boundary conditions $u(0,t)=u(1,t)$ for all $t>0$ for the simple advection equation $u_t + v u_x=0$, given an initial data $u_0(x)$, but I'm not sure if I'm correct.

I need to understand what am I missing in order to implement it in the right way.

---

Say I discretize $[0,1]$ with $N+1$ points, i.e. I have \begin{align}0=x_0, \ldots, x_N=1 \end{align}

In the lecture we say that for periodic boundary conditions we identify with $x_N$ with $x_0$.

Let's consider the FTCS scheme:

\begin{align} u_{i}^{n+1}=u_i^n-\frac{v dt}{2 dx} (u_{i+1}^{n}-u_{i-1}^{n}) \end{align}

Question: should I do the update for $i=0$ to $i=N$, or should I just do the update for $i=0$ to $i=N-1$ and then update manually $u_N=u_0$?

The problems arise when $i=0$, in this case I have

\begin{align} u_0^n=u_0^n - \frac{v dt}{2 dx} (u_1-u_{-1}) \end{align}

Now I've seen so many things on this site and on the internet that I'm really confused on what is the right value to give to $u_{-1}$.

If $x_N=x_0$, then I would say that the right value for this term is \begin{align} u_{-1}=u_{N-1} \end{align}

since $u(x_{-1})=u(x_0-dx)=u(x_N-dx)=u(x_{N-1})$

But, following the first answer here, first he said that $x_0$ is identified with $x_N$ and then he write that the point to the left of $x_0$ is $x_N$, and the I have \begin{align} u_{-1}=u_N \end{align}

but how is this possible that $x_N$ is at the left of $x_0$, if they coincide? I don't know what I'm missing.

Supposing I do the update also for $i=N$, I obtain (drop the time index $n$ to be more clear)

\begin{align} u_{N}=u_N - \frac{v dt}{2 dx} (u_{N+1}-u_{N-1}) \end{align}

Again, I would say that \begin{align} u_{N+1}=u_1 \end{align} since $u(x_{N+1})=u(x_N+h)=u(x_0+h)=u(x_1)$, but the in the linked answer he has that $u_{N+1}=u_0$.

---

EDIT

Here a runable Python code for the FTCS scheme with periodic boundary conditions and initial value $\sin( 2 \pi x)$, is this the right way to implement it?

import numpy as np
import matplotlib.pyplot as plt
from math import pi

def u0(x):
    return np.sin(2*pi*x)

def FTCS(dx,dt,tf):
    #dx: space step size
    #dt: time step size
    #tf=final time
    nx=np.int(np.ceil(1/dx))
    x=np.linspace(0,1,nx+1)
    x=x[0:-1] #leave out x_N, last point. I solve from x_0 to x_N-1
    u=u0(x) #size N, with values from x_0 to x_N-1
    t=0
    c=dt/(2*dx) #multiplication factor

    if(dt>dx**2):
        print('stability issues')
    while(t<tf):
        dt=np.min([tf-t,dt]) #to not exceed final time
        un=u.copy()        
        for i in range(1,nx-1): #from x_1 to x_N-2
            u[i]=un[i]-c*(un[i+1]-un[i-1])
        #update boundary values
        u[0]=un[0]-c*(un[1]-un[nx-1])
        u[nx-1]=un[nx-1]-c*(un[0]-un[nx-2])

        if(t+dt>tf):
            print('dt resized since excedeed final time ') 
            dt=tf-t 
            u=un.copy()
        t=t+dt
    return np.append(u,u[0]) #add also the value at point x_n, which is equal to value at point x_0

nx=100
dx=1/nx
dt=dx**2
tf=1

x=np.linspace(0,1,nx+1)
plt.plot(x,u,'-',x,u0(x-tf),'-')
plt.show()

Answer 1

Two comments.

First, about implementing periodic boundary conditions. There are several ways to implement such boundary conditions. In my opinion, the easiest is to directly incorporate them into your approximation. Say you have $N+1$ points numbered $0$, $1$, ..., $N$. What periodicity means is that the solution at points $0$ and $N$ are identical, so you don't need to store the solution twice. Furthermore, periodicity also means is that if you need the solution to the left of point $0$, i.e., at the fictitious point $-1$, you should actually look to the left of point $N$, i.e. at the point $N-1$. The same goes if you need the solution to the right of point $N$.

So, your method at point $0$ (which is the same as at point $N$),

$$u_0^{n+1} = u_0^n - \frac{v\Delta t}{2\Delta x}(u_{1}^n-u_{-1}^n)$$

should actually be interpreted as

$$u_0^{n+1} = u_0^n - \frac{v\Delta t}{2\Delta x}(u_{1}^n-u_{N-1}^n)$$

Similar, your method at point $N-1$,

$$u_{N-1}^{n+1} = u_{N-1}^n - \frac{v\Delta t}{2\Delta x}(u_{N}^n-u_{N-2}^n)$$

should actually be interpreted as

$$u_{N-1}^{n+1} = u_{N-1}^n - \frac{v\Delta t}{2\Delta x}(u_{0}^n-u_{N-2}^n)$$

At all the other (interior) points, the method is unchanged. Periodic boundaries can be confusing at first. In my experience, it helps to draw this (which I can't do right now) and think through it step-by-step.

The second and more important point is that the FTCS method is unconditionally unstable. What this means is that if you were to implement it, you would never obtain an approximate solution because it would always diverge (grow without bound) irrespective of how small you make your time step. You can easily establish the unconditional instability by applying the von Neumann method, see Section 2.1 in these notes. To get a simple, stable method, consider the FTBS method, see Section 2.2 in the notes. (By the way, these notes also discuss the implementation of periodic boundary conditions for the FTBS method.)