$$ \newcommand{\half}{\frac{1}{2}} \newcommand{\tp}{\thinspace .} \newcommand{\uex}{{u_{\small\mbox{e}}}} \newcommand{\Aex}{{A_{\small\mbox{e}}}} \newcommand{\x}{\boldsymbol{x}} \newcommand{\dfc}{\alpha} % diffusion coefficient \newcommand{\If}{\mathcal{I}_s} % for FEM \newcommand{\Ifb}{{I_b}} % for FEM \newcommand{\basphi}{\varphi} \newcommand{\baspsi}{\psi} \newcommand{\xno}[1]{x_{#1}} \newcommand{\dx}{\, \mathrm{d}x} \newcommand{\ds}{\, \mathrm{d}s} $$

Study guide: Time-dependent problems and variational forms

Hans Petter Langtangen [1, 2]

[1] Center for Biomedical Computing, Simula Research Laboratory [2] Department of Informatics, University of Oslo

Oct 30, 2015

Time-dependent problems
Example: diffusion problem
A Forward Euler scheme; ideas
A Forward Euler scheme; stages in the discretization
A Forward Euler scheme; weighted residual (or Galerkin) principle
A Forward Euler scheme; integration by parts
New notation for the solution at the most recent time levels
Deriving the linear systems
Structure of the linear systems
Computational algorithm
Example using sinusoidal basis functions
Approximating the initial condition
Computing the $ M $ and $ K $ matrices
Solving the equation system
Comparing P1 elements with the finite difference method; ideas
Comparing P1 elements with the finite difference method; results
Discretization in time by a Backward Euler scheme
The variational form of the time-discrete problem
Calculations with P1 elements in 1D
Dirichlet boundary conditions
Boundary function
Finite element basis functions
Modification of the linear system; the raw system
Modification of the linear system; setting Dirichlet conditions
Modification of the linear system; Backward Euler example
Analysis of the discrete equations
Handy formulas
Amplification factor for the Forward Euler method; results
Amplification factor for the Backward Euler method; results
Amplification factors for smaller time steps; Forward Euler
Amplification factors for smaller time steps; Backward Euler

Time-dependent problems

So far: used the finite element framework for discretizing in space

What about $ u_t = u_{xx} + f $?
1. Use finite differences in time to obtain a set of recursive spatial problems
2. Solve the spatial problems by the finite element method

Example: diffusion problem

$$ \begin{align*} \frac{\partial u}{\partial t} &= \dfc\nabla^2 u + f(\x, t),\quad &\x\in\Omega, t\in (0,T]\\ u(\x, 0) & = I(\x),\quad &\x\in\Omega\\ \frac{\partial u}{\partial n} &= 0,\quad &\x\in\partial\Omega,\ t\in (0,T] \end{align*} $$

A Forward Euler scheme; ideas

$$ \begin{equation*} [D_t^+ u = \dfc\nabla^2 u + f]^n,\quad n=1,2,\ldots,N_t-1 \end{equation*} $$

Solving wrt $ u^{n+1} $: $$ \begin{equation*} u^{n+1} = u^n + \Delta t \left( \dfc\nabla^2 u^n + f(\x, t_n)\right) \end{equation*} $$

$ u^n = \sum_jc_j^n\baspsi_j\, \in V $, $ u^{n+1} = \sum_jc_j^{n+1}\baspsi_j\,\in V $

Compute $ u^0 $ from $ I $

Compute $ u^{n+1} $ from $ u^n $ by solving the PDE for $ u^{n+1} $ at each time level

A Forward Euler scheme; stages in the discretization

$ \uex(\x,t) $: exact solution of the PDE problem

$ \uex^n(\x) $: exact solution of time-discrete problem (after applying a finite difference scheme in time)

$ \uex^n(\x)\approx u^n = \sum_{j\in\If}c_j^n\baspsi_j = $ solution of the time- and space-discrete problem (after applying a Galerkin method in space)

$$ \begin{equation*} \frac{\partial \uex}{\partial t} = \dfc\nabla^2 \uex + f(\x, t) \end{equation*} $$ $$ \begin{equation*} \uex^{n+1} = \uex^n + \Delta t \left( \dfc\nabla^2 \uex^n + f(\x, t_n)\right) \end{equation*} $$ $$ \uex^n \approx u^n = \sum_{j=0}^{N} c_j^{n}\baspsi_j(\x),\quad \uex^{n+1} \approx u^{n+1} = \sum_{j=0}^{N} c_j^{n+1}\baspsi_j(\x) $$ $$ R = u^{n+1} - u^n - \Delta t \left( \dfc\nabla^2 u^n + f(\x, t_n)\right)$$

A Forward Euler scheme; weighted residual (or Galerkin) principle

$$ R = u^{n+1} - u^n - \Delta t \left( \dfc\nabla^2 u^n + f(\x, t_n)\right)$$

The weighted residual principle: $$ \int_\Omega Rw\dx = 0,\quad \forall w\in W$$

results in $$ \int_\Omega \left\lbrack u^{n+1} - u^n - \Delta t \left( \dfc\nabla^2 u^n + f(\x, t_n)\right) \right\rbrack w \dx =0, \quad \forall w \in W $$

Galerkin: $ W=V $, $ w=v $

A Forward Euler scheme; integration by parts

Isolating the unknown $ u^{n+1} $ on the left-hand side: $$ \int_{\Omega} u^{n+1}v\dx = \int_{\Omega} \left\lbrack u^n + \Delta t \left( \dfc\nabla^2 u^n + f(\x, t_n)\right) \right\rbrack v\dx $$

Integration by parts of $ \int\dfc(\nabla^2 u^n) v\dx $: $$ \int_{\Omega}\dfc(\nabla^2 u^n)v \dx = -\int_{\Omega}\dfc\nabla u^n\cdot\nabla v\dx + \underbrace{\int_{\partial\Omega}\dfc\frac{\partial u^n}{\partial n}v \dx}_{=0\quad\Leftarrow\quad\partial u^n/\partial n=0} $$

Variational form: $$ \begin{equation*} \int_{\Omega} u^{n+1} v\dx = \int_{\Omega} u^n v\dx - \Delta t \int_{\Omega}\dfc\nabla u^n\cdot\nabla v\dx + \Delta t\int_{\Omega}f^n v\dx,\quad\forall v\in V \end{equation*} $$

New notation for the solution at the most recent time levels

$ u $ and u: the spatial unknown function to be computed

$ u_1 $ and u_1: the spatial function at the previous time level $ t-\Delta t $

$ u_2 $ and u_2: the spatial function at $ t-2\Delta t $

This new notation gives close correspondence between code and math

$$ \begin{equation*} \int_{\Omega} u v\dx = \int_{\Omega} u_1 v\dx - \Delta t \int_{\Omega}\dfc\nabla u_1\cdot\nabla v\dx + \Delta t\int_{\Omega}f^n v\dx \end{equation*} $$

or shorter $$ \begin{equation*} (u,v) = (u_1,v) - \Delta t (\dfc\nabla u_1,\nabla v) + \Delta t (f^n, v) \end{equation*} $$

Deriving the linear systems

$ u = \sum_{j=0}^{N}c_j\baspsi_j(\x) $

$ u_1 = \sum_{j=0}^{N} c_{1,j}\baspsi_j(\x) $

$ \forall v\in V $: for $ v=\baspsi_i $, $ i=0,\ldots,N $

Insert these in $$ (u, \baspsi_i) = (u_1,\baspsi_i) - \Delta t (\dfc\nabla u_1,\nabla\baspsi_i) + \Delta t (f^n,\baspsi_i) $$

and order terms as matrix-vector products ($ i=0,\ldots,N $): $$ \begin{equation*} \sum_{j=0}^{N} \underbrace{(\baspsi_i,\baspsi_j)}_{M_{i,j}} c_j = \sum_{j=0}^{N} \underbrace{(\baspsi_i,\baspsi_j)}_{M_{i,j}} c_{1,j} -\Delta t \sum_{j=0}^{N} \underbrace{(\nabla\baspsi_i,\dfc\nabla\baspsi_j)}_{K_{i,j}} c_{1,j} + \Delta t (f^n,\baspsi_i) \end{equation*} $$

Structure of the linear systems

$$ \begin{equation*} Mc = Mc_1 - \Delta t Kc_1 +\Delta t f \end{equation*} $$ $$ \begin{align*} M &= \{M_{i,j}\},\quad M_{i,j}=(\baspsi_i,\baspsi_j),\quad i,j\in\If\\ K &= \{K_{i,j}\},\quad K_{i,j}=(\nabla\baspsi_i,\dfc\nabla\baspsi_j),\quad i,j\in\If\\ f &= \{(f(\x,t_n),\baspsi_i)\}_{i\in\If}\\ c &= \{c_i\}_{i\in\If}\\ c_1 &= \{c_{1,i}\}_{i\in\If} \end{align*} $$

Computational algorithm

Compute $ M $ and $ K $.

Initialize $ u^0 $ by either interpolation or projection

For $ n=1,2,\ldots,N_t $:
1. compute $ b = Mc_1 - \Delta t Kc_1 + \Delta t f $
2. solve $ Mc = b $
3. set $ c_1 = c $

Initial condition:

Either interpolation: $ c_{1,j} = I(\x_j) $ (finite elements)

Or projection: solve $ \sum_j M_{i,j}c_{1,j} = (I,\baspsi_i) $, $ i\in\If $

Example using sinusoidal basis functions

$$ \begin{align} \frac{\partial u}{\partial t} &= \dfc\frac{\partial^2 u}{\partial x^2},\quad &x\in (0,L),\ t\in (0,T], \label{fem:deq:diffu:pde1D:eq}\\ u(x, 0) & = A\cos(\pi\x/L) + B\cos(10\pi x/L),\quad &x\in[0,L], \label{fem:deq:diffu:pde1D:ic}\\ \frac{\partial u}{\partial x} &= 0,\quad &x=0,L,\ t\in (0,T] \label{fem:deq:diffu:pde1D:bcN} \tp \end{align} $$ $$ \baspsi_i = \cos(i\pi x/L)\tp$$

Approximating the initial condition

$ I(x)\in V $ implies perfect approximation of the initial condition: $$ c_{1,1}=A,\quad c_{1,10}=B,$$ while $ c_{1,i}=0 $ for $ i\neq 1,10 $.

Computing the $ M $ and $ K $ matrices

Note that $ \baspsi_i $ and $ \baspsi_i' $ are orthogonal on $ [0,L] $ such that we only need to compute the diagonal elements $ M_{i,i} $ and $ K_{i,i} $! $$ M_{0,0}=L,\quad M_{i,i}=L/2,\ i>0,\quad K_{0,0}=0,\quad K_{i,i}=\frac{\pi^2 i^2}{2L},\ i>0\tp$$

Solving the equation system

$$ \begin{align*} Lc_0 &= Lc_{1,0} - \Delta t \cdot 0\cdot c_{1,0},\\ \frac{L}{2}c_i &= \frac{L}{2}c_{1,i} - \Delta t \frac{\pi^2 i^2}{2L} c_{1,i},\quad i>0\tp \end{align*} $$ $$ c_i = (1-\Delta t (\frac{\pi i}{L})^2) c_{1,i}\tp $$

We actually get a closed-form discrete solution: $$ u^n_i = A(1-\Delta t (\frac{\pi}{L})^2)^n \cos(\pi x/L) + B(1-\Delta t (\frac{10\pi }{L})^2)^n \cos(10\pi x/L)\tp$$

Comparing P1 elements with the finite difference method; ideas

P1 elements in 1D

Uniform mesh on $ [0,L] $ with cell length $ h $

No Dirichlet conditions: $ \baspsi_i=\basphi_i $, $ i=0,\ldots,N=N_n-1 $

Have found formulas for $ M $ and $ K $ at the element level

Have assembled the global matrices

Have developed corresponding finite difference operator formulas

$ M $: $ h[u + \frac{1}{6}h^2D_xD_x u]^n_i $

$ K $: $ h[\dfc D_xD_x u]^n_i $

Comparing P1 elements with the finite difference method; results

Diffusion equation with finite elements is equivalent to $$ \begin{equation*} [D_t^+(u + \frac{1}{6}h^2D_xD_x u) = \dfc D_xD_x u + f]^n_i \end{equation*} $$

Can lump the mass matrix by Trapezoidal integration and get the standard finite difference scheme $$ \begin{equation*} [D_t^+u = \dfc D_xD_x u + f]^n_i \end{equation*} $$

Discretization in time by a Backward Euler scheme

Backward Euler scheme in time: $$ [D_t^- u = \dfc\nabla^2 u + f(\x, t)]^n $$ $$ \begin{equation*} \uex^{n} - \Delta t \left( \dfc\nabla^2 \uex^n + f(\x, t_{n})\right) = \uex^{n-1} \end{equation*} $$ $$ \uex^n \approx u^n = \sum_{j=0}^{N} c_j^{n}\baspsi_j(\x),\quad \uex^{n+1} \approx u^{n+1} = \sum_{j=0}^{N} c_j^{n+1}\baspsi_j(\x)$$

The variational form of the time-discrete problem

$$ \begin{equation*} \int_{\Omega} \left( u^{n}v + \Delta t \dfc\nabla u^n\cdot\nabla v\right)\dx = \int_{\Omega} u^{n-1} v\dx + \Delta t\int_{\Omega}f^n v\dx,\quad\forall v\in V \end{equation*} $$

or $$ \begin{equation*} (u,v) + \Delta t (\dfc\nabla u,\nabla v) = (u_1,v) + \Delta t (f^n,v) \end{equation*} $$

The linear system: insert $ u=\sum_j c_j\baspsi_i $ and $ u_1=\sum_j c_{1,j}\baspsi_i $, $$ \begin{equation*} (M + \Delta t K)c = Mc_1 + \Delta t f \end{equation*} $$

Calculations with P1 elements in 1D

Can interpret the resulting equation system as $$ \begin{equation*} [D_t^-(u + \frac{1}{6}h^2D_xD_x u) = \dfc D_xD_x u + f]^n_i \end{equation*} $$

Lumped mass matrix (by Trapezoidal integration) gives a standard finite difference method: $$ \begin{equation*} [D_t^- u = \dfc D_xD_x u + f]^n_i \end{equation*} $$

Dirichlet boundary conditions

Dirichlet condition at $ x=0 $ and Neumann condition at $ x=L $: $$ \begin{align*} u(\x,t) &= u_0(\x,t),\quad & \x\in\partial\Omega_D\\ -\dfc\frac{\partial}{\partial n} u(\x,t) &= g(\x,t),\quad & \x\in\partial{\Omega}_N \end{align*} $$

Forward Euler in time, Galerkin's method, and integration by parts: $$ \begin{equation*} \int_\Omega u^{n+1}v\dx = \int_\Omega (u^n - \Delta t\dfc\nabla u^n\cdot\nabla v)\dx + \Delta t\int_\Omega fv \dx - \Delta t\int_{\partial\Omega_N} gv\ds,\quad \forall v\in V \end{equation*} $$

Requirement: $ v=0 $ on $ \partial\Omega_D $

Boundary function

$$ u^n(\x) = u_0(\x,t_n) + \sum_{j\in\If}c_j^n\baspsi_j(\x)$$ $$ \begin{align*} \sum_{j\in\If} & \left(\int_\Omega \baspsi_i\baspsi_j\dx\right) c^{n+1}_j = \sum_{j\in\If} \left(\int_\Omega\left( \baspsi_i\baspsi_j - \Delta t\dfc\nabla \baspsi_i\cdot\nabla\baspsi_j\right)\dx\right) c_j^n - \\ &\quad \int_\Omega\left( \left(u_0(\x,t_{n+1}) - u_0(\x,t_n)\right) \baspsi_i + \Delta t\dfc\nabla u_0(\x,t_n)\cdot\nabla \baspsi_i\right)\dx \\ & \quad + \Delta t\int_\Omega f\baspsi_i\dx - \Delta t\int_{\partial\Omega_N} g\baspsi_i\ds, \quad i\in\If \end{align*} $$

Finite element basis functions

$ B(\x,t_n)=\sum_{j\in\Ifb} U_j^n\basphi_j $

$ \baspsi_i = \basphi_{\nu(j)} $, $ j\in\If $

$ \nu(j) $, $ j\in\If $, are the node numbers corresponding to all nodes without a Dirichlet condition

$$ \begin{align*} u^n &= \sum_{j\in\Ifb} U_j^n\basphi_j + \sum_{j\in\If}c_{1,j}\basphi_{\nu(j)},\\ u^{n+1} &= \sum_{j\in\Ifb} U_j^{n+1}\basphi_j + \sum_{j\in\If}c_{j}\basphi_{\nu(j)} \end{align*} $$ $$ \begin{align*} \sum_{j\in\If} & \left(\int_\Omega \basphi_i\basphi_j\dx\right) c_j = \sum_{j\in\If} \left(\int_\Omega\left( \basphi_i\basphi_j - \Delta t\dfc\nabla \basphi_i\cdot\nabla\basphi_j\right)\dx\right) c_{1,j} - \\ &\quad \sum_{j\in\Ifb}\int_\Omega\left( \basphi_i\basphi_j(U_j^{n+1} - U_j^n) + \Delta t\dfc\nabla \basphi_i\cdot\nabla \basphi_jU_j^n\right)\dx \\ &\quad + \Delta t\int_\Omega f\basphi_i\dx - \Delta t\int_{\partial\Omega_N} g\basphi_i\ds, \quad i\in\If \end{align*} $$

Modification of the linear system; the raw system

Drop boundary function

Compute as if there are not Dirichlet conditions

Modify the linear system to incorporate Dirichlet conditions

$ \If $ holds the indices of all nodes $ \{0,1,\ldots,N=N_n-1\} $

$$ \begin{align*} \sum_{j\in\If} \biggl(\underbrace{\int_\Omega \basphi_i\basphi_j\dx}_{M_{i,j}}\biggr) c_j &= \sum_{j\in\If} \biggl(\underbrace{\int_\Omega \basphi_i\basphi_j \dx}_{M_{i,j}} - \Delta t\underbrace{\int_\Omega \dfc\nabla \basphi_i\cdot\nabla\basphi_j\dx}_{K_{i,j}}\biggr) c_{1,j} \\ &\quad \underbrace{+\Delta t\int_\Omega f\basphi_i\dx - \Delta t\int_{\partial\Omega_N} g\basphi_i\ds}_{f_i},\quad i\in\If \end{align*} $$

Modification of the linear system; setting Dirichlet conditions

$$ \begin{equation*} Mc = b,\quad b = Mc_1 - \Delta t Kc_1 + \Delta t f \end{equation*} $$

For each $ k $ where a Dirichlet condition applies, $ u(\xno{k},t_{n+1})=U_k^{n+1} $,

set row $ k $ in $ M $ to zero and 1 on the diagonal: $ M_{k,j}=0 $, $ j\in\If $, $ M_{k,k}=1 $

$ b_k = U_k^{n+1} $

Or apply the slightly more complicated modification which preserves symmetry of $ M $

Modification of the linear system; Backward Euler example

Backward Euler discretization in time gives a more complicated coefficient matrix: $$ \begin{equation*} Ac=b,\quad A = M + \Delta t K,\quad b = Mc_1 + \Delta t f \end{equation*} $$

Set row $ k $ to zero and 1 on the diagonal: $ M_{k,j}=0 $, $ j\in\If $, $ M_{k,k}=1 $

Set row $ k $ to zero: $ K_{k,j}=0 $, $ j\in\If $

$ b_k = U_k^{n+1} $

Observe: $ A_{k,k} = M_{k,k} + \Delta t K_{k,k} = 1 + 0 $, so $ c_k = U_k^{n+1} $

Analysis of the discrete equations

The diffusion equation $ u_t = \dfc u_{xx} $ allows a (Fourier) wave component $$ \begin{equation*} u = \Aex^n e^{ikx},\quad \Aex = e^{-\dfc k^2\Delta t} \end{equation*} $$

Numerical schemes often allow the similar solution $$ \begin{equation*} u^n_q = A^n e^{ikx} \end{equation*} $$

$ A $: amplification factor to be computed

How good is this $ A $ compared to the exact one?

Handy formulas

$$ \begin{align*} [D_t^+ A^n e^{ikq\Delta x}]^n &= A^n e^{ikq\Delta x}\frac{A-1}{\Delta t},\\ [D_t^- A^n e^{ikq\Delta x}]^n &= A^n e^{ikq\Delta x}\frac{1-A^{-1}}{\Delta t},\\ [D_t A^n e^{ikq\Delta x}]^{n+\half} &= A^{n+\half} e^{ikq\Delta x}\frac{A^{\half}-A^{-\half}}{\Delta t} = A^ne^{ikq\Delta x}\frac{A-1}{\Delta t},\\ [D_xD_x A^ne^{ikq\Delta x}]_q &= -A^n \frac{4}{\Delta x^2}\sin^2\left(\frac{k\Delta x}{2}\right) \end{align*} $$

Amplification factor for the Forward Euler method; results

Introduce $ p=k\Delta x/2 $ and $ C=\dfc\Delta t/\Delta x^2 $: $$ A = 1 - 4C\frac{\sin^2 p}{1 - \underbrace{\frac{2}{3}\sin^2 p}_{\hbox{from }M}}$$

(See notes for details)

Stability: $ |A|\leq 1 $: $$ \begin{equation*} C\leq \frac{1}{6}\quad\Rightarrow\quad \Delta t\leq \frac{\Delta x^2}{6\dfc} \end{equation*} $$

Finite differences: $ C\leq {\half} $, so finite elements give a stricter stability criterion for this PDE!

Amplification factor for the Backward Euler method; results

Coarse meshes: $$ A = \left( 1 + 4C\frac{\sin^2 p}{1 + \frac{2}{3}\sin^2 p}\right)^{-1} \hbox{ (unconditionally stable)} $$