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

% 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;

for j=1:(Ny+1)
y(j) = Ymin + (j-1)*dy;

% Initial condition
for i = 1:(Nx+1)
    for j = 1:(Ny+1)

% Boundary condition
for k=1:(Nt+1)
    u(1,1,k) = 0;
    u(Nx+1,Ny+1,k) = 0;
    t(k) = (k-1)*dt;

% 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));

% Graphical representation of the function
[x,y] = meshgrid(x,y);

for m=1:9
   zlim([0 1]);
   shading interp;

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:

Plot results

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:

  1. 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.

  2. Is there a better way of implementing this algorithm? It runs very slowly for me, even with relatively small amounts of steps.


There was one important mistake in the code, you had +B as the coefficient for the $f_y$ part, not -B. Other than that, I've gone through and removed a lot of for loops by doing vector operations. This speeds the code up quite a bit. I don't have the uzero function so I used mvnpdf.

The stability is quite poor, reducing the time step helps, but a higher-order method for approximating the derivatives might be required if reducing the time step isn't enough, or if the accuracy isn't good enough.

The only for loop required is for the time stepping, since we use results from the previous step to determine the current values, but for the spatial stepping we are using already known values, so it can be done in one go.

Here's my modified code. This should be more accurate and run faster than your code, but it does the same steps:

%// Liouville equation

%// Equation Parameters:
Xmin = -10.0; %// Minimum X
Ymin = -10.0; %// Minimum Y
Xmax = 10.0;  %// Maximum X
Ymax = 10.0;  %// Maximum Y
Tmin=0;       %// Minimum time
Tmax = 5.0;   %// Maximum time
A = 1.0;      %// A parameter
B = 1.0;      %// B parameter

%// Simulation parameters:
Nt = 2000; %// Number of time steps
Nx = 100;  %// Number of X space steps
Ny = 100;  %// Number of X space steps

%// Create time and space vectors with exactly Nt, Nx, Ny elements

%// Make a grid of x and y coordinates    
[X,Y] = ndgrid(x,y);

%// Calculate dt, dx and dy

dtdx = dt/dx;   %// For simplicity
dtdy = dt/dy;

u=zeros(Nx,Ny,Nt); %// Initialise u matrix

%// Initial condition
u(:,:,1)=reshape(mvnpdf([X(:),Y(:)],[1.0 ,0.0],[0.25,0.25]),Nx,Ny);

%// Time stepping algorithm
i=2:Nx-1; %// Interior spatial points
for k=2:Nt
    u(i,j,k) = u(i,j,k-1) - 0.5*A*dtdx*Y(i,j).*(u(i+1,j,k-1)-u(i-1,j,k-1)) ...
        + 0.5*B*dtdy*X(i,j).*(u(i,j+1,k-1)-u(i,j-1,k-1));

%// Animate the solution
for m=1:10:Nt
   shading interp;
  • $\begingroup$ Thanks a lot, David!!! It works very nicely. BTW which higher order method for approximating the derivatives would you recommend? $\endgroup$ Nov 22 '14 at 12:09
  • $\begingroup$ I'm asking cause it seems stability is still a big problem, upon setting B=2, changing the space bounds from 10->13 and the number of steps to 300 for space and 3000 for time, the solution just went to hell at about t=2 :) Ideally I'd want to see the gaussian make a full orbit about the z axis. The B parameter controls the angular frequency, so I'd like the program to work for slightly higher values than 1. $\endgroup$ Nov 22 '14 at 12:16
  • $\begingroup$ I've changed surf to contourf to better see what's going on. With A=1, B=2, X,Y [-10,10], Tmax=5, Nt=5000, Nx=Ny=150 the program holds it together, but chaging Tmax=10 will cause some sort of disco occuring at about t=3. So it seems the amount of time steps is vital... and a sure sign that a more stable algorithm is necessary. $\endgroup$ Nov 22 '14 at 12:36
  • $\begingroup$ A more stable algorithm is probably required. This table: en.wikipedia.org/wiki/… gives you coefficients for higher-order schemes, but you have to be careful with initial conditions. $\endgroup$
    – David
    Nov 22 '14 at 22:26
  • $\begingroup$ So I'd use the first table (central finite difference) on that page for the spacial variables, and the second one (forward) for time, right? And in what sense do you mean to be careful with initial conditions? $\endgroup$ Nov 22 '14 at 22:46

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.