52. Heat Equation

The heat equation is a typical example of parabolic equations. We search for the temperature \(u : \Omega \times[0,T] \rightarrow 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] \]

52.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) \]

52.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], R^N)\).

52.3. Implicit 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{eqnarray*} (M + \tau A) w^j &=& f - A u^{j-1} \\ u^j & = & u^{j-1} + \tau w^j \end{eqnarray*}\]

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:

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

fes = H1(mesh, order=3)
u,v = fes.TnT()
mform = u*v*dx
aform = grad(u)*grad(v)*dx

m = BilinearForm(mform).Assemble()
a = BilinearForm(aform).Assemble()
mstar = BilinearForm(mform+tau*aform).Assemble()
mstarinv = mstar.mat.Inverse()
f = LinearForm(fes).Assemble()

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

scene = Draw (gfu, deformation=True)
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)