Various methods for the Heat Equation

21.2. Various methods for the Heat Equation#

Solve the heat equation from the previous notebook with various methods. Try also the initial condition

\[ u_0 = \chi_{[0.4,0.6]^2}, \]

a characteristic function (see geometry below).

We state the methods for the ODE

\[ y^\prime(t) = f(t, y(t)), \]

which is in our case

\[ y^\prime(t) = - M^{-1} A y(t) \]

21.2.1. Forward Euler method#

\[ \frac{y_{n+1}-y_n}{\tau} = f(y_n) \]

This method is conditionally stable. Find (by experiments) timesteps \(\tau\) such that it is working. How does \(\tau\) depend on the mesh-size and the polynomial order?

The CFL (Courant-Friedrichs-Lewy) condition gives a condition onto the timestep depending on the spectral radius \(\rho (M^{-1} A)\)

21.2.2. Trapezoidal method = Crank Nicolson method#

\[ \frac{y_{n+1} - y_n}{\tau} = \tfrac{1}{2} \left( f(y_n) + f(y_{n+1}) \right) \]

  • You should observe \(2^{nd}\)-order convergence in \(\tau\)

  • Look at \(u(x, 0.5)\) for small time, with discontinuous initial condition. You should observe the so called Crank-Nicolson bumps.

21.2.3. General Runge Kutta methods#

See very short Runge-Kutta methods

or Lecture Runge-Kutta

21.2.4. DIRK methods#

Diagonally implicit Runge-Kutta (DIRK) methods are a sub-class of RK methods, for which the \(a\)-matrix of the Butcher tableau is a lower triangular matrix:

\[\begin{split} \begin{array}{c|cccc} c_1 & a_{1,1} & \dots & 0 \\ \vdots & \vdots & \ddots & \vdots \\ c_{s} & a_{s,1} & \dots & a_{s,s} \\ \hline & b_1 & \dots & b_{s} \end{array} \end{split}\]

An example is the third-order two-stage method

\[\begin{split} \begin{array}{c|cccc} c_1 & c_1 & 0 \\ 1-c_1 & 1-2 c_1 & c_1 \\ \hline & 1/2 & 1/2 \end{array} \qquad \text{with} \qquad c_1 = \frac{3-\sqrt{3}}{6} \end{split}\]

The Runge-Kutta formulas work like that: Determine slopes

\[ k^i = f(y_n + \tau \sum_{l=1}^s a_{il} k^l) \quad \forall \; i \in \{ 1, \ldots , s \} \]

and then update for the next time-step:

\[ u_{n+1} = u_n + \tau \sum_l b_l k^l \]

For our right hand side \(f(u) = - M^{-1} A u\), the Runge-Kutta formulas read as

\[\begin{split} \begin{array}{rcl} k^1 & = & -M^{-1} A (u_n + \tau a_{11} k^1) \\ k^2 & = & -M^{-1} A (u_n + \tau a_{21} k^1 + \tau a_{22} k^2) \end{array} \end{split}\]

To determine \(k^1\) and \(k^2\), we first solve

\[ (M + \tau a_{11} A) k^1 = - A u_n, \]

and then

\[ (M + \tau a_{22} A) k^2 = - A u_n - \tau a_{21} A k^1 \]

Finally, we upate

\[ u_{n+1} = u_n + \frac{\tau}{2} ( k^1 + k^2 ) \]

21.2.5. The stability function#

  • By using an eigenvector-basis (of the evp \(A z = \lambda Mz\)) we can decompose the ODE \(u^\prime = - M^{-1} A u\) to scalar ODEs for \(u_i^\prime = \lambda_i u_i\), with \(\lambda_i \in {\mathbb R}^-\)

  • The stability function \(g(z)\) is defined such that the numerical method applied to \(y^\prime = \lambda y\) leads to the discrete propagation $\( y_{n+1} = g(\tau \lambda) y_n \)$

Examples for stability functions are:

  • Forward Euler method:

    \[ g(z) = 1+z \]

  • Backward Euler method:

    \[ g(z) = \frac{1}{1-z} \]

  • Crank-Nicolson method:

    \[ g(z) = \frac{2+z}{2-z} \]

The method is stable for certain time-steps \(\tau\) if \(\forall \lambda \in \sigma(-M^{-1} A) : | g(\tau \lambda)| \leq 1\).

For parabolic equations, \(\sigma(-M^{-1} A) \subset {\mathbb R}^-\). A method is called A-stable if \(\forall \lambda \in {\mathbb R}^- : | g(\tau \lambda)| \leq 1\).

For large negative \(\lambda\), the solution \(u_i(t) = u_{i,0} e^{-\lambda_i t}\) is falling quickly. A method is called L-stable if it is A-stable, and

\[ \lim_{z \rightarrow \infty} g(z) = 0 \]

Which methods from above are A-stable or L-stable ?

21.2.6. Some resources#

from netgen.occ import *
from ngsolve import *
from ngsolve.webgui import Draw

recto = Rectangle(1,1).Face()
inner = MoveTo(0.4,0.4).Rectangle(0.2,0.2).Face()
outer = recto-inner
outer.faces.name="outer"
inner.faces.name="inner"
inner.faces.maxh=0.05
inner.edges.hpref=1
inner.vertices.hpref=1

shape = Glue( [inner,outer])

mesh = Mesh(OCCGeometry(shape, dim=2).GenerateMesh(maxh=0.1))
mesh.RefineHP(2)
Draw (mesh);

u0 = mesh.MaterialCF( { "inner" : 1 }, default=0)
Draw (u0, mesh, deformation=True, scale=0.3)

BaseWebGuiScene