3.1)
The discrete problem reads
We perturb the solution 
by a small perturbation
by a small perturbation 
which oscillates in space with a single arbitraty wavenumber k, and determine whether the amplitude of perturbation 
grows or decays in time. Substituting the decomposition 
into (1) and knowing that 
is a solution, the evolution of ξ must also satisfy (1)
grows or decays in time. Substituting the decomposition into (1) and knowing that is a solution, the evolution of ξ must also satisfy (1)
We define a growth factor
such that the perturbation grows if 
and decays if 
. Dividing (2) by 
we obtain
and decays if . Dividing (2) by we obtain
The magnitude of the complex growth factor is then
which after re-arrangement leads to the condition of stability
The dimensionless group 
is known as Courant–Friedrichs–Lewy (CFL) number. Note that the critical value 
corresponds to the situation when u is advected by exactly one grid spacing 
after one time step 
.
is known as Courant–Friedrichs–Lewy (CFL) number. Note that the critical value corresponds to the situation when u is advected by exactly one grid spacing after one time step .3.2)
If we compare Lax-Friedrichs scheme
and explicit Euler scheme
we note that (3) is replacing 
by an average of the neighbouring points
by an average of the neighbouring points
in order to suppress oscillations. We can compute the error of the approximation, ε, by expanding 
with Taylor series from
with Taylor series from 
The odd-order terms cancel out, so we can express the error as
We can therefore expand the Lax-Friedrichs scheme as
where the first term is the Euler scheme while the second term adds artificial diffusivity 
to the original equation:
to the original equation:
3.3)
Solution
x1 = 0;
xN = 1;
dx = 0.05;
dt = 0.02;
c = 1;
u_t0 = @(x) exp( -(x -0.5).^2 / (2*0.2^2) );
u_ex = @(x,t) exp( -(x-c*t-0.5).^2 / (2*0.2^2) );
N = (xN-x1)/dx + 1;
s = c*dt/2/dx ;
a = (1/2 + s) * ones(N,1);
b = (1/2 - s) * ones(N,1);
A = spdiags([a b], [-1 1], N,N );
A(1,:) = 0;
A(1,1) = 1;
A(N,N-1) = c*dt/dx;
A(N,N) = 1-c*dt/dx;
x = (x1:dx:xN)';
u = NaN(N,11);
u(:,1) = u_t0(x);
t = 0;
clf;
subplot(1,2,1); hold on; title('Lax-Friedrichs');
xlabel('x'); ylabel('u');
subplot(1,2,2); hold on; title('Analytical');
xlabel('x'); ylabel('u');
for n=1:10
t = t+dt;
u(1,n) = u_ex(0,t);
u(:,n+1) = A * u(:,n);
if mod(n,3)==0
subplot(1,2,1); plot(x,u(:,n+1),'DisplayName',sprintf('t=%g',t));
subplot(1,2,2); plot(x,u_ex(x,t),'DisplayName',sprintf('t=%g',t));
end
end
legend('show','Location','southeast')
Although the Lax-Friedrichs scheme is conditionally stable, it is not wise to use it because it ignores the physical properties of the original equation and just stabilizes the discretization with brute force (AKA artificial diffusivity). If we, on the other hand, want to gain some insight on what is the actual source of the numerical instability, we will analyze the PDE with the method of characteristics.
Recall that characteristic is a curve along which we can describe the evolution of u with total derivatives instead of partial derivatives. Let us choose the material derivative
where 
is the "slope" of the characteristic 
. Comparing (7) with (5) we directly obtain
is the "slope" of the characteristic . Comparing (7) with (5) we directly obtain
which after integration yields a single characteristic
The temporal evolution of u along 
is described by
is described by
Thus, the solution is steady in a co-moving reference frame 
. The inverse transformation to the steady reference frame is simply
. The inverse transformation to the steady reference frame is simply
for every 
. We can therefore compute the exact solution of (5) by replacing x with 
in the initial condition.
. We can therefore compute the exact solution of (5) by replacing x with in the initial condition.We just saw that (5) is a hyperbolic equation with single characteristic (the number of characteristics equals the number of first-order equations), so information only travels in one direction in space and time. This means that we should also reflect this fact in our spatial discretization, e.g. by upwinding. Below you can see a comparison of Lax-Friedrichs and explicit Euler with 1st and 2nd order upwind scheme. You can see that the first order upwind scheme is in this particular example more accurate than the second order Lax-Friedrichs. The second order upwind scheme apparently becomes unstable, but in can be stabilized by applying Crank-Nicolson scheme in time.
Image Credit: Vlad Giurgiu, TU Wien
3.4)
dx = 0.25;
dy = 0.2;
rhs = @(x,y) - pi^2 * sin(pi*x) .* cos(pi*y);
% 2D mesh
x = 0: dx:1;
y = 1:-dy:0;
[x,y] = meshgrid(x,y)
[I,J] = size(x);
N = I*J ;
% Discretization
e = ones(N,1);
Dx2 = spdiags( [ 1*e/dx^2 -2*e/dx^2 1*e/dx^2 ] , [ -I 0 I ] , N,N );
Dy2 = spdiags( [ 1*e/dy^2 -2*e/dy^2 1*e/dy^2 ] , [ -1 0 1 ] , N,N );
A = Dx2 + Dy2;
f = rhs(x(:),y(:));
% Boundary conditions
left = find( x(:)==0 );
right = find( x(:)==1 );
bottom = find( y(:)==0 );
top = find( y(:)==1 );
Id = speye(N);
A(left ,:) = Id(left ,:); f(left ) = 0;
A(right ,:) = Id(right ,:); f(right ) = 0;
A(bottom,:) = Id(bottom,:); f(bottom) = 0.5*sin(pi*x(bottom));
A(top ,:) = Id(top ,:); f(top ) = -0.5*sin(pi*x(top ));
clf; spy(A)
% Computation
u = reshape( A\f , I,J );
subplot(1,2,1); surf(x,y,u);
subplot(1,2,2); fsurf(@(x,y)0.5*sin(pi*x).*cos(pi*y),[0 1 0 1]);
