# 66. 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

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

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

## 66.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:

The initial condition in \(L_2\) is set for \(t = 0\):

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:

## 66.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

which can we expressed by means of the basis as

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

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

By means of the **mass matrix** \(M \in R^{N \times N}\)

and **stiffness matrix** \(A \in R^{N \times N}\)

we obtain the ordinary differential equation (ODE)

including the initial condition

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

## 66.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

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

Often the new step is computed in incremental form:

```
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)
```