In this chapter we consider parabolic partial differential equations that are equipped with initial and boundary conditions, i. e., \[ \begin{align} \frac{\partial u}{\partial t}(x,t) + Au(x,t) +\psi (u(x,t)) &= f(x,t) && \text{in } \Omega \times (0,T) \nonumber \\ u(x,0) &= u_0(x) && \text{in } \overline{\Omega} \\ u(x,t) = b_D (x,t) \quad \text{or} \quad \nabla u (x,t) \cdot \textbf{n}(x) &= b_N (x,t) && \text{on } \partial\Omega \times [0,T),\nonumber \end{align} \] where \( (Av)(x) := \nabla \cdot (K(x)\nabla v(x)) +c(x)\cdot\nabla v(x) \) is a linear elliptic second order differential operator, the maps \(\psi:\mathbb{R}\rightarrow \mathbb{R}\) and \(f: \Omega \times (0,T) \rightarrow \mathbb{R} \) might be non-linear. We assume that \(\psi\) is differentiable. \(u_0\) represents the initial condition at time \(t=0\). Dirichlet or Neumann boundary conditions are given by \(b_D\) or \(b_N\) respectively. If we have Dirichlet boundary conditions, we need in addition that \(u_0(x) = b_D(x,0)\) on \(\partial \Omega\).
For solving parabolic equations numerically one often uses so called semi-discrete methods, i.e.,
The term semi-discrete indicates that there must follow another discretization step to make the discretization scheme complete. The idea of the semi-discretization is to deduce a familiar problem type as an intermediate step. The method of lines in particular leads to a system of ordinary differential equations which can be solved using for example Runge-Kutta methods. Rothe’s method on the other hand results into a sequence of elliptic partial differential equations. Hence one needs to solve an elliptic PDE at each timestep.
There exist also schemes that discretise in space and time simultaneously, e.g. space-time FEM. However the purpose of this section is not to cover every existing time discretization method but to give a short abstract of the basic theory.
The sequence of discretization (first space, then time or first time, then space) does not matter and leads to the same results, as long as the spatial grids are (time-independent and) uniform and the time steps are fixed. For non-uniform grids or time-varying grids and adaptively chosen time-steps, the sequence of discretization does matter. In Kaskade, the sequence first time, then space is chosen such that adaptivity in space and time can be established. Up to a given accuracy, the time discretization yields variable time steps while the resulting elliptic problems can be solved by chosing individual spatial meshes to assure the accuracy required by the temporal discretization. We will introduce some particular semi-discrete schemes and show an implementation of two schemes applied to the heat equation in the example section in the moving source tutorial.
To discretize a parabolic equation like \((1)\) in space first, we need to consider its variational formulation. For this purpose we proceed as follows:
Applying these steps to \((1)\) will lead to the following variational formulation \[ \begin{equation} \left(\frac{d}{dt}u(t), v \right ) + a(u(t),v) +\left(\psi(t),v\right)= \left(f, v \right) + \left( \nabla u(t) \cdot n , v \right)_{\partial \Omega}¸\quad \forall t \in (0,T) \text{ } \forall v \end{equation} \] with \[ a(v,w) = \int_\Omega (K\nabla v \cdot \nabla w + c \cdot \nabla v \text{ }w) dx \quad \text{and} \quad \left(v,w\right) = \int_\Omega vw \text{ } dx \quad \text{and} \quad \left(v,w\right)_{\partial\Omega} = \int_{\partial\Omega} vw \text{ } ds. \] In the following we assume homogeneous Dirichlet or homogeneous Neumann conditions such that the boundary term in \( (2) \) vanishes.
When using the method of lines we can discretize in space using Galerkin discretization with the basis functions \( \phi_1, \dots, \phi_N\) which depend on space only. We can then formulate the semi-discrete problem of \((2)\) in matrix form:
\[ \text{Find } u_h (t)= \sum_{i=1}^N u_i(t) \phi_i \text{ such that} \]
\[ \begin{align} B\dot U +A U + G(U) &=F\quad \forall t\in (0,T) \\ U(0) &= U_0, \nonumber \end{align} \] with \[ \begin{align*} B &= \left( \text{ }(\phi_j,\phi_i) \text{ }\right)^N_{i,j=1} && \text{mass matrix} \\ A &= ( \text{ }a(\phi_j,\phi_i) \text{ })^N_{i,j=1} && \text{stiffness matrix} \\ G(U) &= (\text{ }(\psi(\textstyle\sum_{j=1}^N u_i(t) \phi_j, \phi_i) \text{ })^N_{i=1} && \text{non-linear part} \\ F &= (\text{ }(f(t),\phi_i) \text{ })^N_{i=1} && \text{right-hand side} \\ U &= (\text{ } u_i (t)\text{ } )^N_{i=1} && \text{unknown component vector} \\ U_0 &= ( \text{ }u_i (0) \text{ })^N_{i=1} && \text{initial component vector}. \end{align*} \]
Note that the coefficient vector \(U\) is time dependent such that \((3)\) actually defines a system of ordinary differential equations. There are many well known methods for solving ODEs numerically which can be used to make the semi-discrete problem \((3)\) fully discrete.
JU: our tutorial (in general, i think the whole concept of using time stepping schemes in Kaskade) is based on Rothe’s method In the example section you will find implementations of two time integration schemes, namely
The following two subchapters give further descriptions of the theory they are based on.
Another way of solving instationary problems like \( (1) \) is the Rothe method. Instead of discretizing first in space, the Rothe method discretizes first in time.
A common approach is the so called One-Step \(\Theta\)-Scheme for discretizing the temporal derivative. Given the solution \(u_n\) at time \(t_n\) we can approximate the solution \(u_{n+1}\) at time \(t_{n+1}\) by solving the equation \[ \begin{align*} \frac{u_{n+1}-u_n}{\tau}+\Theta Au_{n+1}(x) +(1-\Theta) Au_{n}(x) + \Theta \psi(u_{n+1}(x)) +(1-\Theta)\psi(u_n(x) ) &=\Theta f_{n+1}(x) +(1-\Theta )f_n(x) && \text{in } \Omega \\ u_{n+1}(x) = b_D (x,t_{n+1}) \quad \text{or} \quad \nabla u_{n+1} (x) \cdot \textbf{n}(x) &= b_N (x,t_{n+1}) && \text{on } \partial\Omega \end{align*} \] with initial conditions \( u_0 (x) = u(x,0) \) in \( \overline{\Omega} \) and stepsize \(\tau = t_{n+1} -t_n\). With this scheme one has to solve a boundary value problem in each timestep which can be discretized in space using Galerkin discretization to obtain a fully discrete scheme.
For \(\Theta = 0\), one obtains the forward Euler scheme, for \(\Theta = 0.5\) the Crank–Nicolson scheme (midpoint rule), and for \(\Theta = 1\) the backward Euler scheme. Note that the explicit Euler is not stable for every stepsize \(\tau\) and that implicit Euler and Crank-Nicolson can become quite expensive if we try to solve more complicated problems, e.g. non-linear ones.
In the literature one also finds that the One-Step \(\Theta\)-Scheme is directly applied onto the variational formulation \( (2)\).
Let us consider the heat equation on a bounded domain \(\Omega\subset \R^n \) with initial condition \(u_0\) and homogeneous Dirichlet boundary conditions \[ \begin{align*} \frac{\partial u}{\partial t}(x,t)-\kappa\Delta u(x,t) &=f(x,t) && \text{in } \Omega \times (0,T) \\ u(0,x) &= u_0(x) && x \in \overline{\Omega} \\ u(t,x) &= 0 && (t,x)\in (0,T)\times \partial \Omega. \end{align*} \] such that the data \(f, u_0\) is continuous and that the consistency condition \( u_0(x) = 0 \text{ } \forall x \in \partial \Omega\) of the initial data on the boundary is satisfied.
The data \(f\) on the right-hand side can be physically interpreted as sinks and sources of heat along the domain \(\Omega\). An important observation is that if there are no sources of heat in the domain, i.e. \(f \leq 0 \) in \((0,T)\times \Omega) \), the temperature \(u\) in \( \overline{(0,T)\times \Omega} \) cannot exceed the highest temperature on the boundary. This can be proven mathematically and is known as the Maximum Principle. It is an important tool for proving uniqueness of the initial boundary value problem.
If we consider a very simple one-dimensional example, where \(\Omega=(0,1)\) and we are in the homogeneous case \(f=0\) with \(u_0 \) continuous, piecewise differentiable and \(g=0\), the unique solution is given by \[ u(t,x) =\sum_{k=0}^\infty a_k e^{-(k\pi)^2t}\text{sin}(k\pi x) \] with the Fourier coefficients \(a_k=2\int_0^1 u_0(\xi)\sin (k\pi\xi)d\xi\).
We will now apply the time integration as described above, starting with the method of lines.
We first derive the semi-discrete variational formulation (see equation \( (2) \) ) of the heat equation. Applying the steps mentioned before leads to \[ \begin{align*} \left (\frac{d}{dt} u(t),v\right) + a(u(t),v) &= (f(t) ,v) \quad \forall t\in (0,T) \forall v \qquad \text{with} \quad a(v,w) = \int_\Omega \kappa \nabla v \cdot \nabla w \text{ }dx \quad \text{and}\quad (v,w)=\int_{\Omega} vw \text{ } dx. \end{align*} \] In analogy to the Galerkin method for elliptic problems we choose a finite dimensional space \(S_h\) and approximate the continuous solution \(u(t)\) by a \(u_h(t)\in S_h\) for all \(t \in [0,T]\).
For this one dimensional example we choose the piecewise linear (P1) finite elements on \(\Omega=[0,1]\) and an equidistant distribution of \(N+2\) nodes \(x_i=ih, \text{ }i=0,\dots, N+1 \) with distance \(h=(N+1)^{-1} \) for some \(N\in\mathbb{N}\).
The basis functions \[ \phi_i(x)=\begin{cases} (x-x_{i-1})/h, & x\in [x_{i-1},x_i] \\ 1-(x-x_{i-1})/h, & x \in [x_i,x_{i+1}]\\ 0, & \text{else.}\end{cases} \quad i=1,\dots, N \] are also known as hat functions. With this basis \( \phi_1,\dots \phi_N\) of \(S_h\) and representating the approximate solution as \( u_h(t)=\sum_{k=1}^N u_k(t)\phi_k \) we can write the semi-discrete problem in matrix form \[ \begin{aligned} B\dot U +A U &=F\quad t\in (0,T) \\ U(0) &= U_0, \end{aligned} \] with \[ \text{mass matrix} \quad B = (\text{ }(\phi_j,\phi_i)\text{ })_{i,j=1}^N = \frac{h}{6}\left( \begin{array}{rrrr} 4 & 1 \\ 1 & 4 & 1 \\ & & \ddots \\ & & 1 & 4 & 1 \\ & & & 1 & 4 \\ \end{array}\right) \qquad \text{and }\qquad \text{stiffness matrix} \quad A = (\text{ }a(\phi_j,\phi_i)\text{ })_{i,j=1}^N = \frac{1}{h}\left( \begin{array}{rrrr} 2 & -1 \\ -1 & 2 & -1 \\ & & \ddots \\ & & -1 & 2 & -1 \\ & & & -1 & 2 \\ \end{array}\right). \] with \[ \begin{aligned} B &= ((\phi_i,\phi_j))_{i,j=1}^N && \text{mass matrix} \\ A &= (a(\phi_i,\phi_j))_{i,j=1}^N && \text{stiffness matrix} \\ F &= ((f(t),\phi_i))_{i=1}^N && \text{right-hand side} \\ U &= ( \mu_i(t) )_{i=1}^N && \text{unknown component vector} \\ U_0 &= (\mu_i(0) )_{i=1}^N && \text{initial component vector} \end{aligned} \]
In practice the matrix \(B/h\) becomes the identity matrix by the so called “lumping”: Basically one simply sums up the row values and writes the sum to the diagonal. This is alike using the trapezoid rule as quadrature rule for the calculation of the values of \(B\). Hereby one saves the calculation of \(M^{-1}\), and astonishingly, this simplification leads to the same approximation quality. This is reasonable as the mass matrix \(B\) leads to a correlation of the time derivatives in different points — a phenomenon which occurs due to the discretization and does not even exist in the continuous problem. In this sense, lumping can be regarded as a rectification of the discrepancy between the discrete and continuous problem.
The matrix \(A/h^2\) has the eigenvalues \( \lambda_i \approx (i\frac{\pi}{2})^2 \) for small \(h\). That means a finer space dicretization leads to a bigger range of eigenvalues. An explicit time discritization would now require a very small step-size restriction. Hence one has to use implicit methods for such problems.
The one-step \(\Theta\)-scheme for the heat equation looks as follows \[ \frac{u_{n+1}-u_n}{\tau}-\Theta \kappa \Delta u_{n+1} - (1-\Theta)\kappa\Delta u_n = \Theta f_{n+1} +(1-\Theta )f_n. \]
To be included into text and then shifted to “references” section
Deuflhard, Weiser: adaptive solutions for PDEs
John: Partial differential euqations
Kornhuber,Schuette: lecture notes Numerics III
P. Knabner, L.Angermann: Numerik partieller Differentialgleichungen