traffic.m

% traffic - Program to solve the generalized Burger  
% equation for the traffic at a stop light problem
clear;  help traffic; % Clear memory and print header
N = input('Enter the number of grid points - ');
L = 400;    % System size (meters)
h = L/N;    % Grid spacing for periodic boundary conditions
v_max = 25;    % Maximum car speed (m/s)
fprintf('Suggested timestep is %g\n',h/v_max);
tau = input('Enter time step (tau) - ');
coeff = tau/(2*h);   % Coefficient used by the schemes
coefflw = tau^2/(2*h^2);  % Coefficient used by Lax-Wendroff
%%%%% Initial condition is a square pulse		%%%%% 
%%%%% from x = -L/4 to x = 0               %%%%%
rho_max = 1;                   % Maximum density
Flow_max = 1/4*rho_max*v_max;  % Maximum Flow
rho = zeros(1,N);
for i=N/4:(N/2-1)
  rho(i) = rho_max;    % Max density in the square pulse
end
rho(N/2) = rho_max/2;  % Try running without this line
%%%%% Set loop and plot variables %%%%%
iplot = 1;
xplot = ((1:N)-1/2)*h - L/2;  % Record x scale for plot
rplot(:,1) = rho(:);          % Record the initial state
tplot(1) = 0;
nstep = ceil((L/4)/(v_max*tau)); % Stop after last car moves
%nstep = 5*nstep;  % Go further for last plots %******%
fprintf('Number of iterations = %g\n',nstep);
%%%%% MAIN LOOP %%%%%
ip(1:N) = (1:N)+1;  ip(N) = 1;   % ip = i+1 with periodic b.c.
im(1:N) = (1:N)-1;  im(1) = N;   % im = i-1 with periodic b.c.
for istep=1:nstep
  Flow = rho .* (v_max*(1 - rho/rho_max));
  %%% FTCS method %%%
  rho_new(1:N) = rho(1:N) - coeff*(Flow(ip)-Flow(im));
  %%% Lax method %%%
%  rho_new(1:N) = .5*(rho(ip)+rho(im)) ...
%                  - coeff*(Flow(ip)-Flow(im));
  %%% Lax-Wendroff method %%%
%  cp = v_max*(1 - (rho(ip)+rho(1:N))/rho_max);
%  cm = v_max*(1 - (rho(1:N)+rho(im))/rho_max);
%  rho_new(1:N) = rho(1:N) - coeff*(Flow(ip)-Flow(im)) ...
%   + coefflw*(cp.*(Flow(ip)-Flow(1:N)) ...
%             - cm.*(Flow(1:N)-Flow(im)));
  %%% ----------- %%%
  rho = rho_new;
  iplot = iplot+1;
  rplot(:,iplot) = rho(:);   % Record density for plot
  tplot(iplot) = tau*istep;  % Record time for plot
  plot(xplot,rho,'-',xplot,Flow/Flow_max,'--');
  title('Density (solid) and flow (dashed) versus position');
end
mesh(rplot,[-90 30])  % Wire-mesh plot of density
xlabel('Position'); ylabel('Time');
title('Density versus x and t');
pause;    % Pause between plots; strike any key to continue
% Use rot90 function to graph t vs x since
% contour(rplot) graphs x vs t.
levels = 0:.1:1;    % Contour levels in plot
cs = contour(rot90(rplot),levels,xplot,tplot); 
xlabel('Position');  ylabel('Time');              
clabel(cs,0:0.1:1);       % Put labels on contour levels            
title('Density contours');