Advection equation with discontinuous initial condition

I'm trying to solve in MatLab/Octave the advection equation, with $a > 0$ :

$$\begin{cases} u_t + au_x=0\\ u(x,0)=u^{0}(x) \end{cases}$$

with $u^{0}(x)=\begin{cases} 1.5 \quad x<0 \\0.5 \quad x>0 \end{cases}$

I have to use Lax-Wendroff's method, and then, after visualizing that the numerical solution is not 'good', use the conservative Lax-Friedrichs' method.

For a equally spaced grid both in $x$ and $t$ -axis, where the step sizes are respectively $h$ and $k$, we set $u(x_j,t_n) \approx U_j^{N}$ and use finite differences schemes in order to get, for the advection equation:

- Lax-Wendroff

$$U_j^{N+1}=U_j^N - \frac{k}{2h}a(U_{j+1}^{N}-U_{j-1}^{N}) + \frac{k^2}{2h^2}a^2(U_{j+1}^{N} - 2U_{j}^{N} + U_{j-1}^{N})$$

- Conservative Lax-Friedrichs' method:

$$U_j^{N+1}=\frac{1}{2}(U_{j+1}^{N}+U_{j-1}^{N}) - \frac{k}{2h}a(U_{j+1}^{N} - U_{j-1}^{N})$$

My doubts are basically on the computation of $U_{0}^{N}$ and $u(x_\text{final},t_n)$ with $N=1,2,\ldots$ and

I thought to use forward and backward finite differences, using the following matlab code and getting the following figure. Please let me know if I am wrong in here

     %% Lax-Wendroff. a is a given positive number
% the coefficients  h,k are taken in order to satisfy CFL condition
% mt=number of time nodes in the t-axis
% mx= number of nodes in the x-axis
% u0 is defined as in the text.

U=zeros(mt,mx);
U(1,:)=u0;
K=(k*a)/(2*h);
KK=((k*a)^2)/(2*h^2); %re-define cfl coeff used in the iteration

for i=1:mt-1
for j=2:mx-1
U(i+1,j)=U(i,j) - K*(U(i,j+1)-U(i,j-1)) + KK*(U(i,j+1)-2*U(i,j)+U(i,j-1));
end
j=mx;
U(i+1,1)=U(i,2) - K*(U(i,2)-U(i,1)); %first node for every time
U(i+1,j)=(1-K^2)*U(i,j)+0.5*K*((K-1)*(2*U(i,j)-U(i,j-1))+(K+1)*U(i,j-1)); %last

plot(x,U(i,:),'b--')
pause(0.01)
end


end

From this image it seems that the numerical solution is not accurate enough, as expected.

\ So I use conservative Lax-Friedrichs' method. Here the problem is, as before, the computation of the values at the first and last node: I used finite differences (again) but I'm not so sure:

This is coded as follows

U=zeros(mt,mx);
U(1,:)=u0;
CFL=k/(2*h);

for i=1:mt-1
for j=2:mx-1
U(i+1,j)= 0.5*(U(i,j+1) + U(i,j-1)) - CFL*(U(i,j+1) - U(i,j-1));
end
j=mx;
U(i+1,1)=U(i,1)-0.5*(k/h)*(U(i,2)-U(i,1)); %first x-axes node
U(i+1,j)=U(i,j)-0.5*(k/h)*(U(i,j)-U(i,j-1)); %last x-axis node
plot(x,U(i,:),'b--',x,u0,'r')
legend('numerical solution','initial condition u(x,0)')
pause(0.01)
end


Here I got the following graphic:

Now, as you can see, the numerical solution is a little better from the previous, but I think it should not be so 'smooth' in the discontinuity.

I guess this is due to my approach in computing the values in the first and last node. How should I approximate those nodes?

EDIT In order to see how the initial condition is trasported, using the code of @Harry49, I should write:

while (t+k)<tf
for j=2:Mx
% Lax-Wendroff
Utemp(1,j) = U(1,j) - 0.5*k/h*a*(U(1,j+1)-U(1,j-1)) + 0.5*(k/h*a)^2*(U(1,j+1)-2*U(1,j)+U(1,j-1));
% Lax-Friedrichs
Utemp(2,j) = 0.5*(U(2,j-1)+U(2,j+1)) - 0.5*k/h*a*(U(2,j+1)-U(2,j-1));
end
% Ghost cells
Utemp(:,1)    = Utemp(:,2);
Utemp(:,Mx+1) = Utemp(:,Mx);
% Update
U = Utemp;
t = t + k;
plot(x,U(1,:),'bo',x,U(2,:),'r--',[-2,a*t,a*t,2],[u0(-1),u0(-1),u0(1),u0(1)],'k-')
xlabel('x')
ylabel('u(t(i),x)')
legend('Lax Wendroff','Lax Friedrichs','Initial conditon')
axis([-2,3,0,2])
pause(0.01)

end


Indeed, both schemes have the same stencil. The values $U_{j-1}^{N}$ and $U_{j+1}^{N}$ are required to compute $U_{j}^{N+1}$. The classical way of handling the boundaries of the numerical domain is to consider $U_{0}^{N}$ and $U_{M}^{N}$ as ghost cells [1], which are not updated using the scheme. They are updated according to different rules, which may account for various boundary conditions. Here, to represent an infinite domain, the simplest is to use a piecewise constant extrapolation by setting \begin{aligned} U_{0}^{N} &= U_{1}^{N} \\ U_{M}^{N} &= U_{M-1}^{N} \end{aligned} at each time step.

Here is a runable Matlab code:

% parameters
Mx = 200;
CFL = 0.95; % Courant number
a = 1;
tend = 1;

u0 = @(x) 1.5*(x<0) + 0.5*(x>=0);
x = linspace(-2,2,Mx+1);
h = 4/Mx;
k = CFL*h/a;

t = 0;
U = zeros(2,Mx+1);
U = [u0(x); u0(x)];
Utemp = U;

figure;
plot([-2,0,0,2],[u0(-1),u0(-1),u0(1),u0(1)],'k--');
hold on
hth = plot([-2,a*t,a*t,2],[u0(-1),u0(-1),u0(1),u0(1)],'k-');
hLW = plot(x,U(1,:),'bo');
hLF = plot(x,U(2,:),'r.');
hold off
xlim([-2 2]);
ylim([0.2 1.8]);
legend('init','theo','LW','LF');
ht = title(sprintf('t = %0.3f',t));
xlabel('x');
ylabel('u');

while (t+k)<tend
for j=2:Mx
% Lax-Wendroff
Utemp(1,j) = U(1,j) - 0.5*k/h*a*(U(1,j+1)-U(1,j-1)) + 0.5*(k/h*a)^2*(U(1,j+1)-2*U(1,j)+U(1,j-1));
% Lax-Friedrichs
Utemp(2,j) = 0.5*(U(2,j-1)+U(2,j+1)) - 0.5*k/h*a*(U(2,j+1)-U(2,j-1));
end
% Ghost cells
Utemp(:,1)    = Utemp(:,2);
Utemp(:,Mx+1) = Utemp(:,Mx);
% Update
U = Utemp;
t = t + k;
set(hth,'XData',[-2,a*t,a*t,2]);
set(hLW,'YData',U(1,:));
set(hLF,'YData',U(2,:));
set(ht,'String',sprintf('t = %0.3f',t));
drawnow;
end


And the corresponding graphics output:

[1] R.J. LeVeque, Finite Volume Methods for Hyperbolic Problems, Cambridge University Press, 2002, doi:10.1017/CBO9780511791253

• Yes, that was the problem: according to your MatLab code, if I would like to see how the initial condition is trasported along the x-axis, I could write the code I posted in my edit, isn't it? – VoB Apr 28 '18 at 15:20
• yes, that was the idea ! – VoB Apr 28 '18 at 16:14