Project 5
This assignment can be worked on and submitted in groups of up to three
persons.
Goals. In this project, we will work with parabolic and hype
olic partial
differential equations, and combine elementary techniques from 2p-BVPs
with time-stepping methods from IVPs. The goal is to gain experience with
the method of lines and stability constraints of Courant-Friedrichs-Lewy
(CFL) type on the time-step ∆t.
In the first part we shall primarily work with the parabolic diffusion
equation ut = uxx. In the second part we will work with conservation laws such
as the hype
olic advection equation ut = −ux.
The objective is to gain an elementary understanding of and experience
with “evolution equations” (time dependent partial differential equations),
and the many varying properties one encounters in such problems. In
addition, we will work with different types of boundary conditions, and it is
important to implement boundary conditions co
ectly in order to get the
co
ect results.
We will only work in one space dimension, and on equidistant grids. We
will use explicit and implicit time-stepping methods, in order to study the
properties of some important elementary finite difference methods. We will
also observe some undesirable phenomena such as numerical instability and
numerical damping, and learn how to cope with them.
Part 1. The diffusion equation
The diffusion equation is
ut = uxx
u(t, 0) = u(t, 1) = 0
u(0, x) = g(x)
and its method-of-lines semi-discretization is
u̇ = T∆xu
1
on an equidistant grid on [0, 1], with ∆x = 1/(N + 1) and initial condition
g(x) sampled on the grid. Generate the grid using linspace, so that you have
N internal points on (0, 1). The matrix T∆x is the usual Toeplitz symmetric
tridiagonal matrix that approximates ∂2/∂x2 to second order accuracy on the
grid. In Project 3 you have generated this (and similar operators. Go back to
your old programs and retrieve the code segments you need in order to create
T∆x. You will need (some representation of) the operator T∆x in order to solve
equation systems involving this matrix.
Next, write your Matlab file eulerstep and make sure that it takes one
time step, and that it works in the matrix–vector case. It will probably be
necessary to make some minor modifications so that you can call eulerstep
with a command of the type
unew = eulerstep(Tdx,uold,dt)
in order to solve u̇ = T∆xu by the Euler method, um+1 = um + ∆t · T∆xum.
Task 1.1 (1p) Write a script that constructs T∆x, assuming N interior
points on the grid and solves the semidiscretization u̇ = T∆xu by the Euler
method.
1. Visualize the computational procedure by plotting the numerical solu-
tion u(m∆t, n∆x) over the t, x plane. Basically, you should apply the
following procedure. You determine your spatial grid size ∆x and time
step ∆t from the number N of internal mesh points on the x-interval
[0, 1] and the number M of time steps to go from t = 0 to t = tend,
espectively. Always run your experiments by varying N and M , not
y varying ∆x and ∆t directly. Then construct a vecto
xx = linspace(0,1,N+2)
and a vecto
tt = linspace(0,tend,M+1)
and run [T,X]=meshgrid(tt,xx) to generate the independent vari-
ables for 3D graphics in Matlab. If you are not familiar with these
commands, check Matlab’s help.
When you solve the PDE, save the solution vectors over the meshed
grid. Don’t forget to insert the boundary values! When the computa-
tion is done, make a 3D color graph using either the mesh or the surf
command. (Run help to find out how they work.)
2
2. Determine experimentally the CFL condition on the time step ∆t fo
the explicit Euler method applied to the semi-discretization. Choose
your initial condition g(x) as you like, but make sure that they are com-
patible with the boundary conditions at x = 0 and x = 1, respectively,
and make sure that g(x) is not an eigenfunction of ∂2/∂x2.
3. Try to obtain representative plots of the solution both when you fulfill
the CFL condition and when you violate it. In the latter case, don’t
overdo it; try to see what happens when you violate the CFL condition
y a small amount only, so that the numerical instability is clearly
visible.
Task 1.2 (1p) Next implement the Crank–Nicolson method. This is the
same as applying the trapezoidal rule to the equation,
un+1 = un +
∆t
2
(T∆xu
n + T∆xu
n+1).
Implement a function
unew = TRstep(Tdx,uold,dt)
Solve the same diffusion equation as in Task 1. Verify that you can now
obtain stable solutions where you previously observed instability due to a
violated CFL condition. Try to verify that you can use much larger time
steps ∆t than before, but remember to vary time step and Courant numbe
y varying N and M . Check the the results are the same no matter what
initial conditions you use.
Part 2. The advection equation
The linear advection equation is
ut + aux = 0,
with appropriate initial and boundary conditions. We will only work with
periodic boundary conditions. This means that the solution satisfies u(t, 0) =
u(t, 1) for all times t. That is, the solution “wraps around” from x = 1 to
x = 0, and vice versa.
With periodic boundary conditions, we represent the solution u on the
entire grid, including the boundaries throughout the entire computation. Let
the solution vector at time t be an N -vector and introduce an equidistant
3
grid on [0, 1], such that 0 = x1, . . . , xj = (j−1)∆x, . . . , xN = (N−1)∆x and
xN+1 = 1.
This implies that ∆x = 1/N , and that we are going to compute the
solution at the N points x1, . . . xN . The solution at the right boundary,
xN+1 = 1, is uN+1 and is identified with u1 at the left boundary, x1 = 0.
Make sure you have understood exactly what this implies (make a sketch of
the grid!), as your code will differ from what you did in the parabolic case.
Initial conditions are always needed, and have the form
u(0, x) = g(x).
Note that the initial condition must satisfy g(0) = g(1) and g′(0) = g′(1)
when periodic boundary conditions are used. Do not use a trigonometric
initial condition such as sin 3πx.
One of the simplest methods of order 2 for the advection equation is the
Lax–Wendroff scheme. It can be derived using a Taylor series expansion.
Thus,
u(t+ ∆t, x) = u(t, x) + ∆tut +
∆t2
2!
utt + O(∆t
3).
From the differential equation ut + aux we get utt = −auxt, as well as
utx = −auxx. Therefore we have utt = a2uxx provided that the solution
u is sufficiently differentiable. Inserting this into the Taylor series expansion
we get
u(t+ ∆t, x) = u(t, x)− a∆tux +
a2∆t2
2!
uxx + O(∆t
3).
This gives the Lax–Wendroff scheme,
un+1j =
aµ
2
(1 + aµ)unj−1 + (1− a2µ2)unj −
aµ
2
(1− aµ)unj+1.
Note that the method coefficients are not symmetric – instead, they change
with the flow direction as determined by the sign of a.
Task 2.1 (1p) Implement a Lax–Wendroff solver for the scalar advection
equation with periodic boundary conditions. Write a function
function unew = LaxWen(u, amu)
that takes a step of size ∆t starting from the “initial condition” u. The
parameter amu is the product a∆t/∆x. Write a script for running your Lax–
Wendroff solver. Test your solver using some pulse-like initial data, e.g.
g(x) = e−100(x−0.5)
2
.
4
Check how this pulse propagates for 5 units of time, and verify that the
method works both for positive and negative a. Pay close attention to how
you implement the periodic boundary conditions (sketch the grid and the
computational stencil before you write the program).
Produce plots of the solution for at least two different CFL numbers. In
addition, plot the L2 norm (RMS norm) of the solution vs time t ∈ [0, 5]
for a∆t/∆x = 1 as well as for a∆t/∆x = 0.9, and explain the difference.
Is the Lax-Wendroff method conserving the norm of the solution? Does the
amplitude remain constant, or is it damped? Motivate why one actually
wants to run the code at the CFL stability limit.
5
Matlab Script Codes
Prj5p1
clear all;close all;
%Part 1a: The Diffusion Equation
%
% ut= uxx
% u(t,0)= u(t,1)= 0
% u(0,x)= g(x)
%
% dx= 1/(N+1)
% u_t= Tdx*u ~ u^(m+1)= u^m + delta_t*Tdx*u^m <--- This is the Explicit Euler Form
%
% Exact solution
ga = @(x) sin(pi*x) - sin(3*pi*x);
ua = @(x,t) exp(-(pi)^2*t)*sin(pi*x) - exp(-(3*pi)^2*t)*sin(3*pi*x);
% Spatial Finite Difference Discretization
N = 30;
dx = 1/(N+1); % equidistance grid on [0,1] with N+1 subintervals
T = diag(-2*ones(N,1),0) + diag(ones(N-1,1),-1 )+ diag(ones(N-1,1),1);
Tdx = 1/(dx^2)*T;
xx= linspace(0,1,N+2); % equidistance grid on [0,1] with N+1 points
xx = xx';
% Temporal Euler Explicit
tend = 0.1;
Mg = [200, 50];
for i = 1:length(Mg)
M = Mg(i)
dt = tend/M % total M time stepping
tt= linspace(0,tend,M+1); % time instants
Z= zeros(N+2,M+1); % soution u(x_i,t_m)
Za = zeros(N+2,M+1); % exact solution
e
= zeros(M+1,1);
Z(:,1)= ga(xx); % initial condition u(x_i,0) = g(x_i)
Za(:,1) = ua(xx,tt(1));
e
(1) = norm(Za(:,1)-Z(:,1));
for m = 1:M
Z(:,m+1) = eulerstep(Tdx,Z(:,m),dt); % explicit Eule
Za(:,m+1) = ua(xx,tt(m+1));
e
(m+1) = norm(Za(:,m+1)-Z(:,m+1));
end
figure(i)
[T,X]= meshgrid(tt,xx);
subplot(3,1,1)
surf(T,X,Za) % This is the 3D plot of the exact solution
subplot(3,1,2)
plot(tt,e
,'o-') % This is the 2D plot of the e
ors between the exact solution and computed ones.
subplot(3,1,3)
surf(T,X,Z) %This is going to give us the 3D plot
CFL = dt/dx^2
nrm0 = norm(Z(:,M));
nrm1 = norm(Z(:,M+1));
if (CFL < 1/2)
title('CFL Condition are met and solutions should be stable')
elseif (CFL>1/2 & nrm1
title('CFL Condition are not met and solutions are stable')
elseif (CFL>1/2 & nrm1>nrm0)
title('CFL Condition are not met and solutions are not stable')
end
end
Prj5p1
clear all;close all;
%Part 1.2: Solving the Diffusion Equation by the Crank-Nicolson Method
%
% ut= uxx
% u(t,0)= u(t,1)= 0
% u(0,x)= g(x)
%
% Spatial Finite Difference Discretization