%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% Integrate the one-dimensional
% advection equation
%
% du/dt + c * du/dx = 0
%
% using the leapfrog method and
% assuming periodic boundary conditions
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clear; % Clear the memory
clf; % Clear the graphics window
L = 10000; % Domain size (km)
M = 100; % Number of grid steps
DX = L/M; % Grid length (km)
T = 5; % Time length of integration (days)
N = 100; % Number of time steps
DT = T/N; % Time step (days)
dx = 1000 * DX; % Grid length (m)
dt = (24*60*60) * DT; % Time step (s)
c = 10 ; % Constant phase speed (m/s)
lambda = c*dt/dx; % Courant number
fprintf(' CFL number: %g \n', lambda)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Arrays for the solution
un = zeros(1,M); % Solution at time n*dt
unm1 = zeros(1,M); % Solution at time (n-1)*dt
unp1 = zeros(1,M); % Solution at time (n+1)*dt
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Indices for periodic domain
M0 = [1:M ]; % Indices for centre points
MM = [M 1:M-1]; % Indices for points to left
MP = [2:M 1 ]; % Indices for points to right
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Define the initial conditions
k = 2; % Number of wavelengths
un = 10*sin(2*pi*k*M0/M);
plot(M0,un,'b',M0,un,'rx');
axis([0 M -12 +12]);
hold on; drawnow; pause
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% MAIN LOOP
% First time step is forward
dudx = un(MP)-un(MM) ;
unp1 = un - 0.5*(lambda)*dudx;
% shuffle the solutions
unm1 = un;
un = unp1;
plot(M0,un); drawnow; pause
for n=2:N
dudx = un(MP)-un(MM) ;
unp1 = unm1 - lambda*dudx;
% shuffle the solutions
unm1 = un;
un = unp1;
plot(M0,un); drawnow; pause(0.1)
end
hold off; % Release hold on graphics window
return
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%