18 Transport and Heat Equations
Lecture 18
Note
I encourage you to play around with the juptyer notebook for this lecture - you can copy the code with the this notebook button on the side of this page.
Warning
These notes are mainly a record of what we discussed and are not a substitute for attending the lectures and reading books! If anything is unclear/wrong, let me know and I will update the notes.
Tip
This chapter is mostly based on
- Burden, Faires, and Burden (2015) Chapter 11,
- Driscoll and Braun (2018) Chapter 10,
- Lecture 16: Driscoll and Braun (2018) Galerkin method,
- Lecture 17: Chapter 14.4 of Suli and Mayers (2012): Theorem 14.6.
18.1 Recap
We have seen how to approximate the solution to general initial value problems: Find u:[0,T]\to \mathbb R such that
\begin{align} &u'(t) = f(t, u(t) ) \nonumber\\ &u(0) = u_0. \end{align}
(e.g exponential growth/decay, population models, ….) We saw Euler’s method, general one-step methods, Runge-Kutta methods, multistep methods (and their associated error estimates). We then saw linear boundary value problems of the form: find u : [a,b] \to \mathbb R such that
\begin{align} -u''(x) &+ p(x) u'(x) + q(x) u(x) = f(x) \nonumber\\ &u(a) = A, \quad u(b) = B \tag{Dirichlet} \\ \text{or} \qquad &u(a) = u(b), \quad u'(a) = u'(b) \tag{Periodic} \end{align}
(e.g. semiconductors, pollutants, drug delivery, asset prices, heat transfer in fluids, …) We saw finite differences and finite elements and some assosciated error estimates. Now, we’ll take what we have learnt so far and apply it to Partial Differential Equations (PDEs).
18.2 Some examples…
PDEs can model a variety of phenomena. Here, we list some examples to introduce some of the notation:
Laplace’s equation:
\begin{align} - \Delta u := \sum_{i=1}^d \frac{\partial^2 u}{\partial x_i^2} = 0 \end{align}
Here, \Delta is the Laplacian. Solutions to this equation are called Harmonic. Poisson’s equation is the inhomogenous version: - \Delta u = f. This models heat conduction, elastic deformation, appears in electrostatics, Newton’s law of gravitation, ….
Helmholtz’s equation:
\begin{align} - \Delta u = \lambda u \end{align}
This is an eigenvalue equation which appears in the study of acoustics, optics, electromagnetism, it’s related to time-independent Schr"odinger equation, ….
Heat (diffusion) equation:
\begin{align} u_t - \Delta u = 0 \end{align}
where u_t := \frac{\partial u}{\partial t}. Models the dissipation of some quantity as a function of time. (Related to the Black-Scholes equation in finance).
Schrödinger equation:
\begin{align} i \hbar \frac{\partial \Psi}{\partial t} = \hat{H} \Psi \end{align}
where \hat{H} is the Hamiltonian and \Psi is the wave function representing the distribution of electrons in a quantum system.
Wave equation:
\begin{align} u_{tt} - \Delta u = 0. \end{align}
So far, we have only seen linear equations (the equations are linear in u and it’s derivatives). Nonlinear equations are also of importance in many fields:
Minimal surfaces:
\begin{align} \nabla \cdot \left( \frac{\nabla u}{\sqrt{1 + |\nabla u|^2}} \right) = 0. \end{align}
This equation describes surfaces with zero mean curvature - this is the shape of soap films stretched across a wire boundary. Here, \nabla \cdot v := \sum_{i=1}^d \frac{\partial v_{i}}{\partial x_i} is the divergence of v.
Monge-Ampère equation:
\begin{align} \det D^2 u = f. \end{align}
Here, D^2 u := (u_{x_i x_j})_{i,j=1}^n denotes the Hessian. This equation arises in optimal transport (moving mass efficiently), differential geometry (prescribing the Gauss curvature of a surface), and meteorology.
Burgers’:
\begin{align} u_t + (u \cdot \nabla) u = \nu \Delta u \end{align}
(\nu=0: inviscid Burgers’)
Korteweg-de Vries (KdV):
\begin{align} u_t + 6 u u_x + u_{xxx} = 0 \end{align}
Waves in shallow water.
In this course, we will look at linear PDEs which can be categorised as Elliptic, Parabolic, or Hyperbolic: consider the general linear equation
\begin{align} A u_{\xi\xi} + B u_{\xi\zeta} + C u_{\zeta\zeta} + \mathrm{l.o.t} =0 \end{align}
where \mathrm{l.o.t} represents “lower order terms” (functions involving derivatives of lower order). Here, we are thinking of either (\xi,\zeta) = (x,t) (one spatial dimension) or (\xi, \zeta) = (x,y) (two spatial dimensions and no time dependence). To this equation, we may assign the conic section A \xi^2 + B \xi \zeta + C \zeta^2 = 0. The equation is classed as Elliptic, Parabolic, or Hyperbolic if the corresponding conic section is Elliptic, Parabolic, or Hyperbolic, respectively. That is,
\begin{align} &B^2 - 4AC < 0 &&\text{Elliptic} \\ &B^2 - 4AC = 0 &&\text{Parabolic} \\ &B^2 - 4AC > 0 &&\text{Hyperbolic} \end{align}
Roughly speaking, Elliptic equations describe steady-state systems, Parabolic equations describe diffusion (systems that evolve towards equilibrium), and Hyperbolic equations describe wave propagation.
Exercise 18.1 Show that Laplace’s equation (in 2d) is Elliptic, the wave equation is Hyperbolic, and the heat equation is Parabolic
18.3 Laplace equation (Elliptic)
See presentation on finite differences for Laplace equation in 2d.
18.4 Transport/advection equation (Hyperbolic)
The simplest possible equation of hyperbolic type is the linear transport equation with constant coefficient:
\begin{align} & u_t + b \cdot \nabla u = 0 \nonumber\\ & u(x,0) = u_0(x). \end{align}
Note
We use the notation \nabla u for the vector of first derivatives with respect to the spatial dimensions: (\nabla u)_i = u_{x_i} = \frac{\partial u}{\partial x_i} = \partial_{x_i} u = \partial_i u = \cdots. The dot is the Euclidean dot product so that b \cdot \nabla u := \sum_{i=1}^d b_i u_{x_i}.
Exercise 18.2 Show that u(x,t) := u_0(x - bt) is a solution.
Answer. You can show this directly but here we look at a more general method.
We will use the so-called method of characteristics - the idea is to reduce the PDE to a system of ODEs by trying to parametrise the solution as z(s) := u\big( x(s), t(s) \big) where x(0) = \xi and t(0) = 0. Then, z satisfies the ODE:
\begin{align} z(0) &= u( \xi, 0 ) = u_0(\xi) \nonumber\\ z'(s) &= t'(s) u_t\big( x(s), t(s) \big) + x'(s) \cdot \nabla u\big(x(s),t(s)\big) \end{align}
Here, we have applied the chain rule. Therefore, if
\begin{align} &t(0) = 0, &&t'(s) = 1 \nonumber\\ &x(0) = \xi &&x'(s) = b \end{align}
(these are known as the characteristic equations), then z'(s) = 0 and z(s) is constant (and equal to z(0) = u_0(\xi)). The solution to the characteristic equations are given by x(s) = \xi + sb and t(s) = s, and so z(s) = u(\xi + sb, s) = u_0(\xi) = z(0) and thus u(x,t) = u_0(x-tb). \square
We have also shown that for the linear transport equation, the solution is constant along the characteristic curves (\xi+sb, s):
Here, we look at what these solutions do:
[ Info: Saved animation to c:\Users\math5\Math 5485\Math5486\pics\transport-1.gif
Inhomogeneous equation:
\begin{align} & u_t + b \cdot \nabla u = f \nonumber\\ & u(x,0) = u_0(x) \end{align}
Exercise 18.3 Write down an integral form of the solution to this equation.
Transport with decay: take \mu > 0 and consider
\begin{align} &u_t + b \cdot \nabla u + \mu u = 0\nonumber\\ &u(x,0) = u_0(x) \end{align}
Exercise 18.4 Applying the exact same argument as before, we find that z'(s) + a z(s) = 0 (where z(s) = u(\xi + sb, s)). Use the integrating factor to conclude that the exact solution is
\begin{align} u(x,t) = u_0(x-bt) e^{-\mu t} \end{align}
[ Info: Saved animation to c:\Users\math5\Math 5485\Math5486\pics\transport-2.gif
18.5 Heat equation (Parabolic)
The heat equation (or the diffusion equation) is
\begin{align} u_t - \alpha^2 \Delta u = f \end{align}
This models e.g. how temperature density u distributes (diffuses) accross a domain \Omega \subset \mathbb R^d. Here, \alpha^2 is the thermal diffusivity.
Exercise 18.5 Take d = 1. Explain why the term -\alpha^2 u'' makes the solution “smooth out” (diffuses) over time.
Exercise 18.6 Consider the equation
\begin{align} &u_t - \alpha^2 u_{xx} = 0 \qquad \text{on }(0,1)\nonumber\\ \tag{IC} &u(x,0) = \sin \pi x \nonumber\\ \tag{BC} &u(0,t) = u(1,t) = 0. \end{align}
- Show that u(x,t) = e^{-\pi^2 t} \sin \pi x is a solution to this equation.
Here, we look at a finite difference implementation of this equation:
function fdm_heat( m, n, α, u0; a=0., b=1., f=0., g_a = 0., g_b = 0., T=.2, warn=true )
# space discretisation
h = (b-a)/m
x = collect(a:h:b)
# time discretisation
τ = T/n
t = collect(0:τ:T)
c = τ*(α/h)^2
if ((c > 1/2) && (warn))
@warn "c = $c, unstable!"
end
u = zeros( m+1, n+1 )
u[:,1] = u0.(x)
# boundary condition @ t=0
u[1,1] = typeof(g_a) <: Function ? g_a(t[1]) : g_a
u[m+1,1] = typeof(g_b) <: Function ? g_b(t[1]) : g_b
for j ∈ 1:n
for i ∈ 2:m
source = typeof(f) <: Function ? f(x[i],t[j]) : f
u[i,j+1] = u[i,j] + c * (u[i-1,j]-2u[i,j]+u[i+1,j]) + τ*source
end
u[1,j+1] = typeof(g_a) <: Function ? g_a(t[j+1]) : g_a
u[m+1,j+1] = typeof(g_b) <: Function ? g_b(t[j+1]) : g_b
end
return x,t,u
end
# initial condition
u0(x) = sin(π*x)
# exact solution
u_e(x,t) = exp( -(π^2) * t ) * sin( π*x )
x,t,u = fdm_heat(50, 1000, 1., u0)
anim = @animate for j ∈ 1:10:length(t)
scatter(x, u[:, j];
title = L"Heat equation $u_t - u_{xx} = 0$",
label = "t = $(t[j])",
ylims = [0,1], m=2)
plot!( x->u_e(x,t[j]);
label="exact solution", linestyle=:dot)
end
gif(anim, "pics/heat-1.gif")[ Info: Saved animation to c:\Users\math5\Math 5485\Math5486\pics\heat-1.gif
Here, we run the same code with different mesh sizes:
x,t,u = fdm_heat(60, 1000, 1., u0)
anim = @animate for j ∈ 1:5:length(t)
plot(x, u[:, j];
title = L"Heat equation $u_t - u_{xx} = 0$",
label = "t = $(t[j])", m=2)
plot!( x->u_e(x,t[j]);
label="exact solution", linestyle=:dot)
end
gif(anim, "pics/heat-2.gif", fps=3)┌ Warning: c = 0.7200000000000001, unstable! └ @ Main In[5]:12 [ Info: Saved animation to c:\Users\math5\Math 5485\Math5486\pics\heat-2.gif
