Linear advection: analysis of common numerical methods
In this project we will use Finite Difference schemes to solve (or not!) the linear advection equation in 1 spatial dimension.
We will start by 3 intuitive schemes that actually fail (most of them) and we will explain why. Then, we will show some of the most famous schemes used in industry and academia (which of course work! at least under certain conditions).
The linear advection equation
We are going to numerically solve the following partial differential equation:
`(del u)/(del t) + c (del u)/(del x) = 0`
The parameter c is a real non-zero value (positive or negative) with units of velocity.
The physical meaning of this equation is the transport of the scalar space-dependant quantity "u" at a uniform speed c.
The linear advection is the most basic problem in CFD, which serves as a baseline for lots of numerical schemes for more advanced PDEs such as Euler equations or Navier-Stokes equations.
The exact solution is widely-known, very easy and intuitive. If we have an initial condition given by:
`u(t=0, x) = u_0(x)`
Then the exact solution in function of time is:
`u(t,x) = u_0(x - c t)`
For a spatial domain of finite length, depending on c being positive or negative, we will also need to specify a left or right boundary condition, respectively, as the information coming from left (or right) is unknown at t>0 and needs to be provided.
The fact that the exact solution for arbitrary initial conditions is well-known makes this problem so popular for baseline to compare numerical schemes to be used in more difficult problems.
A type of numerical methods to solve this equation are the Finite Difference methods, which are based on approximating the time and space derivatives by discrete derivatives based on Taylor series.
However, we are going to prove that not any arbitrary choice or combination of discrete derivatives will work, as the resultant scheme needs to be stable and consistent.
In fact, the three most obvious choices (forward, backward, and central derivative) will fail most of times.
Section 1: How NOT to solve the linear advection equation
1. Forward difference scheme (1st order)
This scheme would be the first anybody would try. Both time and space derivatives are 1st order forward difference:
`(u_i^(n+1) - u_i^(n) )/dt - c (u_(i+1)^n - u_i^n )/dx = 0`
With index "i" indicating the spatial point and index "n" the time point. With this scheme we can easily isolate the spatial values of the next time step (n+1) in function of the value from the previous time-step (n):
`u_i^(n+1) = u_i^n - (c*dt)/dx *( u_(i+1)^(n) - u_(i)^(n) )`
Let's solve this in a domain of length 1m, for an initial profile with discontinuities, u being equal to 2 between x= [0.1, 0.3] and equal to 1 in the rest of the domain. The value of c will be positive and equal to 1 m/s.
The following MATLAB code implements the solution:
%% 1. Forward difference scheme (1st order)
% u^(n+1)_i = u^n_i - c*dt/dx*(u^n_(i+1) - u^n_i)
%Parameters
c = 1; %Speed of advection [m/s]
L = 1; %Length of domain [m]
Nx = 100; %Number of points in space
T = 0.1; %Duration of simulation [s]
Nt = 20; %Numer of points in time
%Initial condition
x = linspace(0,L,Nx);
u0 = 1.*ones(1,Nx);
u0(x>=0.1 & x<=0.3) = 2;
u = nan(Nt, Nx);
u(1, :) = u0;
%Left B.C
u(:,1) = u0(1);
%Solve
dx = L/(Nx-1);
dt = T/(Nt-1);
for n = 1:Nt-1
u(n+1,1:Nx-1) = u(n,1:Nx-1) - c*dt/dx.*(u(n,2:Nx) - u(n,1:Nx-1));
end
%Plot
colormapgreen = [ linspace(0, 0, 12)', linspace(0.5, 0.9, 12)', linspace(0.1, 0.5, 12)'];
set(groot,'defaultAxesColorOrder', colormapgreen)
figure();
lstart = plot(x,u(1,:), 'LineWidth', 2);
hold on
plot(x, u(floor(Nt/10):floor(Nt/10):end-floor(Nt/10), :), 'LineWidth', 0.1);
lend = plot(x, u(end, :), 'LineWidth', 2);
title('Forward difference scheme (c>0)');
xlabel('x [m]');
legend([lstart, lend], {'t=0', 't=0.1'});The resulting plot, which shows the profile at 10 equally-spaced times (superposed) looks like this:
This solution is completely wrong. We can observe strong oscillations appearing in the edges of the step and propagating throught the domain. This behavior is a classical unstable one.
The reason for this scheme to be unstable is that in the physical problem the information is propagating from left to right, and a forward spatial derivative is using information that has not been reached yet.
Conversely, for c<0 there's apriori no reason for the scheme to be wrong, let's see:
% Make c<0 ------------------
c = -1;
%Initial condition
u0 = 1.*ones(1,Nx);
u0(x>=0.6 & x<=0.8) = 2;
u = nan(Nt, Nx);
u(1, :) = u0;
%Right B.C
u(:,end) = u0(end);
%Solve
dx = L/(Nx-1);
dt = T/(Nt-1);
for n = 1:Nt-1
u(n+1,1:Nx-1) = u(n,1:Nx-1) - c*dt/dx.*(u(n,2:Nx) - u(n,1:Nx-1));
end
%Plot
figure();
lstart = plot(x,u(1,:), 'LineWidth', 2);
hold on
plot(x, u(floor(Nt/10):floor(Nt/10):end-floor(Nt/10), :), 'LineWidth', 0.1);
lend = plot(x, u(end, :), 'LineWidth', 2);
ylim([0.5 2.5]);
title('Forward difference scheme (c<0)');
xlabel('x [m]');
legend([lstart, lend], {'t=0', 't=0.1'});We can see that this time the solution is stable, propagating left in the correct speed (with some diffusion being introduced by the discretization error)
2. Backward difference scheme (1st order)
This scheme is exactly the same than the previous one but using a spatial backward difference:
`u_i^(n+1) = u_i^n - (c*dt)/dx *( u_(i)^(n) - u_(i-1)^(n) )`
For the same reasons than before, this scheme will work for c>0 but not for c<0. Let's see the case c>0, this time with a longer time-span:
%% 2. Backward difference scheme (1st order)
% u^(n+1)_i = u^n_i - c*dt/dx*(u^n_(i+1) - u^n_i)
%Parameters
c = 1; %Speed of advection [m/s]
L = 1; %Length of domain [m]
Nx = 100; %Number of points in space
T = 0.3; %Duration of simulation [s]
Nt = 100; %Numer of points in time
%Initial condition
x = linspace(0,L,Nx);
u0 = 1.*ones(1,Nx);
u0(x>=0.1 & x<=0.3) = 2;
u = nan(Nt, Nx);
u(1, :) = u0;
%Left B.C
u(:,1) = u0(1);
%Solve
dx = L/(Nx-1);
dt = T/(Nt-1);
for n = 1:Nt-1
u(n+1,2:Nx) = u(n,2:Nx) - c*dt/dx*(u(n,2:Nx) - u(n,1:Nx-1));
end
%Plot
colormapblue = [linspace(0, 0, 12)', linspace(0, 1, 12)', linspace(1, 1, 12)'];
set(groot,'defaultAxesColorOrder', colormapblue)
figure();
lstart = plot(x,u(1,:), 'LineWidth', 2);
hold on
plot(x, u(floor(Nt/10):floor(Nt/10):end-floor(Nt/10), :), 'LineWidth', 0.1);
lend = plot(x, u(end, :), 'LineWidth', 2);
ylim([0.5 2.5]);
title('Backward difference scheme');
xlabel('x [m]');
legend([lstart, lend], {'t=0', 't=0.3'});
We can see it works as expected.
3. Central difference scheme (2nd order)
Some people would try to get a higher-order of accuracy by using the central difference for the spatial 1st derivative (which is of second order in space). However, this would terribly fail for both cases c>0 and c<0 as it is unwisely using information from both left and right side, and one or the other are physically not yet propagated. Let's see it:
%% 3. Central difference scheme
% u^(n+1)_i = u^n_i - 0.5*c*dt/dx*(u^n_(i+1) - u^n_(i-1))
%Parameters
c = 1; %Speed of advection [m/s]
L = 1; %Length of domain [m]
Nx = 100; %Number of points in space
T = 0.1; %Duration of simulation [s]
Nt = 20; %Numer of points in time
%Initial condition
x = linspace(0,L,Nx);
u0 = 1.*ones(1,Nx);
u0(x>=0.1 & x<=0.3) = 2;
u = nan(Nt, Nx);
u(1, :) = u0;
%Left B.C
u(:,1) = u0(1);
%Solve
dx = L/(Nx-1);
dt = T/(Nt-1);
for n = 1:Nt-1
u(n+1,2:Nx-1) = u(n,2:Nx-1) - 0.5*c*dt/dx*(u(n,3:Nx) - u(n,1:Nx-2));
end
%Plot
colormapred = [ linspace(0.5, 1, 7)', linspace(0, 0, 7)', linspace(0, 0, 7)'];
set(groot,'defaultAxesColorOrder', colormapred)
figure();
lstart = plot(x,u(1,:), 'LineWidth', 2);
hold on
plot(x, u(floor(Nt/10):floor(Nt/10):end-floor(Nt/10), :), 'LineWidth', 0.1);
lend = plot(x, u(end, :), 'LineWidth', 2);
title('Central difference scheme');
xlabel('x [m]');
legend([lstart, lend], {'t=0', 't=0.1'});
Again, unstabilities are produced in the edges of the step, causing the solution to fail.
We have seen in this section that choosing the most obvious discretizations for this problem will make us fail miserably for unstability reasons.
Fortunately, there exists wise choices of discretizations that wise people have found in the past. These schemes are stable always or under certain conditions, and we can find a variety of orders of accuracy.
Section 2: Literature schemes (or how to solve the linear advection equation)
In this section we are going to show some of the most well-known schemes that appear in the literature. Most of these are used in commercial and research CFD codes.
4. Upwind scheme
This scheme is explicit, conditionally-stable, and of 1st order in both space and time.
The idea of this scheme is to use the forward or backward difference depending on the sign of c (intelligent choice, right?).
A compact form of the scheme reads as follows:
`u_i^(n+1) = u_i^n - (c*dt)/dx *( max(c,0)*( u_(i)^(n) - u_(i-1)^(n) ) + min(c,0)*( u_(i+1)^(n) - u_(i)^(n) ) )`
As both c>0 and c<0 are considered, we need to provide both left B.C. and right B.C.
Let's prove that it works both for c>0 and c<0:
%% 4. Upwind scheme
% u^(n+1)_i = u^n_i - dt/dx*( max(c,0)*(u^n_i - u^n_(i-1)) + min(c,0)*(u^n_(i+1) - u^n_i))
% 5.1 Case c>0
%Parameters
c = 1; %Speed of advection [m/s]
L = 1; %Length of domain [m]
Nx = 200; %Number of points in space
T = 0.3; %Duration of simulation [s]
Nt = 100; %Numer of points in time
%Initial condition
x = linspace(0,L,Nx);
u0 = 1.*ones(1,Nx);
u0(x>=0.4 & x<=0.6) = 2;
u = nan(Nt, Nx);
u(1, :) = u0;
%Left B.C
u(:,1) = u0(1);
%Right B.C
u(:,end) = u0(end);
%Solve
dx = L/(Nx-1);
dt = T/(Nt-1);
for n = 1:Nt-1
u(n+1,2:Nx-1) = u(n,2:Nx-1) ...
- dt/dx*( max(c,0)*(u(n,2:Nx-1) - u(n,1:Nx-2) ...
+ min(c,0)*(u(n,3:Nx) - u(n,2:Nx-1) )));
end
%Plot
colormapora = [ linspace(1, 0.929, 7)', linspace(0.4118, 0.6941, 7)', linspace(0.1608, 0.1255, 7)'];
set(groot,'defaultAxesColorOrder', colormapora)
figure();
lstart = plot(x,u(1,:), 'LineWidth', 2);
hold on
plot(x, u(floor(Nt/5):floor(Nt/5):end-floor(Nt/5), :), 'LineWidth', 0.1);
lend = plot(x, u(end, :), 'LineWidth', 2);
ylim([0.5 2.5]);
title('Upwind scheme (c>0)');
xlabel('x [m]');
legend([lstart, lend], {'t=0', 't=0.3'});
% 5.2 Case c<0
c = -1;
%Initial condition
u = nan(Nt, Nx);
u(1,:) = u0;
%Left B.C
u(:,1) = u0(1);
%Right B.C
u(:,end) = u0(end);
%Solve
dx = L/(Nx-1);
dt = T/(Nt-1);
for n = 1:Nt-1
u(n+1,2:Nx-1) = u(n,2:Nx-1) ...
- dt/dx*( max(c,0)*(u(n,2:Nx-1) - u(n,1:Nx-2)) ...
+ min(c,0)*(u(n,3:Nx) - u(n,2:Nx-1)) );
end
%Plot
figure();
lstart = plot(x,u(1,:), 'LineWidth', 2);
hold on
plot(x, u(floor(Nt/5):floor(Nt/5):end-floor(Nt/5), :), 'LineWidth', 0.1);
lend = plot(x, u(end, :), 'LineWidth', 2);
ylim([0.5 2.5]);
title('Upwind scheme (c<0)');
xlabel('x [m]');
legend([lstart, lend], {'t=0', 't=0.3'});The upwind scheme is stable depending on the step size of time and space discretizations, according to the following criteria:
`(abs(c)dt)/dx <= 1`
If the last condition is met, then the upwind scheme is stable and the solution converges in 1st order on space and time.
This condition is called Courant-Friedrichs-Lewy condition (or CFL condition for short) and it is very similar for most the explicit schemes.
5. Lax-Friedrichs scheme
This scheme is explicit, conditionally stable and first order accurate both in space and time. It has a higher dissipative effect than the Upwind scheme
The scheme is:
`u_(i)^(n+1) = 1/2(u_(i+1)^(n) + u_(i-1)^(n)) - 1/2(c*dt)/dx(u_(i+1)^(n) + u_(i-1)^(n))`
%% 5. Lax-Friedrichs scheme
%Parameters
c = 1; %Speed of advection [m/s]
L = 1; %Length of domain [m]
Nx = 200; %Number of points in space
T = 0.3; %Duration of simulation [s]
Nt = 100; %Numer of points in time
%Initial condition
x = linspace(0,L,Nx);
u0 = 1.*ones(1,Nx);
u0(x>=0.1 & x<=0.3) = 2;
u = nan(Nt, Nx);
u(1, :) = u0;
%Left B.C
u(:,1) = u0(1);
%Right B.C
u(:,end) = u0(end);
%Solve
dx = L/(Nx-1);
dt = T/(Nt-1);
for n = 1:Nt-1
u(n+1,2:Nx-1) = 0.5*( u(n,3:Nx) + u(n,1:Nx-2) ) ...
- 0.5*c*dt/dx*( u(n,3:Nx) - u(n,1:Nx-2) );
end
%Plot
colormappurp = [ linspace(0.7176, 1, 7)', linspace(0.2745, 0, 7)', linspace(1, 1, 7)'];
set(groot,'defaultAxesColorOrder', colormappurp)
figure();
lstart = plot(x,u(1,:), 'LineWidth', 2);
hold on
plot(x, u(floor(Nt/5):floor(Nt/5):end-floor(Nt/5), :), 'LineWidth', 0.1);
lend = plot(x, u(end, :), 'LineWidth', 2);
ylim([0.5 2.5]);
title('Lax-Friedrichs scheme');
xlabel('x [m]');
legend([lstart, lend], {'t=0', 't=0.3'});
We can see the solution contains more diffusion and is a bit step-wised, but still quite accurate.
The Lax-Friedrichs scheme has the same CFL condition than the Upwind:
`(abs(c)dt)/dx <= 1`
6. Lax-Wendroff scheme
This scheme is also explicit, conditionally stable, but 2nd order accurate both in space and time, which gives it a strong advantage over the previous ones. This means that the error is `O(Deltax^2 + Deltat^2)` instead of `O(Deltax + Deltat)`
%% 6. Lax-Wendroff scheme
%Parameters
c = 1; %Speed of advection [m/s]
L = 1; %Length of domain [m]
Nx = 200; %Number of points in space
T = 0.3; %Duration of simulation [s]
Nt = 100; %Numer of points in time
%Initial condition
x = linspace(0,L,Nx);
u0 = 1.*ones(1,Nx);
u0(x>=0.1 & x<=0.3) = 2;
u = nan(Nt, Nx);
u(1, :) = u0;
%Left B.C
u(:,1) = u0(1);
%Right B.C
u(:,end) = u0(end);
%Solve
dx = L/(Nx-1);
dt = T/(Nt-1);
for n = 1:Nt-1
u(n+1,2:Nx-1) = u(n, 2:Nx-1) ...
- 0.5*c*dt/dx*( u(n,3:Nx) - u(n,1:Nx-2) ) ...
+ 0.5*(c*dt/dx)^2*( u(n,3:Nx) - 2*u(n,2:Nx-1) + u(n,1:Nx-2) );
end
%Plot
colormapblue = [linspace(0, 0, 12)', linspace(0, 1, 12)', linspace(1, 1, 12)'];
set(groot,'defaultAxesColorOrder', colormapblue)
figure();
lstart = plot(x,u(1,:), 'LineWidth', 2);
hold on
plot(x, u(floor(Nt/5):floor(Nt/5):end-floor(Nt/5), :), 'LineWidth', 0.1);
lend = plot(x, u(end, :), 'LineWidth', 2);
ylim([0.5 2.5]);
title('Lax-Wendroff scheme');
xlabel('x [m]');
legend([lstart, lend], {'t=0', 't=0.3'});
We can see that the solution has less dissipative effects (it is more squared) but it contains some perturbations in the edges of the step. This is the disadvantage of the LW scheme, as discontinuities in the solution introduce some not-so-beautiful perturbations.
The LW scheme CFL condition is also `(abs(c)dt)/dx <= 1`
7. Leapfrog scheme
This scheme is explicit, and it has a big peculiarity (apart from its funny name): it uses 2 time steps to calculate the following one:
`u_(i)^(n+1) = u_(i)^(n-1) - (c*dt)/dx*(u_(i+1)^(n)-u_(i-1)^(n))`
This means that in order to calculate the solution on the next time step (n+1) we need to use the solution in (n) and in (n-1).
The consequence is that we will need to use a 1 time-step scheme for the initial iteration to find out n=2 from n=1 and then use them both to calculate n=3 and the rest using the Leapfrog scheme.
In our case, we choose the Upwind scheme to initialize the problem:
%% 7. Leapfrog scheme
%Parameters
c = 1; %Speed of advection [m/s]
L = 1; %Length of domain [m]
Nx = 200; %Number of points in space
T = 0.3; %Duration of simulation [s]
Nt = 100; %Numer of points in time
%Initial condition
x = linspace(0,L,Nx);
u0 = 1.*ones(1,Nx);
u0(x>=0.1 & x<=0.3) = 2;
u = nan(Nt, Nx);
u(1, :) = u0;
%Left B.C
u(:,1) = u0(1);
%Right B.C
u(:,end) = u0(end);
%Solve
dx = L/(Nx-1);
dt = T/(Nt-1);
%1. First iteration to get u(2,:) -> Use upwind
u(2,2:Nx-1) = u(1,2:Nx-1) ...
- dt/dx*( max(c,0)*(u(1,2:Nx-1) - u(1,1:Nx-2) ...
+ min(c,0)*(u(1,3:Nx) - u(1,2:Nx-1) )));
%2. Rest of iterations
for n = 2:Nt-1
u(n+1,2:Nx-1) = u(n-1, 2:Nx-1) ...
- c*dt/dx*( u(n,3:Nx) - u(n,1:Nx-2) );
end
%Plot
colormapfrog = [linspace(0.4667, 0.3922, 12)', linspace(0.6745, 0.8314, 12)', linspace(0.1882, 0.0745, 12)'];
set(groot,'defaultAxesColorOrder', colormapfrog)
figure();
lstart = plot(x,u(1,:), 'LineWidth', 2);
hold on
plot(x, u(floor(Nt/5):floor(Nt/5):end-floor(Nt/5), :), 'LineWidth', 0.1);
lend = plot(x, u(end, :), 'LineWidth', 2);
ylim([0.5 2.5]);
title('Leapfrog scheme');
xlabel('x [m]');
legend([lstart, lend], {'t=0', 't=0.3'});We can see that the solution has lots of perturbations that quickly propagate backwards. This effects may be changed by initializing with different schemes.
As a summary of this 4 explicit schemes, let's watch an animation of them all:
The code used to generate the animation is:
%% Animation
%Parameters
c = 1; %Speed of advection [m/s]
L = 1; %Length of domain [m]
Nx = 200; %Number of points in space
T = 0.6; %Duration of simulation [s]
Nt = 200; %Numer of points in time
%Initial condition
x = linspace(0,L,Nx);
u0 = 1.*ones(1,Nx);
u0(x>=0.1 & x<=0.3) = 2;
dx = L/(Nx-1);
dt = T/(Nt-1);
% 1. Upwind
u_uw = nan(Nt, Nx);
u_uw(1, :) = u0;
%Left B.C
u_uw(:,1) = u0(1);
%Right B.C
u_uw(:,end) = u0(end);
for n = 1:Nt-1
u_uw(n+1,2:Nx-1) = u_uw(n,2:Nx-1) ...
- dt/dx*( max(c,0)*(u_uw(n,2:Nx-1) - u_uw(n,1:Nx-2) ...
+ min(c,0)*(u_uw(n,3:Nx) - u_uw(n,2:Nx-1) )));
end
% 2. Lax-Friedrichs
u_lf = nan(Nt, Nx);
u_lf(1, :) = u0;
%Left B.C
u_lf(:,1) = u0(1);
%Right B.C
u_lf(:,end) = u0(end);
for n = 1:Nt-1
u_lf(n+1,2:Nx-1) = 0.5*( u_lf(n,3:Nx) + u_lf(n,1:Nx-2) ) ...
- 0.5*c*dt/dx*( u_lf(n,3:Nx) - u_lf(n,1:Nx-2) );
end
% 3. Lax-Wendroff
u_lw = nan(Nt, Nx);
u_lw(1, :) = u0;
%Left B.C
u_lw(:,1) = u0(1);
%Right B.C
u_lw(:,end) = u0(end);
for n = 1:Nt-1
u_lw(n+1,2:Nx-1) = u_lw(n, 2:Nx-1) ...
- 0.5*c*dt/dx*( u_lw(n,3:Nx) - u_lw(n,1:Nx-2) ) ...
+ 0.5*(c*dt/dx)^2*( u_lw(n,3:Nx) - 2*u_lw(n,2:Nx-1) + u_lw(n,1:Nx-2) );
end
% 4. Leapfrog
u_lpf = nan(Nt, Nx);
u_lpf(1, :) = u0;
%Left B.C
u_lpf(:,1) = u0(1);
%Right B.C
u_lpf(:,end) = u0(end);
u_lpf(2,2:Nx-1) = u_lpf(1,2:Nx-1) ...
- dt/dx*( max(c,0)*(u_lpf(1,2:Nx-1) - u_lpf(1,1:Nx-2) ...
+ min(c,0)*(u_lpf(1,3:Nx) - u_lpf(1,2:Nx-1) )));
for n = 2:Nt-1
u_lpf(n+1,2:Nx-1) = u_lpf(n-1, 2:Nx-1) ...
- c*dt/dx*( u_lpf(n,3:Nx) - u_lpf(n,1:Nx-2) );
end
fig = figure('Position', [200, 200, 1024, 768]);
for i = 1:Nt
subplot(2,2,1);
set(groot,'defaultAxesColorOrder', colormapora)
plot(x, u_uw(i,:), 'LineWidth', 2);
ylim([0.5 2.5]);
title('Upwind')
xlabel('x [m]')
subplot(2,2,2);
set(groot,'defaultAxesColorOrder', colormappurp)
plot(x, u_lf(i,:), 'LineWidth', 2);
ylim([0.5 2.5]);
title('Lax-Friedrichs')
xlabel('x [m]')
subplot(2,2,3);
set(groot,'defaultAxesColorOrder', colormapblue)
plot(x, u_lw(i,:), 'LineWidth', 2);
ylim([0.5 2.5]);
title('Lax-Wendroff')
xlabel('x [m]')
subplot(2,2,4);
set(groot,'defaultAxesColorOrder', colormapfrog)
plot(x, u_lpf(i,:), 'LineWidth', 2);
ylim([0.5 2.5]);
title('Leapfrog')
xlabel('x [m]')
pause(3*dt);
M(i) = getframe(fig);
end
In this project we have focused only on explicit methods which are conditionally-stable most of times, but there exists also lots of implicit methods which have the advantage of being unconditionally stable. However, they are computationally more costly as they require to solve a linear system of Nx equations and Nx variables for each time step, so Nt times.
As a conclusion: solving CFD problems can become a real headache if we don't use the right scheme!