using LinearAlgebra, Plots, LaTeXStrings38 A6: Transport Equation
Assignment 6: Transport Equation and the Courant–Friedrichs–Lewy (CFL) condition
- Complete the following and submit to Canvas before April 29, 23:59
- Late work will recieve 0%,
- Each assignment is worth the same,
- Please get in contact with the TA/instructor in plenty of time if you need help,
- Before submitting your work, make sure to check everything runs as expected. Click Kernel > Restart Kernel and Run All Cells.
- Feel free to add more cells to experiment or test your answers,
- I encourage you to discuss the course material and assignment questions with your classmates. However, unless otherwise explicitly stated on the assignment, you must complete and write up your solutions on your own,
- The use of GenAI is prohibited as outlined in the course syllabus. If I suspect you of cheating, you may be asked to complete a written or oral exam on the content of this assignment.
Name:
Approx. time spent on this assignment:
In lectures, we have looked at
\begin{align} \tag{IVP} &u'(t) = f(t, u(t)) \nonumber\\ \tag{BVP} &-u'' + p u' + qu = f\nonumber\\ \tag{heat} &u_t - \alpha^2 \nabla u = f \end{align}
and various methods to approximate the solution (IVP: Euler, one-step, Runge-Kutta, multistep; BVP: finite differences (including with FFT), finite elements; Heat equation: finite differences).
In this assignment, you’ll analyse simple finite difference approximations to the linear transport equation:
\begin{align} &u_t + c u_x = 0 \nonumber\\ &u(x,0) = u_0(x) \tag{transport} \end{align}
with constant velocity c \in \mathbb R. In lectures, we saw that the solution is given by u(x,t) = u_0(x-ct). While in this case the exact solution is known, it is still useful to analyse methods for approximating this solution numerically. These methods apply more generally to more complicated hyperbolic equations.
Consider the spatial domain [a,b] and temporal domain [0,T] with meshes given as follows: x_i := a + \frac{b-a}{m} i for i = 0,\dots,m and t_j := j\frac{T}{n} for j = 0,\dots,n. Again, we define h := \frac{b-a}{m} and \tau := \frac{T}{n}. We aim to approximate u(x_i, t_j) with u_{i,j}.
Exercise 1. Suppose one approximates the linear transport equation with simple forward differences for both the space and time derivatives. Show that this approximation gives
\begin{align} u_{i,j+1} = (1+\sigma) u_{i,j} - \sigma u_{i+1,j} \tag{transport$_{\tau, h}$} \end{align}
where \sigma:= \frac{c \tau}{h}.
You answer here…
Exercise 2. We consider the linear transport equation with periodic boundary conditions: u(a,t) = u(b,t) and u_x(a,t) = u_x(b,t). Complete the following code sample to implement the finite difference scheme from exercise 1:
function fdm_transport(m, n, c, u0; a=0., b=1., T=5)
h, τ = (b-a)/m, T/n
x, t = collect(a:h:b), collect(0:τ:T)
u = zeros(m+1, n+1)
# Initial condition
u[:, 1] = u0.(x)
# CLF
σ = c * τ / h
println("σ = $(round(σ, digits=2))")
for j ∈ 1:n
for i ∈ 1:m
# periodic neighbour of m is 1:
i′ = (i%m) + 1
# finite difference rule here:
u[i,j+1] = # your code here
end
# periodic boundary
u[m+1,j+1] = # your code here
end
return x, t, u
end
a, b, T = -10, 10, 20
m, n, c = 100, 100, -1.
# initial condition (w/ PBCs)
g(x) = sin(x) * exp(-x^2)
u0(x) = g(a + mod(x - a, b-a))
# exact solution
u_e(x, t) = u0(x - c*t)
x,t,u = fdm_transport(m, n, c, u0; a, b, T)
anim = @animate for j ∈ 1:length(t)
scatter(x, u[:, j];
title = L"Transport equation $u_t + c u_{x} = 0$",
label = "t = $(t[j])", legend=:topright,
ylims = [-1/2,1/2], m=2)
plot!( x->u_e(x,t[j]);
label="exact solution", linestyle=:dot)
end
gif(anim, "pics/transport-1.gif", fps=10)σ = -1.0
[ Info: Saved animation to c:\Users\math5\Math 5485\Math5486\pics\transport-1.gif
Exercise 3. Experiment with the above code. What happens when
- (i) \sigma < -1,
- (ii) \sigma = -1,
- (iii) -1 < \sigma < 0,
- (iv) \sigma = 0, and
- (v) \sigma>0.
Either plot the numerical solution in each case, or write down the parameters you used.
You answer here…
The domain of dependence of the solution at (x,t) is the set of all y for which changing u_0(y) can affect the solution u at (x,t).
Exercise 4. Write down the domain of dependence of the solution to the linear transport equation at (x,t).
You answer here…
The numerical domain of dependence at (x_i,t_j) is the set of all \{x_k\} such that the initial data u_{k,0} can affect u_{i,j}.
Exercise 5. Write down the numerical domain of dependence of \text{(transport}_{h,\tau}\text{)}.
You answer here…
The Courant–Friedrichs–Lewy (CFL) condition states that, for a numerical method to converge, the exact domain of dependence must be contained in the numerical domain of dependence as h,\tau \to 0.
Exercise 6. Under what condition on \sigma does our numerical method for the linear transport equation satisfy the CLF condition.
You answer here…
Remark: If the CLF condition is violated, then one can change some of the initial data that is both outside the domain of numerical dependence (so does not affect the numerical solution) and inside the domain of exact dependence (so does affect the exact solution). As a result, the numerical scheme cannot approximate the true solution. The CLF condition is necessary but not sufficient for a numerical scheme to converge (as we will now see):
Consider the following centered difference approximation to the linear transport equation:
\begin{align} u_{i,j+1} = -\frac12 \sigma u_{i+1,j} + u_{ij} + \frac{1}2\sigma u_{i-1,j} \end{align}
where \sigma = \frac{c\tau}{h}.
Exercise 7. Show that, for |\sigma| \leq 1, this scheme satisfies the CLF condition.
You answer here…
Exercise 8. Show that this scheme is unstable for all \sigma\not=0 via a von Neumann stability analysis.
You answer here…
Here, we implement this method (and we can see that it is unstable for our choice of parameters):
Code
function fdm_transport_unstable(m, n, c, u0; a=0., b=1., T=5)
h, τ = (b-a)/m, T/n
x, t = collect(a:h:b), collect(0:τ:T)
u = zeros(m+1, n+1)
u[:, 1] = u0.(x)
# CLF
σ = c * τ / h
println("σ = $(round(σ, digits=2))")
for j ∈ 1:n
for i ∈ 1:m
i₋ = (i==1) ? m : i-1
i₊ = (i==m) ? 1 : i+1
u[i,j+1] = (1/2)*σ*u[i₋,j] + u[i,j] + (-1/2)*σ*u[i₊,j]
end
u[m+1,j+1] = u[1,j+1]
end
return x, t, u
end
a, b, T = -10, 10, 20
m, n, c = 100, 100, .5
# initial condition (w/ PBCs)
g(x) = sin(x) * exp(-x^2)
u0(x) = g(a + mod(x - a, b-a))
# exact solution
u_e(x, t) = u0(x - c*t)
x,t,u = fdm_transport_unstable(m, n, c, u0; a, b, T)
anim = @animate for j ∈ 1:length(t)
plot(x, u[:, j];
title = L"Transport equation $u_t + c u_{x} = 0$ (unstable)",
label = "t = $(t[j])", legend=:topright,
ylims = [-1/2,1/2], m=2)
plot!( x->u_e(x,t[j]);
label="exact solution", linestyle=:dot)
end
gif(anim, "pics/transport-2.gif", fps=10)σ = 0.5
[ Info: Saved animation to c:\Users\math5\Math 5485\Math5486\pics\transport-2.gif