$$ \newcommand{\uex}{{u_{\small\mbox{e}}}} \newcommand{\uexd}[1]{{u_{\small\mbox{e}, #1}}} \newcommand{\vex}{{v_{\small\mbox{e}}}} \newcommand{\vexd}[1]{{v_{\small\mbox{e}, #1}}} \newcommand{\Aex}{{A_{\small\mbox{e}}}} \newcommand{\half}{\frac{1}{2}} \newcommand{\halfi}{{1/2}} \newcommand{\tp}{\thinspace .} \newcommand{\Ddt}[1]{\frac{D #1}{dt}} \newcommand{\E}[1]{\hbox{E}\lbrack #1 \rbrack} \newcommand{\Var}[1]{\hbox{Var}\lbrack #1 \rbrack} \newcommand{\Std}[1]{\hbox{Std}\lbrack #1 \rbrack} \newcommand{\xpoint}{\boldsymbol{x}} \newcommand{\normalvec}{\boldsymbol{n}} \newcommand{\Oof}[1]{\mathcal{O}(#1)} \newcommand{\x}{\boldsymbol{x}} \newcommand{\X}{\boldsymbol{X}} \renewcommand{\u}{\boldsymbol{u}} \renewcommand{\v}{\boldsymbol{v}} \newcommand{\w}{\boldsymbol{w}} \newcommand{\V}{\boldsymbol{V}} \newcommand{\e}{\boldsymbol{e}} \newcommand{\f}{\boldsymbol{f}} \newcommand{\F}{\boldsymbol{F}} \newcommand{\stress}{\boldsymbol{\sigma}} \newcommand{\strain}{\boldsymbol{\varepsilon}} \newcommand{\stressc}{{\sigma}} \newcommand{\strainc}{{\varepsilon}} \newcommand{\I}{\boldsymbol{I}} \newcommand{\T}{\boldsymbol{T}} \newcommand{\dfc}{\alpha} % diffusion coefficient \newcommand{\ii}{\boldsymbol{i}} \newcommand{\jj}{\boldsymbol{j}} \newcommand{\kk}{\boldsymbol{k}} \newcommand{\ir}{\boldsymbol{i}_r} \newcommand{\ith}{\boldsymbol{i}_{\theta}} \newcommand{\iz}{\boldsymbol{i}_z} \newcommand{\Ix}{\mathcal{I}_x} \newcommand{\Iy}{\mathcal{I}_y} \newcommand{\Iz}{\mathcal{I}_z} \newcommand{\It}{\mathcal{I}_t} \newcommand{\If}{\mathcal{I}_s} % for FEM \newcommand{\Ifd}{{I_d}} % for FEM \newcommand{\Ifb}{{I_b}} % for FEM \newcommand{\setb}[1]{#1^0} % set begin \newcommand{\sete}[1]{#1^{-1}} % set end \newcommand{\setl}[1]{#1^-} \newcommand{\setr}[1]{#1^+} \newcommand{\seti}[1]{#1^i} \newcommand{\sequencei}[1]{\left\{ {#1}_i \right\}_{i\in\If}} \newcommand{\basphi}{\varphi} \newcommand{\baspsi}{\psi} \newcommand{\refphi}{\tilde\basphi} \newcommand{\psib}{\boldsymbol{\psi}} \newcommand{\sinL}[1]{\sin\left((#1+1)\pi\frac{x}{L}\right)} \newcommand{\xno}[1]{x_{#1}} \newcommand{\Xno}[1]{X_{(#1)}} \newcommand{\yno}[1]{y_{#1}} \newcommand{\Yno}[1]{Y_{(#1)}} \newcommand{\xdno}[1]{\boldsymbol{x}_{#1}} \newcommand{\dX}{\, \mathrm{d}X} \newcommand{\dx}{\, \mathrm{d}x} \newcommand{\ds}{\, \mathrm{d}s} \newcommand{\Real}{\mathbb{R}} \newcommand{\Integerp}{\mathbb{N}} \newcommand{\Integer}{\mathbb{Z}} $$
INF5620 Lecture: Analysis of finite difference schemes for diffusion processes

INF5620 Lecture: Analysis of finite difference schemes for diffusion processes

Hans Petter Langtangen [1, 2]

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

Dec 14, 2013

Table of contents

Analysis of schemes for the diffusion equation
      Properties of the solution
      Example
      Visualization of the damping in the diffusion equation
      Damping of a discontinuity; problem and model
      Damping of a discontinuity; Backward Euler simulation
      Damping of a discontinuity; Forward Euler simulation
      Damping of a discontinuity; Crank-Nicolson simulation
      Fourier representation
      Analysis of the finite difference schemes
      Analysis of the Forward Euler scheme
      Results for stability
      Analysis of the Backward Euler scheme
      Stability
      Analysis of the Crank-Nicolson scheme
      Stability
      Summary of accuracy of amplification factors; large time steps
      Summary of accuracy of amplification factors; time steps around the Forward Euler stability limit
      Summary of accuracy of amplification factors; small time steps
      Observations

Analysis of schemes for the diffusion equation

Properties of the solution

The PDE

 
$$ \begin{equation} u_t = \dfc u_{xx} \tag{1} \end{equation} $$

 
admits solutions

 
$$ \begin{equation} u(x,t) = Qe^{-\dfc k^2 t}\sin\left( kx\right) \tag{2} \end{equation} $$

 

Observations from this solution:

Example

Test problem:

 
$$ \begin{align*} u_t &= u_{xx},\quad &x\in (0,1),\ t\in (0,T]\\ u(0,t) &= u(1,t) = 0,\quad &t\in (0,T]\\ u(x,0) & = \sin (\pi x) + 0.1\sin(100\pi x) \end{align*} $$

 

Exact solution:

 
$$ \begin{equation} u(x,t) = e^{-\pi^2 t}\sin (\pi x) + 0.1e^{-\pi^2 10^4 t}\sin (100\pi x) \tag{3} \end{equation} $$

 

Visualization of the damping in the diffusion equation

Damping of a discontinuity; problem and model

Problem. Two pieces of a material, at different temperatures, are brought in contact at \( t=0 \). Assume the end points of the pieces are kept at the initial temperature. How does the heat flow from the hot to the cold piece?

Solution.

Assume a 1D model is sufficient (insulated rod):

 
$$ u(x,0)=\left\lbrace \begin{array}{ll} U_L, & x < L/2\\ U_R,& x\geq L/2 \end{array}\right. $$

 

 
$$ \frac{\partial u}{\partial t} = \dfc \frac{\partial^2 u}{\partial x^2},\quad u(0,t)=U_L,\ u(L,t)=U_R$$

 

Damping of a discontinuity; Backward Euler simulation

Damping of a discontinuity; Forward Euler simulation

Damping of a discontinuity; Crank-Nicolson simulation

Fourier representation

Represent \( I(x) \) as a Fourier series

 
$$ \begin{equation} I(x) \approx \sum_{k\in K} b_k e^{ikx} \end{equation} $$

 

The corresponding sum for \( u \) is

 
$$ \begin{equation} u(x,t) \approx \sum_{k\in K} b_k e^{-\dfc k^2t}e^{ikx} \tag{4} \tp \end{equation} $$

 

Such solutions are also accepted by the numerical schemes, but with an amplification factor \( A \) different from \( \exp{({-\dfc k^2t})} \):

 
$$ \begin{equation} u^n_q = A^n e^{ikq\Delta x} = A^ne^{ikx} \tag{5} \end{equation} $$

 

Analysis of the finite difference schemes

Stability:

Accuracy:

$$

Analysis of the Forward Euler scheme

 
$$ \begin{equation*} [D_t^+ u = \dfc D_xD_x u]^n_q \end{equation*} $$

 

Inserting

 
$$ u^n_q = A^n e^{ikq\Delta x}$$

 
leads to

 
$$ \begin{equation} A = 1 -4C\sin^2\left( \frac{k\Delta x}{2}\right),\quad C = \frac{\dfc\Delta t}{\Delta x^2} \end{equation} $$

 

The complete numerical solution is

 
$$ \begin{equation} u^n_q = (1 -4C\sin^2 p)^ne^{ikq\Delta x},\quad p = k\Delta x/2 \end{equation} $$

 

Results for stability

We always have \( A\leq 1 \). The condition \( A\geq -1 \) implies

 
$$ 4C\sin^2p\leq 2 $$

 
The worst case is when \( \sin^2 p=1 \), so a sufficient criterion for stability is

 
$$ \begin{equation} C\leq {\half} \end{equation} $$

 
or:

 
$$ \begin{equation} \Delta t\leq \frac{\Delta x^2}{2\dfc} \end{equation} $$

 

Implications of the stability result. Less favorable criterion than for \( u_{tt}=c^2u_{xx} \): halving \( \Delta x \) implies time step \( \frac{1}{4}\Delta t \) (not just \( \half\Delta t \) as in a wave equation). Need very small time steps for fine spatial meshes!

Analysis of the Backward Euler scheme

 
$$ \begin{equation*} [D_t^- u = \dfc D_xD_x u]^n_q\end{equation*} $$

 

 
$$ u^n_q = A^n e^{ikq\Delta x}$$

 

 
$$ \begin{equation} A = (1 + 4C\sin^2p)^{-1} \tag{6} \end{equation} $$

 

 
$$ \begin{equation} u^n_q = (1 + 4C\sin^2p)^{-n}e^{ikq\Delta x} \end{equation} $$

 

Stability

We see from (6) that \( |A|<1 \) for all \( \Delta t>0 \) and that \( A>0 \) (no oscillations).

Analysis of the Crank-Nicolson scheme

The scheme

 
$$ [D_t u = \dfc D_xD_x \overline{u}^x]^{n+\half}_q$$

 
leads to

 
$$ \begin{equation} A = \frac{ 1 - 2C\sin^2p}{1 + 2C\sin^2p} \end{equation} $$

 

 
$$ \begin{equation} u^n_q = \left(\frac{ 1 - 2C\sin^2p}{1 + 2C\sin^2p}\right)^ne^{ikp\Delta x} \end{equation} $$

 

Stability

The criteria \( A>-1 \) and \( A<1 \) are fulfilled for any \( \Delta t >0 \).

Summary of accuracy of amplification factors; large time steps

Summary of accuracy of amplification factors; time steps around the Forward Euler stability limit

Summary of accuracy of amplification factors; small time steps

Observations

The problems of correct damping for \( u_t = u_{xx} \) is partially manifested in the similar time discretization schemes for \( u'(t)=-\dfc u(t) \).

← →