Linear advection eqution and periodic b.c. implementation

I want to implement the forward in time, centered in space scheme for the linear advection \begin{align} u_t+ a u_x=0 \end{align} with periodic boundary conditions and initial datum $$u(x,0)$$. I know that this scheme is unconditionally unstable and the theory about that, but my question is just on the implementation part

I'm doing this in the interval $$x \in [0,1]$$ and $$t \in [0,1]$$.

The scheme is

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

Say I discretize $$[0,1]$$ with $$M$$ points. So $$dx=\frac{1}{M-1}$$. My main problem is the implementation of the boundary conditions. I overlap the nodes $$x_1=0$$ and $$x_M=1$$.

So, at the first iteration, my scheme is (I omit the time index)

\begin{align} u_1=u_1- \frac{a dt}{2 dx} (u_2-u_0) \end{align}

Now, since $$u_0$$ is not known, I would use the fact that $$x_1=x_M$$, and hence the point at the left of $$x_1$$ is $$x_{M-1}$$ and hence $$u_0=u_{M-1}$$.

Now I would solve this up to node $$x_{M-1}$$, so I will have

\begin{align} u_{M-1}=u_{M-1} - \frac{a dt}{2 dx} (u_{1} - u_{M-2}) \end{align} where I used the fact that $$u_M=u_1$$.

Then I will update $$u_{M}=u_{1}$$.

• @Harry49 If my implementation of the periodic boundary conditions is right and if there's anything wrong in what I've written, since indices are making me crazy
– VoB
Oct 7 '19 at 16:17
• Many thanks! @Harry49
– VoB
Oct 7 '19 at 17:52
• The main problem is that I've seen different implementation, and I still don't understand how to treat the fact that $x_1=x_M$ (i.e. overlapping first and last node)
– VoB
Oct 7 '19 at 19:42
• Why do you solve up to the point $x_{M-1}$? It should be $x_M$. And the point $i=M-1$, you should have $u_{M-1}=u_{M-1} - \frac{a dt}{2 dx} (u_{M} - u_{M-2})$. But the point that you interested is $i=M$, in this case: we have $u_{M}=u_{M} - \frac{a dt}{2 dx} (u_{M+1} - u_{M-1})$ but $u_{M+1} = u_{2}$ Oct 7 '19 at 20:02
• @Sesame I solve up to $M-1$ and the impose $u_M=U_1$. You're right for the expression for $u_{M-1}$, it was a typo, now I fixed it
– VoB
Oct 7 '19 at 20:08

We give the following space discretization : $$\forall i\in\{1,\dots,M\}$$: $$\{x_1=0,x_2,\dots,x_M=l\}$$ (I started to count from 1 on purpose as you code in Matlab...). Suppose that we want to solve the above PDE by using an Euler scheme: \begin{align} u_{i}^{n+1}= u_{i}^n-\gamma u_{i+1}^n + \gamma u_{i-1}^n \end{align} where $$\gamma = \frac{a dt}{2dx}$$. We impose periodic boundary conditions, i.e. $$u(0,t)=u(l,t)$$. This means $$\forall n$$, \begin{align} u_0^n&=u^n_M\\ u_{M+1}^n &= u_1^n\\ \end{align} Let's look at know what happen at $$i=1$$ and $$i=M$$: \begin{align} i=1 , \quad u_{1}^{n+1} &= u_{1}^n-\gamma u_{2}^n + \gamma u_{0}^n \\ &= u_{1}^n-\gamma u_{2}^n + \gamma u_{M}^n \\ i=M, \quad u_{M}^{n+1} &= u_{M}^n-\gamma u_{M+1}^n + \gamma u_{M-1}^n \\ &= u_{M}^n-\gamma u_{1}^n + \gamma u_{M-1}^n \\ \end{align} Writting the above in matrix form, you have for $$n\in \{1,\dots, N\}$$ : $$\begin{equation} \mathbf{U}^{n+1} = \mathbf{A}\mathbf{U}^{n} \end{equation}$$ Which can rewritten as for $$M=5$$ \begin{align} \begin{pmatrix} u_1^{n+1}\\ u_2^{n+1}\\ u_3^{n+1}\\ u_4^{n+1}\\ u_5^{n+1}\\ \end{pmatrix}= \begin{pmatrix} &1 &-\gamma &0 &0 &\gamma \\ &\gamma &1 &-\gamma &0 &0 \\ &0 &\gamma &1 &-\gamma &0 \\ &0 &0 &\gamma &1 &-\gamma \\ &-\gamma &0 &0 &\gamma &1 \\ \end{pmatrix}\begin{pmatrix} u_1^n\\ u_2^n\\ u_3^n\\ u_4^n\\ u_5^n\\ \end{pmatrix} \end{align}

• Thanks, this makes sense. Just one more question. With this discretization, is it true that $u_1=u_M$?
– VoB
Oct 7 '19 at 20:54
• Yes. I think in last row I have a coefficient of 1 for $u_1$. So that $u_{M}^{n+1} = u_{1}^n-\gamma u_{2}^n + \gamma u_{M-1}^n$ Oct 7 '19 at 21:01
• Sorry but I do not agree, because that expression does not imply that $u_{M}=u_1$. Moreover, in your matrix the $1$ should be in entry (m,m), as the scheme told us
– VoB
Oct 7 '19 at 21:14
• @VoB: Sorry I was not thinking straight on this one. I made some modifications at my initial post. Hope it is clearer now. Oct 8 '19 at 0:03
• so, the periodicity is imposed by using a value outside of the domain, i.e. $u_0$, not imposing that the first and last value nodes $u_1,u_M$ are equal, right?
– VoB
Oct 8 '19 at 7:23