The heat equation

Authors: Anders Logg and Hans Petter Langtangen

Minor modifications by: Jørgen S. Dokken

As a first extension of the Poisson problem from the previous chapter, we consider the time-dependent heat equation, or the time-dependent diffusion equation. This is the natural extension of the Poisson equation describing the stationary distribution of heat in a body to a time-dependent problem. We will see that by discretizing time into small time intervals and applying standard time-stepping methods, we can solve the heat equation by solving a sequence of variational problems, much like the one we encountered for the Poisson equation.

The PDE problem

The model problem for the time-dependent PDE reads

(17)\[\begin{align} \frac{\partial u}{\partial t}&=\nabla^2 u + f && \text{in } \Omega \times (0, T],\\ u &= u_D && \text{n } \partial\Omega \times (0,T],\\ u &= u_0 && \text{at } t=0. \end{align}\]

Here \(u\) varies with space and time, e.g. \(u=u(x,y,t)\) if the spatial domain \(\Omega\) is two-dimensional. The source function \(f\) and the boundary values \(u_D\) may also vary with space and time. The initial condition \(u_0\) is a function of space only.

The variational formulation

A straightforward approach to solving time-dependent PDEs by the finite element method is to first discretize the time derivative by a finite difference approximation, which yields a sequence of stationary problems, and then turn each stationary problem into a variational formulation. We will let the superscript \(n\) denote a quantity at time \(t_n\), where \(n\) is an integer counting time levels. For example, \(u^n\) means \(u\) at time level \(n\). The first step of a finite difference discretization in time consists of sampling the PDE at some time level, for instance \(t_{n+1}\)

(18)\[\begin{align} \left(\frac{\partial u }{\partial t}\right)^{n+1}= \nabla^2 u^{n+1}+ f^{n+1}. \end{align}\]

The time-derivative can be approximated by a difference quotient. For simplicity and stability reasons, we choose a simple backward difference:

(19)\[\begin{align} \left(\frac{\partial u }{\partial t}\right)^{n+1}\approx \frac{u^{n+1}-u^n}{\Delta t}, \end{align}\]

where \(\Delta t\) is the time discretization parameter. Inserting the latter expression into our equation at time step \(n+1\) yields

(20)\[\begin{align} \frac{u^{n+1}-u^n}{\Delta t}= \nabla^2 u^{n+1}+ f^{n+1}. \end{align}\]

This is our time-discrete version of the heat equation. It is called a backward Euler or a implicit Euler discretization.

We reorder the equation such that the left-hand side contains the terms with only the unknown \(u^{n+1}\) and right-hand side contains only computed terms. The resulting equation is a sequence of stationary problems for \(u^{n+1}\), assuming \(u^{n}\) is known from the previous time step:

(21)\[\begin{align} u^0&=u_0 &&\\ u^{n+1}-\Delta t \nabla^2 u^{n+1}&= u^{n} + \Delta t f^{n+1}, && n = 0,1,2,\dots \end{align}\]

Given \(u_0\), we can solve for \(u^0, u^1, u^2\) and so on.

We then in turn use the finite element method. This means that we have to turn the equation into its weak formulation. We multiply by the test-function of \(v\in \hat{V}\) and integrate second-order derivatives by parts. we now introduce the symbol \(u\) for \(u^{n+1}\) and we write the resulting weak formulation as

(22)\[\begin{align} a(u,v)&=L_{n+1}(v), \end{align}\]

where

(23)\[\begin{align} a(u,v)&=\int_{\Omega}(uv + \Delta t \nabla u \cdot \nabla v )~\mathrm{d} x\\ L_{n+1}(v)&=\int_{\Omega} (u^n+\Delta t f^{n+1}) \cdot v~\mathrm{d} x. \end{align}\]

Projection or interpolation of the initial condition

In addition to the variational problem to be solved in each time step, we also need to approximate the initial condition. This equation can also be turned into a variational problem

(24)\[\begin{align} a_0(u,v)&=L_0(V), \end{align}\]

with

(25)\[\begin{align} a_0(u,v)&=\int_{\Omega}uv~\mathrm{d} x,\\ L_0(v)&=\int_{\Omega}u_0v~\mathrm{d} x. \end{align}\]

When solving this variational problem \(u^0\) becomes the \(L^2\)-projection of the given initial value \(u_0\) into the finite element space.

The alternative is to construct \(u^0\) by just interpolating the initial value \(u_0\). We covered how to use interpolation in DOLFINx in the membrane chapter.

We can use DOLFINx to either project or interpolate the initial condition. The most common choice is to use an projection, which computes an approximation to \(u_0\). However, in some applications where we want to verify the code by reproducing exact solutions, one must use interpolate. In this chapter, we will use such a problem.