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