Follow up to my previous question: MATLAB: solving 1st order hyperbolic equation in 2 spacial dimensions
The equation I'm solving has the form: $$f_t + A y f_x - B x f_y =0$$
I wrote the following code according to the comments from the previous question:
% Liouville equation clear; % Equation Parameters: Xmin = -10.0; % Minimum X Ymin = -10.0; % Minimum Y Xmax = 10.0; % Maximum X Ymax = 10.0; % Maximum Y Tmax = 1.0; % Maximum time A = 1.0; % A parameter B = 1.0; % B parameter % Simulation parameters: Nt = 1000; % Number of time steps dt = Tmax/Nt; Nx = 200; % Number of X space steps dx = (Xmax-Xmin)/Nx; Ny = 200; % Number of X space steps dy = (Ymax-Ymin)/Ny; dtdx = dt/dx; % For simplicity dtdy = dt/dy; % Filling the x,y steps for i=1:(Nx+1) x(i) = Xmin + (i-1)*dx; end for j=1:(Ny+1) y(j) = Ymin + (j-1)*dy; end % Initial condition for i = 1:(Nx+1) for j = 1:(Ny+1) u(i,j,1)=uzero(x(i),y(j),1.0,0.0,0.25,0.25); end end % Boundary condition for k=1:(Nt+1) u(1,1,k) = 0; u(Nx+1,Ny+1,k) = 0; t(k) = (k-1)*dt; end % Time stepping algorithm for k=1:Nt % Time loop for i=2:Nx % Space loop for j=2:Ny % Space loop u(i,j,k+1) = u(i,j,k)- 0.5*A*dtdx*y(j)*(u(i+1,j,k)-u(i-1,j,k)) - 0.5*B*dtdy*x(i)*(u(i,j+1,k)-u(i,j-1,k)); end end end % Graphical representation of the function [x,y] = meshgrid(x,y); for m=1:9 subplot(3,3,m); surf(x,y,u(:,:,round(m*Nt/9))); zlim([0 1]); shading interp; end
The uzero function is a gaussian. I'm using one centered at $(1,0)$ with standard deviations $0.25$. This is the result of the plots at 9 different points in time:
The theoretical behaviour of this should be this: (http://upload.wikimedia.org/wikipedia/commons/d/d6/DisplacedGaussianWF.gif). Clearly this isn't the case, the Gaussian loses its shape (becomes flat). I've checked the scheme and I don't see any mistakes in implementation. So my questions are:
What's wrong? Is it some mistake in the code I didn't catch or is the used solution scheme too primitive? Even small changes in the paramaters sometimes produce completely bizzare results, so it looks to me like a stability problem.
Is there a better way of implementing this algorithm? It runs very slowly for me, even with relatively small amounts of steps.