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