Heat Equation

21.1. Heat Equation#

The heat equation is a typical example of parabolic equations. We search for the temperature $u : \Omega \times[0,T] \rightarrow {\mathbb R}$ satisfying

\[ \frac{\partial u}{\partial t}(x,t) - \Delta u(x,t) = f(x,t) \qquad \forall \, x \in \Omega, \forall t \in (0,T), \]

where $f$ is a given heat source density. Here, $\Omega \subset R^d$ is the spatial domain, and $[0,T]$ is the time interval. In addition, we need initial conditions

\[ u(x,0) = u_0(x) \qquad \forall \, x \in \Omega, \]

and boundary conditions at the spatial boundary. For example Dirichlet boundary conditions read as

\[ u(x,t) = u_D(x,t) \qquad \forall \, x \in \partial \Omega, \; \forall t \in(0,T] \]

21.1.1. Variational formulation in space#

Similar as for elliptic equations we pass over to the weak formulations. We do this only for the spatial variables $x$, for all $t \in (0,T)$. The test functions $v$ depend only on $x$. We multiply by $v$, integration over $\Omega$, and perform integration by parts to obtain:

\[ \int_\Omega \frac{\partial u}{\partial t}(x,t) v(x) \, dx + \int \nabla u(x,t) \nabla v(x) \, dx = \int_\Omega f(x,t) v(x) \, dx \qquad \forall \, v \in H_0^1(\Omega) \; \forall t \, \in (0,T) \]

The initial condition in $L_2$ is set for $t = 0$: $$ u(\cdot, 0) = u_0 \qquad \text{in} \; L_2(\Omega) $$

We consider the unknown function $u$ as a function with values in the Hilbert space $H^1$, e.g. $u : [0,T] \rightarrow H^1(\Omega)$. This allows to express the equation as Hilbert-space valued ordinary differential equation: $$ \frac{d u }{dt }(t) + A u(t) = f(t) $$

21.1.2. Galerkin method in space#

As a next step we apply the finite element method in space. Let $V_h$ be a finite element sub-space of $H^1(\Omega)$, with the basis functions $\{ p_1(x), \ldots p_N(x) \}$. We approximate the time-dependent function by

\[ u_h : [0,T] \rightarrow V_h, \]

which can we expressed by means of the basis as

\[ u_h(t) = \sum_{i=1}^N u_i(t) p_i(x) \]

Inserting this expansion into the variational formulation, and using spatial test-functions $v = p_j$ we obtain the semi-discrete equation:

\[ \int_\Omega \frac{\partial}{\partial t} \Big\{ \sum_{i=1}^N u_i(t) p_i(x) \Big\} p_j(x) dx + \int_\Omega \nabla \Big\{ \sum_{i=1}^N u_i(t) p_i(x) \Big\} \; \nabla p_j(x) = \int_\Omega f(x,t) p_j(x) \, dx \qquad \forall \, j = 1 \ldots N, \forall \, t \in (0,T) \]

The initial condition for $u_h$ is obtained by some kind of projection, typically we want the $L_2$ projection

\[ \int u_h(x,0) p_j(x) \, dx = \int u_0(x) p_j(x) \, dx \qquad \forall \, j \in \{ 1, \ldots N \} \]

By means of the mass matrix $M \in R^{N \times N}$ $$ M = \left\{ \int_\Omega p_i p_j \, dx \right\}_{i,j = 1, \ldots N} $$ and **stiffness matrix** $A \in R^{N \times N}$ A = \left\{ \int_\Omega \nabla p_i \nabla p_j \, dx \right\}_{i,j = 1, \ldots N} $$ we obtain the ordinary differential equation (ODE)

\[ M \dot u(t) + A u(t) = f(t) \qquad \forall \, t \]

including the initial condition

\[ u(0) = u_0 \]

for the unknown coefficient function $u(\cdot) = (u_1(\cdot), \ldots u_N(\cdot)) \in C^1([0,T], {\mathbb R}^N)$.

21.1.3. Backward Euler time-stepping#

We approximate the function $u$ at time-steps $t_j = j \tau$ by $u^j$, with the time-step $\tau$ and $j \in \{0, \ldots m\}$. By replacing the time-derivative $\dot u$ by the backward difference quotient we obtain

\[ M \frac{u^{j}-u^{j-1}}{\tau} + A u^j = f^j \]

The initial value $u^0 := u_0$ is given by the initial condition. If the old time-step solution $u^{j-1}$ is known, we can compute the next step by the linear system of equations

\[ (M + \tau A) u^j = \tau f + M u^{j-1}. \]

Often the new step is computed in incremental form:

\[\begin{split} \begin{array}{rcl} (M + \tau A) w^j &=& f - A u^{j-1} \\ u^j & = & u^{j-1} + \tau w^j \end{array} \end{split}\]

from ngsolve import *
from ngsolve.webgui import Draw
from netgen.occ import unit_square
from time import sleep

mesh = Mesh(unit_square.GenerateMesh(maxh=0.05))

setting up matrices:

tend = 0.1
u0 = exp(-100*( (x-0.5)**2 + (y-0.5)**2))

fes = H1(mesh, order=3)
u,v = fes.TnT()

m = BilinearForm(u*v*dx).Assemble()
a = BilinearForm(grad(u)*grad(v)*dx).Assemble()
f = LinearForm(fes).Assemble()

tau = 0.005
mstarinv = (m.mat + tau*a.mat).CreateSparseMatrix().Inverse()

gfu = GridFunction(fes)
gfu.Set(u0)

scene = Draw (gfu, deformation=True, euler_angles=[-60,-7,-15])
sleep(3)

for j in range(int(tend/tau)):
    res = f.vec - a.mat * gfu.vec
    w = mstarinv * res
    gfu.vec.data += tau*w
    scene.Redraw()
    sleep(0.2)

21.1.3.1. Linear functionals to evaluate total energy, and point values:#

ftot = LinearForm(1*v*dx).Assemble()
print ("total heat energy: ", InnerProduct(ftot.vec, gfu.vec))

print ("gfu(0.5, 0.5) =", gfu(0.5, 0.5))

fpoint = LinearForm(fes)
fpoint += v(0.5, 0.5)
fpoint.Assemble()
print ("point value: ", InnerProduct(fpoint.vec, gfu.vec))

total heat energy:  0.03141604953852029
gfu(0.5, 0.5) = 0.034651659201342375
point value:  0.034651659201342375

21.1.3.2. Exercises:#

plot total heat energy as a function of time. Try also with (partial) homogeneous Dirichlet boundary conditions.
plot temperature $u(0.6,0.5)$ as a function of time. Plot for various time-steps
plot $u(x,0.5)$ for various times