Analysis of Finite Difference Schemes

Reminder on Fourier Analysis on ${\mathbb{R}}$

1. Continuous Case

   Fourier Transform Pair $\left\{\begin{array}{ll} \hat {u}(\omega)=\displaystyle \frac{1}{\sqrt{2\pi}}\... ...}\int^{\infty}_{-\infty}e^{i\omega x} \hat{u}(\omega)d\omega \end{array}\right.$

2. on a grid of integers ${\mathbb{Z}}$ or $h{\mathbb{Z}}=\{hm:m \in{\mathbb{Z}}\}$

$$\left\{ \begin{array}{ll} \hat{u}(\xi)=\frac{1}{\sqrt{2\pi}}\... ...}\int^{\pi}_{-\pi}e^{im\xi} \hat{u}(\xi)d\xi \end{array}\right.$$

Sometimes more convenient to express as

$$\hat{u}(\xi)=\frac{1}{\sqrt{2\pi}}\displaystyle \sum^{\infty}_{m=-\infty}e^{-im h \xi} v_m h \quad\quad \xi\in [-\pi/h, \pi/h]$$

$$u_m=\displaystyle \frac{1}{\sqrt{2\pi}}\int^{\frac{\pi}{h}}_{\frac{-\pi}{h}}e^{imu\xi}\hat{u}(\xi)d\xi$$

definition $\vert\vert u\vert\vert _2 =\displaystyle \sqrt{\int^{\infty}_{-\infty}\vert\vert u(x)\vert^2dx}$, is the $L^2$-norm.

Parseval's Theorem states that

$$\int^{\infty}_{-\infty}\vert u(x)\vert^2dx=\int^{\infty}_{-\infty}\vert\hat{u} (\omega)\vert^2d\omega$$

and for the grid functions:

$$\vert\vert\hat{u}\vert\vert^2_h=\int^{\pi/ h}_{-\pi/h}\vert\hat{u... ...vert^2d\xi=\sum^{\infty}_{-\infty}\vert u_m\vert^2h=\vert\vert u\vert\vert^2_h$$

Von-Neumann Analysis of Finite Difference Schemes

Use of Fourier methods in the analysis of finite difference schemes. Fourier methods for the analysis of finite difference schemes are very useful due to their simplicity. They are applicable on all linear problems and somewhat applicable for nonlinear problems. Briefly, the idea is as follows: the finite-dimensional approximation $u^n_m$ on the lattice $(mh, nk)$ of the function $u(mh, nk)$ is decomposed into a superposition of normalized sines and cosines with wave numbers $\xi$ in the range $$\frac{-\pi}{h}$$ to $$\frac{\pi}{h}$$. Thus, each sine/cosine wave is of the form

$$\hat u(\xi)e^{ix_m\xi}$$

where $\hat u(\xi)$ is the complex amplitude of the $\xi^{th}$ wave. If $u$ depends on both space $x_m$ and on time $t_n$, then the complex amplitude $\hat u(\xi)$ depends on time, hence $\hat u^n(\xi)e^{ix_m\xi}$ is the $\xi^{th}$ wave component of $u$ at time $t_n=nk$. Note that the phase $e^{ix_m\xi}$ is of magnitude 1. Hence, for stability all we have to study is each time-dependent complex amplitude $\hat u^n(\xi)$.

Take

$$\frac{u_m^{n+1}-u^n_m}{k}+a\frac{u^n_m-u^n_{m-1}}{h}=0$$

$$\mbox {or } u_m^{n+1}=(1-a\lambda) u^n_m+a\lambda u^n_{m-1}\quad \lambda =k/ h$$

   take $$u^n_m=\frac{1}{\sqrt{2\pi}}\int^{\pi/h}_{-\pi/h}e^{imh\xi} \hat{u}^n(\xi)d\xi$$

   and substitute in the above equation:

$$u^{n+1}_{m}=\frac{1}{\sqrt{2\pi}}\in... ...ambda)+a\lambda e^{-ih\xi}]\hat{u}^n(\xi)d\xi} _{\therefore\quad \hat{u}^{n+1}}$$

$$\left.\begin{array}{c} \hat{u}^{n+1}(\xi)=\underbrace{\left[(1-a\... ...{\mbox {\lq\lq amplification}\atop\mbox{factor''}} \hat u^n(\xi) \end{array}\right\}$$

amplification factor because its magnitude is the amount that the complex amplitude of each wave in the solution, is amplified in advancing one step in time.

In fact, from the above expression we obtain

$$\hat{u}^n(\xi)=g(h\xi)^n\hat{u}^0(\xi)$$

$$\mbox {So returning to } \frac{u^{n+1}_m-u^n_m}{k}+a\frac{u_m^n-u^n_{m-1}}{h}=0$$ (125)

First, by Perseval's Theorem:

$$h\sum^{\infty}_{m=-a}\vert u^n_m\vert^2 = \int^{\pi/h}_{-\pi/h}\vert\hat{u}^n(\xi)\vert^2d \xi$$

$$= \int^{\pi/h}_{-\pi/h}\vert g(h\xi)\vert^{2n}\vert \hat {u}^o(\xi)\vert^2d\xi$$ (126)

Recall that for first order single step scheme we require that

$$h\sum^{\infty}_{m=-\infty}\vert u^n_m\vert^2\le C_Th\sum^{\infty}_{m=-\infty}\vert u^0_m\vert^2$$ (127)

$$\mbox {for } 0\le nk\le T\qquad \mbox {where }\qquad t_n=0, k, 2k,\ldots$$

for stability.

Comparison of (127) and (128) implies that $\vert g(h\xi)\vert^{2n}$ must be suitably bounded. For (126):

$$u^{n+1}_m=u^n_m-a\lambda(u^n_m-u^n_{m-1})$$

enough to consider

$$\hat u^n_m=e^{imh\xi} \hat u^n(\xi)$$$$\quad \mbox {(is a wave component of solution)}$$

since linear

$$\hat{u}^{n+1}(\xi)e^{imh\xi}= \hat{u}^n(\xi)e^{imh\xi}-a\lambda \left(\hat {u}^n_ne^{im h\xi} -\hat{u}^m e^{i(m-1)h\xi}\right)$$

$$u^{n+1}(\xi)=\underbrace{[1-a\lambda(1-e^{-i\theta})]}_{g(h\xi)} \hat{u}^n(\xi)\qquad \theta =h\xi$$

$$\vert g(\theta)\vert^2=1-4a\lambda(1-a\lambda)\sin^2\left(\frac12\theta\right)$$

$$\mbox { if }\vert g(\theta)\vert^2\le 1\quad \mbox {then }\quad 0\le a \lambda\le 1\quad \therefore$$

$$h\sum^{\infty}_{m=-\infty}\vert u^n_m\vert^2\le \int^{\pi/h}_{-\pi/h}\vert u^0(\xi)\vert^ 2d\xi=h\sum^{\infty}_{m=-\infty}\vert u^0_m\vert^2$$

$\Box$

For this particular example the amplification factor $g$ depended on $\theta=h\xi$ only, but in general it can depend on $h$ and $k$.

Theorem (Stability, Von Newmann): A one-step constant coefficient scheme is stable if and only if $\exists$ a constant $K$ (independent $\theta$, $k$, and $h$) and some positive grid spacings $k_0$ and $h_0$ such that

$$\vert g(\theta, k, h)\vert\le 1+Kk$$ (128)

     $\forall \theta, 0<k<k_0, 0<h<h_0$. If $g$ is independent of $h$, and $k$ then (129) is replaced by

$$\vert g(\theta)\vert\le 1$$

$\Box$

Example

$$U_t+aU_x-U=0\qquad ( \mbox {has solutions that grow with }t) \mbox { using Lax-Friedrichs:}$$

$$\frac{u^{n+1}_m-\frac12(u^n_{m+1}+u^n_{m-1})}{k}+a\frac{u^n_{m+1}-u^n _{m-1}}{2h}-u^n_m=0$$

$$g(\theta, k, h)=\cos\theta -ia\lambda \cos\theta +k$$

$$\vert g\vert^2=(\cos\theta+k)^2+a^2\lambda^2\sin^2\theta \le (1+k)^2 \vert a\lambda\vert\le 1$$

$$\vert g\vert^2 1 \theta.$$

Proof: By Parseval's Theorem

$$\vert\vert u^n\vert\vert^2_h=\int^{\pi/h}_{-\pi/2}\vert g(h\xi, k, h)\vert^{2n} \vert u^0(\xi)\vert^2d\xi$$

$$\vert g(h\xi, k,h)\vert\le 1+Kk\therefore$$

$$\vert\vert u^n\vert\vert^2_h\le \int^{\pi/h}_{-\pi/h}(1+Kk)^{2n}\vert\hat{u}^0(\xi)\vert^2d\xi = (1+Kk)^{2n}\vert\vert u^0\vert\vert^2_h$$

Now $n\le T/k$ so      $(1+Kk)^n\le (1+Kk)^{T/k}\le e^{kT}$

$\therefore \vert\vert u^n\vert\vert _n\le e^{Kt}\vert\vert u^0\vert\vert _h$ which is the condition

$$h\sum^{\infty}_{m=-\infty}\vert u^n_m\vert^2\le C_Th\sum^{J}_{j=0} \sum^{\infty}_{m=\infty}\vert u^j_m\vert^2 0\le nk\le T$$ (129)

Now we prove that if inequality (129) cannot be satisfied for any value of $K\Rightarrow$ scheme is not stable. To do so, we can show that any amount of growth in the solution, that is, we show that the stability inequality (130) cannot hold.

If for some positive value $C$ there's an interval of $\theta$'s, $\theta\in [\theta_1, \theta_0]$ and $h\in(0,h_0]$ and $k\in(0,k_0]$ with $\vert g(\theta,k,h)\vert\ge1+Ck$ then we construct a function $v^0_m$ as

$$\hat{u}^0(\xi)=\left\{\begin{array}{ll} 0 & \mbox {if } h\xi\noti... ..._2-\theta_1)^{-1}} & \mbox{if } h\xi\in [\theta_1, \theta_2] \end{array}\right.$$

then

$$\vert\vert u^n\vert\vert^2_h = \int^{\pi/L}_{-\pi/h}\vert g(h\xi, k, h)\vert^{2n}\vert\hat{u}^0(\xi)\vert^2 d\xi$$

$$= \int^{\theta_2/h}_{\theta_1/h}\vert g(h\xi, k, h)\vert^{2n}\frac{h}{\theta_2-\theta_1}d\xi\ge (1+Ck)^{2n}$$

$$\ge \frac12 e^{2TC}\vert\vert u^0\vert\vert^2_h \mbox { for }n \mbox{ near } T/k.$$

This shows that the scheme to be unstable if $C$ can be arbitrarily large.

$\Box$

Corollary: If a scheme as in previous theorem is modified so that the modifications result only in the addition to the amplification factor of terms that are ${\cal O}(k)$ uniformly in $\xi$, then the modified scheme is stable if and only if the original scheme is stable.

Proof: If $g$ is the amplification factor for the scheme and satisifies $\vert g\vert\le1+Kk$, then the amplification factor or the modified scheme $g'$ satisfies

$$\vert g'\vert=\vert g+{\cal O}(k)\vert\le 1+Kk+Ck=1+K'k$$ (130)

Hence the modified scheme is stable if the original scheme is stable and vice versa. $\Box$

Stability for variable coefficients

Take $U_t+a(x,t)U_x=0$ as an example.

The general procedure is to consider the problem with $a$ as a frozen coefficient for each $x,t$ values in question. If each frozen coefficient case is stable then the scheme is stable. For example, the CFL condition would require

$$\vert a(t_n, x_m)\vert\lambda\le 1$$

for all $t_n, x_m$ in computational lattice.

Remark: Numerical vs Dynamic Stability $\rightarrow$ Numerical stability refers to the behavior of approximations to a grid projected equation over a finite time interval as the lattice is refined. Dynamic stability refers to the behavior of solutions of PDE as $t\rightarrow\infty$.

Example

$$\begin{array}{ll} u_t+au_x+bu=0 & x\in {\mathbb{R}^{1}}\\ & , t>0 \end{array}$$

for $b>0$, solution is dynamically stable (bounded) since solution decays as $t\rightarrow\infty$. (Show this using Fourier methods).

for $b>0$, solution is unstable.

For above equation a stable numerical scheme would be one in which the approximation converges to exact solution for any $b$ as $k,h\rightarrow 0$.

$\Box$

Comments on Instabilities in Hyperbolic Equation Approximation

1)

As always, non-convergent schemes are useless.

2)

Take Lax-Friedrichs: $\vert g(\theta)\vert^2=\cos^2\theta+a^2\lambda^2\sin^3\theta$. The maximum value of $g$ is attained when $\theta=\displaystyle \frac{\pi}{2}$, where $\vert g\vert=1.6\Rightarrow$ so instabilities are related to high frequency oscillations.

3)

In general (not a theorem) $\Rightarrow$ instabilities will manifest themselves as rapid growth of high wave numbers and would then be more evident if the initial data contains high amplitude high wave number modes (e.g. non-smooth data). Also, the instabilities will have wavelengths that are comparable or commensurate to the grid space in $x$. Also, the instability phenomenom will be local in nature but propagate in time. Of course, it will eventually swamp the solution, but for close times after the onset of the instability, they are local.

4)

Having a good feel for this allows one to discern between programming errors and improper schemes.

5)

The other considerations related to stability are that schemes with very restrictive conditions on step sizes, in particular, on time steps, will require more computational effort in calculation. This is ok, but aside from stability, we also need to consider the dissipation and dispersion in the scheme.

6)

Stability of a scheme requires that this issue be checked carefully at the boundaries. One of the most common sources of instabilities comes from using inappropriate boundary conditions.

Example

Take Lax-Friedrichs on

$$U_t+aU_{xxx}=f.$$

The finite difference approximation is

$$u^{n+1}_{m}=\frac12(u^n_{m+1}+u^n_{m-1})-\frac12akh^{-3}(u^n_{m+2}-2 u ^{n}_{m+1}+ 2 u^{n}_{m-1}-u^{n}_{m-2})+kf^n_m$$

this scheme is consistent if $h^2/k\rightarrow 0$ as $h$ and $k$ tend to 0.

The amplification factor in this case is

$$g(\theta)=\cos\theta+i4akh^{-3}\sin \theta\sin^2\frac{1}{2}\theta$$

so scheme is stable if $$\frac{4\vert a\vert k}{h^3}$$ is bounded.

However the consistency condition $\frac{h^2}{k}\rightarrow 0$ as $h$ and $k \rightarrow 0$ and stability condition $$\frac{4\vert a\vert k}{h^3}$$ bounded cannot be both satisfied.

$\therefore$    Scheme is not convergent. $\Box$

Truncation Error and Order of Accuracy for FD Schemes

definition: A scheme $P_{k,h}u=R_{k,h}f$ that is consistent with $PU=f$ is accurate of order $p$ in time and $q$ in space if for any smooth $\phi(t,x)$

$$P_{k,h}\phi-R_{k,h}P\psi={\cal 0}(k^p)+{\cal O}(h^q)$$

We say scheme is order ($p,q$). $\Box$

Example Crank-Nicholson

Take

$$U_t=\frac{U(t+k,x)-U(t,x)}{k}+{\cal O}(k^2)$$

for

$$U_t(t+\frac12k, x)$$

Take

$$U_x=(t+\frac12k,x) = \frac{U_x(t+\frac12k,x)+U_x(t,x)}{2}+{\cal O}{k^2}$$

$$= \frac12\left[\frac{U(t+k,x+h)-U(t+k, x-h)}{2h}+ \frac{U(t,x+h) -U(t,x-h)}{2h}\right]$$

$$+ {\cal O}(k^2)+{\cal O}(h^2)$$

Using these approximate $U_t+a U_x=f$ about $\left(t+\frac12k, x\right)$ we get

$$\frac{u^{n+1}_m-u^{n}_m}{k}+a\frac{u^{n+1}_{m+1} -u^{n+1}_{m-1}+u^{n}_{m+1}-u^{n}_{m-1}}{4h}= \frac{f_m^{n+1}+f^n_m}{2}$$

And is an order $(2,2)$ scheme.

Exercise

For $U_t + a U_x =0$ show that Crank-Nicholson has an amplification factor

$$g(\theta)=\frac{1-i\frac12a\lambda\sin\theta} {1+i\frac12a\lambda\sin\theta}$$$$\mbox { and is unconditionally stable.}$$

Exercise

Show that Lax-Wendroff is consitent with $U_t+a U_x=f$:

$$g(\theta)=1-2a^2\lambda^2\sin^2\frac12\theta-ia\lambda\sin\theta$$

$$\vert a\lambda\vert\le 1 \mbox { and of order } (1,2)$$

$$\frac{u_m^{n+1}-u^n_m}{k}+a\frac{u^n_{m+1}-u^n _{m-1}}{2h}-\frac{a^2k}{2h^2}(u^n_{m+1}-2u^n_m+u^n_{m-1})$$

$$=\frac12\left(f^{n+1}_m+f^n_m\right)-\frac{ak}{4h} \left(f^{n}_{m+1}-f^n_{m-1}\right)$$

$\Box$

Order of Accuracy: The choice of norm is problem-dependent. Could use our grid $L_2$ norm. Then

Error $(t)=\vert\vert U(t,\cdot)-u^n\vert\vert _h=\left(h\sum_m\vert U(t,x_m)-u^n_m\vert^2\right)^{\frac12}=0(h^r)$ gives the accuracy of the ``solution'' $u^n$ as an approximation to exact solution $U(t,x)$ on the grid. The usefulness of the above norm is that we should get the order of accuracy to be equal to order of truncation of scheme FOR SMOOTH DATA.