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


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
    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
    c=dt/(2*dx) #multiplication factor

        print('stability issues')
        dt=np.min([tf-t,dt]) #to not exceed final time
        for i in range(1,nx-1): #from x_1 to x_N-2
        #update boundary values

            print('dt resized since excedeed final time ') 
    return np.append(u,u[0]) #add also the value at point x_n, which is equal to value at point x_0



  • $\begingroup$ What do you use for programming? If you use Python, numpy.gradient will impose periodic boundary condition for you automatically when you discretize $u_{x}$. $\endgroup$ Oct 23 '19 at 23:10
  • $\begingroup$ @AloneProgrammer I use python, but my aim is to write my own function, so I want to understand if my argument above is correct. $\endgroup$ Oct 23 '19 at 23:16
  • 1
    $\begingroup$ OK, so your question seems more about mathematics of PBC not how to implement it right? If I got you correctly you want to understand the concept of PBC right? Otherwise, If you use Python, I don't the see the point to not using numpy.gradient unless you want to develop and implement something new that is more efficient/accurate/robust/reliable. $\endgroup$ Oct 23 '19 at 23:18
  • $\begingroup$ The implementation will be the final part. First I want to understand what's behind that conditions, as you said ! $\endgroup$ Oct 23 '19 at 23:38
  • $\begingroup$ Is my answer bad posed? Because I don't understand why it has been downvoted $\endgroup$ Oct 24 '19 at 16:25

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

  • $\begingroup$ Thanks for the answer, I already consulted that link but lots of other answers on the net really confused me! Yes I'm aware of Von Neumann stab. analysis for this method but also for the others in the link, my problem was really to understand how to impose the boundary values! So, if I discretize with $N+1$ points ($x_0, \ldots, x_{N}$), and I actually solve from $x_0$ up to $x_{N−1}$ (included), then I will have the solution for all points except the last one, since it is equal to the value at $x_0$. When I plot the solution, should I also add the last point $x_N=1$ and the value $u_N$? $\endgroup$ Oct 25 '19 at 10:07
  • $\begingroup$ @bobinthebox: You ask whether you should add the lat point when you plot the solution. What difference would it make? The solution up to $x_{N-1}$ is going to look exactly the same with as it does without it. So it makes little difference whether you plot it or not. $\endgroup$ Oct 25 '19 at 17:12
  • $\begingroup$ Alright! I just wanted to plot it to see the solution in the whole domain! I just edited my question with a bunch of code that summarize what you said! I know there are better ways (as roll function) to implement it, but I just need to be sure that the update is done correctly. Numerical experiments seems to confirm, but I'd like to have another confirmation! $\endgroup$ Oct 25 '19 at 17:17
  • $\begingroup$ @bobinthebox: I think you misunderstand what it means for a method to be unconditionally unstable. It makes no sense to write a code that implements such a method. You will never get a solution that is not growing exponentially. Looking at your code, I don't understand why you have the check whether $\Delta t > \Delta x^2$, this check makes no sense from a dimensional point of view. $\endgroup$ Oct 25 '19 at 18:01
  • $\begingroup$ Actually that solutions is not growing exponentially because from a stability analysis we have the restriction $dt <dx^2$. The FTCS scheme in unconditionally unstable in the Von Neumann sense, which is an indicator (in this case) of a really severe restriction $\endgroup$ Oct 25 '19 at 18:32

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.