Why is this implementation of the Lax-Friedrichs scheme so dissipative?

Question

I needed a numerical solver for the advection equation in 1D. So I decided to implement the Lax-Friedrichs scheme in Matlab like so:

function [output] = advectionSolver(velocity,initial_Condition,boundary_Condition,xmin,xmax,tmin,tmax)
  N = length(initial_Condition);
  M = length(boundary_Condition);
  dx = (xmax-xmin)/(N-1);
  dt = (tmax-tmin)/(M-1);
  data = zeros(M,N+2);
  data(1,2:N+1) = initial_Condition;
  data(:,2) = boundary_Condition;
  data(:,N+1) = boundary_Condition;
  CFL = 0.5*velocity*dt/dx;
  for t = 1:M-1
    data(:,N+2) = data(:,2);
    data(:,1) = data(:,N+1);
    for x = 2:N+1
        data(t+1,x) = 0.5*(data(t,x+1)+data(t,x-1))-CFL*(data(t,x+1)-data(t,x-1));
    end
  end
  output = data(:,2:N+1);
end

This is for a problem with periodic boundary conditions. It takes as input a constant velocity, two arrays for the conditions and the bounds of the domain.

Now the problem I stumble upon is that this seems very dissipative. For example when inputting $100$ values of the function $\frac{1}{\sqrt{2\pi\cdot 0.1^2}}e^{-\frac{(x-0.5)^2}{0.1^2}}$ where $x$ is between $0$ and $1$ and taking $100$ zero values for the boundary condition with $t$ between $0$ and $0.5$, I get the following result:

Here, the blue graph is at $t=0$ and the green graph is at $t=0.5$. This is not what I would expect to happen, indeed the graph moves to the right as the time increases, but why does it dissipate so much? I would expect the graph to move to the right and have the same height everywhere or maybe decrease a little bit due to numerical errors, but the final graph is almost half the size of the first one.

I would like to know why this is exactly, are the numerical errors indeed that large? Is this an effect obtained by using the Lax-Friedrichs scheme? Or is there an error in my reasoning?

Answer 1

By discretizing a PDE, you will introduce an error, which then again might introduce numerical diffusion, which is what you observe here.

If you look at the Taylor expansion that you use to construct the Lax-Friedrichs scheme, you see that you cut off the expansion at the second order term, introducing an error of the form $h^2\frac{\partial^2 u}{\partial x^2}$, which looks like a diffusion term, and as you might expect this is what causes the numerical diffusion you observe as well.

This can be controlled however by increasing the number of point you use in $x$, but note that you might need to use smaller $\Delta t$ as well in order to preserve stability.