I'm trying to solve the advection equation $$m_t+(\alpha m)_x=0$$ with $m(0,\cdot)=m_0$ numerically using the first order Godunov scheme. Hence I write $$m_j^{i+1}=m_j^i-\frac{dt}{dx}(m_{j+1/2}^i\alpha_{j+1/2}^i-m_{j-1/2}^i\alpha_{j-1/2}^i)$$ where $m_{j+1/2}^i=m_{j+1}^i$ if $\alpha_{j+1/2}^i<0$ and $m_j^i$ otherwise.
Here's a Matlab code with $\alpha=1/2$ everywhere:
%% parameters
N=99; M=99;
%Allocate memory
x=linspace(0,1,M+1);
mk=zeros(N+1,M+1); a_ph=ones(N,M)*0.5;
%% Set initial conditions
sigma=0.08;
mk(1,:)=1/sqrt(2*pi*sigma^2)*exp(-(x-0.5).^2/(2*sigma^2));
q=dt/h;
for i=1:N
ai=ak_ph(i,2:M); aim=a_ph(i,1:M-1); %local var
B=diag([0,ai.*(ai>=0),0]) + diag([0,ai.*(ai<0)],1) ...
- diag([0,aim.*(aim<0),0]) - diag([aim.*(aim>=0),0],-1);
mk(i+1,:)=(eye(M+1)-q*B)*mk(i,:)';
end
[X1,X2]=meshgrid(x1,x2);
surf(X1,X2,mk','FaceColor','interp','EdgeColor','none')
When I compute the solution this is what I get:
I expect a wave travelling at a uniform velocity but the amplitude gets smaller in the meantime which is strange. There are also ripples at the peak of the wave. Does anyone have a suggestion?