Lecture 3: Finite Difference Method

To be learned:

• Approximation + Stability = Convergence
• 2D elliptic equation (a.k.a. ${\rm div}(k \nabla u) = f$)
• 2D diffusion equation
  • explicit method
  • implicit method
• Application to Image Processing

Approximation + Stability = Convergence

• Exact problem: $A u = f$
• Approximate problem: $A_h u_h = f_h$
• Approximation: $\|A_h u - f_h\| = \epsilon_h$
  • we had earlier that $\epsilon_h = O(h^2)$

• Stability (hard to prove for nontrivial equations): $\|A_h^{-1}\| \leq C$ as $h\to 0$.

• Then $\|u_h - u\| \leq C \|A_h u_h - A_h u\| = C \|f - A_h u\| = C \epsilon_h~~~$ ($=O(h^2)$ in our case)

• Exercise: how did we do the first transition?

How does it work for a 1D ODE?

• Exact problem: $ u'' = f, \quad u(0)=u(1) = 0 $

• Approximate problem: $$ \begin{align*} h^{-2} (u _h(x+h)-2u _h(x)+u _h(x-h)) &= f(x), \\ u _h(0) = u _h(1) &= 0 \end{align*} $$ on a grid with $h=1/N$

• Approximation (we need to substitute $u$ into the approximate problem): $$ \begin{align*} & h^{-2} (u(x+h)-2u(x)+u(x-h)) - f(x) \\&= u''(x) - f(x) + O(h^2) = O(h^2) \end{align*} $$ (easy!)

• Stability... gotta work hard...

Proving stability

• $A_h$ is symmetric (disagreers will be suspected of not doing the homework! :) )

• $\|A_h^{-1}\| \leq C$ $\Leftrightarrow$ $|\lambda_\min(A_h)| \geq C^{-1}$, where $\lambda_\min$ is the eigenvalue with the smallest absolute value

• Exercise: Do some math:

  • Prove that $v^{(k)}(x) := \sin(k \pi x)$ is an eigenvector
  • Find its eigenvalue $\lambda^{(k)}$
  • What is an easy way to prove that $v^{(k)}$ are linearly independent for $k=1,...,N-1$?
  • Can you now prove that $\|A_h^{-1}\| \leq C$?

2D Elliptic Equation

• A general diffusion equation is $$ {\rm div}(k\, \nabla u) = f $$
• Here $k$ is typically a scalar diffusion coefficient
• Sometimes $k$ is a matrix (diffusion tensor)

Discretization

• If $k={\rm const}$ then we know how to do it: $$ {\rm div}(k\, \nabla u) = k \Delta u = k \Delta_h u $$

• In a tabular form we write it as $$ k \Delta_h \sim \begin{pmatrix} 0 & k & 0 \\ k & -4k & k \\ 0 & k & 0 \end{pmatrix} $$

• If $k\ne{\rm const}$ then we shall be more careful: $$\scriptstyle {\rm div}(k\, \nabla u)(x,y) = \left(\begin{array}{ccc} 0 & k(x,y-h/2) & 0 \\ k(x-h/2,y) & \text{guess} & k(x+h/2,y) \\ 0 & k(x,y+h/2) & 0 \end{array}\right) $$

• Let's see how it works in 1D:

  • $k(x+h/2) u'(x+h/2) = k(x+h/2) (u(x+h)-u(x)) h^{-1} + O(h^2)$ (Exercise: prove)
  • Likewise for $k(x-h/2) u'(x-h/2)$
  • Hence $$\textstyle \frac{k(x+h/2) (u(x+h/2)-u(x)) - k(x-h/2) (u(x)-u(x-h))}{h^2} = (k u')'+ O(h^2) $$
• Exercise: work it out in 2D
• Exercise: attempt to work it out in 2D when $k=k(x)$ is a 2x2 matrix

Nicer notation

• If we denote the FD operators $D_+ u = u(x+h) - u(x)$ and $D_- u = u(x) - u(x-h)$, then

$$ {\rm div}(k\nabla u)(x) = D_- \big(k(x+h/2) \, D_+ u(x)\big) $$

Time-dependent: 2D diffusion

$$ \begin{align*} u _t &= k (u _{xx} + u _{yy})\\ u| _{t=0} &= u _0 \\ u| _{\Gamma} &= 0 \end{align*} $$

• Approximation is typically easy: just replace with discrete operators:

$$ \begin{align*} \frac{u _h(t+\tau)-u _h(t)}{\tau} &= k \Delta _h u _h(t) \qquad \text{for }(x,y)\in\Omega \\ u _h| _{t=0} &= u _0 \\ u _h| _{\Gamma} &= 0 \end{align*} $$

• This is an explicit discretization: we don't need to invert $\Delta_h$ (also called the "forward difference scheme")

• Exercise: Order of approximation: $O(\tau + h^2)$

• Stability: again not easy

Stability

• We want the result in the form $\|u_h\| \leq C \|u_0\|$
• We know that $\|u(t)\|$ decreases in time (energy diffuses $\Rightarrow$ not growing)
• Hence we require $\|u_h(t+\tau)\| \leq \|u_h(t)\|$
• In other words, eigenvalues of $I + \tau k\Delta_h$ are between $-1$ and $1$
• Exercise: derive the stability condition on $\tau$ (and $h$ and $k$)

So that you know: Implicit discretization

$$ \begin{align*} \frac{u_h(t+\tau)-u_h(t)}{\tau} &= k \frac{\Delta_h u_h(t+\tau) + \Delta_h u_h(t)}{2} \qquad \text{for }(x,y)\in\Omega \\ u_h|_{t=0} &= u_0 \\ u_h|_{\Gamma} &= 0 \end{align*} $$

• This is called the "Crank-Nicolson" scheme

• Here we do need to invert $\Delta_h$: $$ u_h(t+\tau) = (I - \frac{1}{2} \tau \Delta_h)^{-1} (I + \frac{1}{2} \tau \Delta_h) u_h(t) $$

• Exercise: approximation error is $O(\tau^2 + h^2)$

• Exercise: stability is unconditional (i.e., always stable!)

• If we can efficiently invert $(I - \frac{1}{2} \tau \Delta_h)^{-1}$ then it is better than the explicit!

Cool application: Image processing

• Let's solve $$ \begin{align*} \frac{u_h(t+\tau)-u_h(t)}{\tau} &= k \Delta_h u_h(t) \qquad \text{for }(x,y)\in\Omega \\ u_h|_{t=0} &= u_0 \\ u_h|_{\Gamma} &= u_0|_{\Gamma} \end{align*} $$ Here $u_0(x,y)$ is the image we want to process

"Anisotropic" diffusion

• Let's solve $$ \begin{align*} u' &= {\rm div} (k \nabla u) \qquad \text{for }(x,y)\in\Omega \\ u|_{t=0} &= u_0 \\ u|_{\Gamma} &= u_0|_{\Gamma} \end{align*}, $$ where $k=k(\nabla u) = 1/\sqrt{1+|\nabla u|^2}$ ($|v|$ is the length of $v$)

• We use the explicit finite difference method:

• We discretize it with $\tau=0.25$ below:

To install the JSAnimation library

In order to view animation of noise removal you need to install JSAnimation library. The easiset way to do this is to run from command line:

conda install --channel https://conda.binstar.org/IOOS jsanimation

Win-64 users:

conda install --channel https://conda.binstar.org/melund JSAnimation

You can explore all available packages of that library by running:

binstar search -t conda jsanimation