Backup of this question:
Plotting Graph - 1D Advection
clear all; clc; close all; % clear workspace and editor and close figures, respectively
hold on
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
xlength = 2; % length of computational domain
n = 1000; % number of grid points
h = xlength / (n-1); % dx
cfl = 0.9 ; % cfl = U dx/dt with stability between 0 and 1
U = 1; % advective speed
dt = h * cfl / U; % dt - which is now function of cfl, h & dx
tout = 1; % desired output time
time = 0;
x = zeros (1,n); % initialising grid points
fn = zeros (1,n); % initialising solution of first-order Upwind scheme
fnlw = zeros (1,n); % initialising solution of Lax-Wendroff scheme
fnlf = zeros (1,n); % initialising solution of Lax-Friedrichs scheme
f = zeros (1,n);
flw = zeros (1,n);
freal = zeros (1,n);
x(1) = 0;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
for time = [0 0.5 1]
z = 0.75 * exp(-(((x-0.5)/0.1)).^2);
for i = 2:n % making the mesh
x(i) = x(i-1) + h;
end
for i = 1:n
f(i) = 0.75 * exp(-(((x(i)-0.5)/0.1)).^2); % first-order Upwind scheme
flw(i) = 0.75 * exp(-(((x(i)-0.5)/0.1)).^2); % Lax-Wendroff scheme
flf(i) = 0.75 * exp(-(((x(i)-0.5)/0.1)).^2); % Lax-Friedrichs scheme
end
nt = time/dt;
for k = 1:nt
for i = 2:n % first-order Upwind scheme
flux = U * (f(i)-f(i-1));
fn(i) = f(i)-(dt/h) * flux;
end
fn(1) = fn(n); % boundary condition for first-order Upwind scheme
f = fn; % boundary condition for first-order Upwind scheme
for i = 2:n-1 % Lax-Wendroff scheme
l0 = (dt/(2*h)) * U *(flw(i+1)-flw(i-1));
h0 = (dt^2/(2*h^2)) * U.^2 * (flw(i+1)-(2*flw(i))+flw(i-1));
fnlw(i) = flw(i) - l0 + h0;
end
fnlw(1) = fn(n); % boundary condition for Lax-Wendroff scheme
fnlw(n) = fn(1); % boundary condition for Lax-Wendroff scheme
flw = fnlw; % boundary condition for Lax-Wendroff scheme
for i = 2:n-1 % Lax-Friedrichs scheme
fnlf(i) = 0.5 * (flf(i-1)+flf(i+1))-(dt/(2*h)) * U * (flf(i+1)-flf(i-1));
end
fnlf(1) = fn(n); % boundary condition for Lax-Friedrichs scheme
fnlf(n) = fn(1); % boundary condition for Lax-Friedrichs scheme
fnlf(i) = fn(i); % boundary condition for Lax-Friedrichs scheme
fnlf(n) = fnlf(1); % boundary condition for Lax-Friedrichs scheme
flf = fnlf; % boundary condition for Lax-Friedrichs scheme
time = nt*dt;
freal(:) = 0.75 * exp(-(((x(:)-0.5)/0.1)).^2);
for i = 2:n-1
freal(i) = 0.75 * exp(-(((x(i)-0.5-U*time)/0.1)).^2);
end
freal(1) = freal(n);
freal(n) = freal(1);
end
figure(1)
plot (x,z,'b') % plot initial
hold on
plot (x,freal,'g') % plot real
hold on
plot (x,f,'k') % plot Upwind
hold on
figure(2)
plot (x,z,'b') % plot initial
hold on
plot (x,freal,'g') % plot real
hold on
plot (x,flf,'r') % plot Lax-Friedrichs
hold on
figure(3)
plot (x,z,'b') % plot initial
hold on
plot (x,freal,'g') % plot real
hold on
plot (x,flw,'c') % plot Lax-Wendroff
hold on
title('\color{black}\fontsize{12}\bf Concentration Profile of Plume')
legend('Initial','Real','First-Order Upwind','Lax-Friedrichs','Lax-Wendroff')
set (gca, 'fontsize', 14)
xlabel('\color{black}\fontsize{12}\bf Displacement, x')
xlim([0,2])
ylabel('\color{black}\fontsize{12}\bf Concentration, c')
ErrFOUp = h * sqrt(sum((fn-freal).^2)) % error in first-order Upwind scheme
ErrLaxW = h * sqrt(sum((fnlw-freal).^2)) % error in Lax-Wendroff scheme
ErrLaxF = h * sqrt(sum((fnlf-freal).^2)) % error in Lax-Friedrichs scheme
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
