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)$ which is a a gaussian. 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_{M-2} - u_{1}) \end{align} where I used the fact that $u_M=u_1$.

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

The following MatLab is an implementatin of what I've written:

clear all close all M=100; dx=1/(M-1); x=linspace(0,1,M); N=500; dt=1/(N-1); u0=@(x) exp(-50*(x-1/2).^2); T=1; #final time t=0; U=u0(x); v=1; #velocity C=(v*dt)/dx; while t<T U(1:M-1)=U(1:M-1) - (v*dt)/(2*dx)*( U([2:M-1,1]) - U([M-1,1:M-2]) ); U(M)=U(1); t+=dt; plot(x,U,'*',x,u0(x-v*t),'r'); legend('numerical solution','exact solution') title(sprintf('Soluzions a t = %0.2f',t)); pause(0.001) end