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}
\label{diffu:pde1:eq}
\end{equation}
$$
admits solutions
$$
\begin{equation}
u(x,t) = Qe^{-\dfc k^2 t}\sin\left( kx\right)
\label{diffu:pde1:sol1}
\end{equation}
$$
Observations from this solution:
The initial shape \( I(x)=Q\sin kx \)
undergoes a damping \( \exp{(-\dfc k^2t)} \)
The damping is very strong for short waves (large \( k \))
The damping is weak for long waves (small \( k \))
Consequence: \( u \) is smoothened with time
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)
\label{diffu:pde1:sol2}
\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}
\label{diffu:pde1:u:Fourier}
\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}
\label{diffu:pde1:analysis:uni}
\end{equation}
$$
Analysis of the finite difference schemes
Stability:
\( |A|< 1 \): decaying numerical solutions (as we want)
\( A<0 \): oscillating numerical solutions (as we do not want)
Accuracy:
Compare numerical and exact amplification factor: \( A \) vs \( \Aex = \exp{(-\dfc k^2 \Delta t)} \)
$$
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}
\label{diffu:pde1:analysis:BE:A}
\end{equation}
$$
$$
\begin{equation}
u^n_q = (1 + 4C\sin^2p)^{-n}e^{ikq\Delta x}
\end{equation}
$$
Stability
We see from \eqref{diffu:pde1:analysis:BE:A} 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
Crank-Nicolson gives oscillations and not much damping of short waves
for increasing \( C \).
These waves will manifest themselves as high frequency
oscillatory noise in the solution.
All schemes fail to dampen short waves enough
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) \).
