INF5620 Lecture: Analysis of finite difference schemes for diffusion processes
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:
- 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)
\tag{3}
\end{equation}
$$
Visualization of the damping in the diffusion equation
Damping of a discontinuity; problem and model
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:
- \( |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}
$$
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
- 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) \).