PDEs, Stability¶

• FTCS method
• Example
• Von Neumann stability analysis

Example¶

• See Excel example.
• PDE and FTCS method: $$\frac{\partial f}{\partial t} = \alpha\frac{\partial^2f}{\partial x^2},$$ $$f_i^{n+1} = f_i^n + \underbrace{\frac{\alpha\Delta t}{\Delta x^2}}_{d}(f_{i-1}^n-2f_i^n+f_{i+1}^n).$$

• Vary $d$, and examine the results.
• Start with a triangular profile.
  • How does the peak decay?

• Vary the parameter $d$:

• For even larger $d$, the profile will go negative, which may not even be physical.
• Note, in the PDE: $$\frac{\partial f}{\partial t} = \alpha\frac{\partial^2f}{\partial x^2} = \alpha\cdot(\mbox{curvature}).$$
  • Where there is no curvature, there is no change from one step to the next.
  • So, in the above image $f_{i-1}^{n+1} = f_{i-1}^n$, and $f_{i+1}^{n+1} = f_{i+1}^n$. $$f_{i-1}^{n+1} = f_{i-1}^n + d(f_{i-2}^n - 2f_{i-1}^n + f_i^n) = f_{i-1}^n + d(1\cdot 0 - 2\cdot 1 + 1\cdot 2) = f_{i-1}^n + d(0) = f_{i-1}^n.$$
• For large $d$, on the next timestep, the pattern repeates itself on the newly formed twin peaks to continue the oscillatory cycle.
  • If $d$ is large enough, the processes is unstable and the oscillations grow in amplitude.

Timescale¶

$$\frac{\partial f}{\partial t} = \alpha\frac{\partial^2f}{\partial x^2},$$

• Nondimensionalize this PDE.
  • $f^* = f/f_r \rightarrow f = f^*f_r.$
  • $x^* = x/L \rightarrow x = x^*L.$
  • $t^* = t/\tau \rightarrow t = t^*\tau.$
  • Insert these into the PDE $$\left[\frac{1}{\tau}\right]\left(\frac{\partial f^*}{\partial t^*}\right) = \left[\frac{\alpha}{L^2}\right]\left(\frac{\partial^2 f^*}{\partial x^{*2}}\right).$$
• The units of this equation imply that $$\tau = \frac{L^2}{\alpha}.$$

• This is the characteristic timescale for diffusion.
• If $L=\Delta x$, then $\tau=\Delta x^2/\alpha$.
• Now, we see that our parameter $d$, is a ratio of timescales: $$ d= \frac{\alpha\Delta t}{\Delta x^2} = \frac{\Delta t}{\tau}.$$
• For stability $$d\le \frac{1}{2},$$

$$\Delta t \le \frac{1}{2}\frac{\Delta x^2}{\alpha}.$$

• For stability, don't step larger than a factor (1/2) of the characteristic diffusion time across a cell of width $\Delta x$.
• Note how $\Delta t$ is tied to $\Delta x^2$.
  • If I use 10$\times$ more grid points, then I have to take 100$\times$ smaller timesteps!

Von Neumann Stability¶

• Here we show that $d\le 1/2$ is the stability criterion.
• Linear PDE with constant coefficients (as above): $$f_t = \alpha f_{xx}.$$
• Express the solution $f(x,t_n)$ as a Fourier Series. $$f(x,t_n) = \sum_{\eta=-\infty}^{\infty}c_{\eta} e^{ik_{\eta}x},$$ $$k_{\eta} = \frac{2\pi\eta}{2L},$$ where $k_{\eta}$ is the wavenumber (inverse wavelength), and $c_{\eta}$ is a complex coefficient.

• The PDE is linear, so we can consider each mode independently,
  • e.g., if $f=g+h$ then $f_t=\alpha f_{xx}$ implies $g_t+h_t = \alpha g_{xx}+\alpha h_{xx}.$
  • And our above Fourier Series gives $f$ as a sum.
• Now, our FTCS method is

$$f_i^{n+1} = f_i^n + d(f_{i-1}^n-2f_i^n+f_{i+1}^n),$$

• The Fourier Series term implies:

$$f_i = c_{\eta}e^{ik_{\eta}x},$$

$$f_{i+1} = c_{\eta}e^{ik_{\eta}(x+\Delta x)} = \underbrace{c_{\eta}e^{ik_{\eta}x}}_{f_i}e^{ik_{\eta}\Delta x}= f_ie^{ik_{\eta}\Delta x},$$

$$f_{i-1} = c_{\eta}e^{ik_{\eta}(x-\Delta x)} = \underbrace{c_{\eta}e^{ik_{\eta}x}}_{f_i}e^{-ik_{\eta}\Delta x}= f_ie^{-ik_{\eta}\Delta x},$$

• Now, subsitute these into the FTCS equation: $$f_i^{n+1} = f_i^n(1+d(e^{-ik_{\eta}\Delta x} - 2 + e^{ik_{\eta}\Delta x})),$$
• Use $$ e^{ik_{\eta}\Delta x} = \cos(k_{\eta}\Delta x) + i\sin(k_{\eta}\Delta x),$$ to get $$\frac{f_i^{n+1}}{f_i^n} = 1+2d(\cos(k_{\eta}\Delta x)-1) = G,\,\,\mbox{say}.$$

• This will be stable for $|G|\le 1$. This implies $$-1\le 1 + \underbrace{2d}_{>0}(\underbrace{\cos(k_{\eta}\Delta x)}_{-1\,\mbox{to}\,1}-1) \le 1.$$
  • The upper bound in this equation (right-most inequality) is always satisfied.
  • For the lower bound, rearranging, we have $$ -2\le2d(\cos(k_{\eta}\Delta x)-1),$$ or $$d\le\frac{1}{1-\cos(k_{\eta}\Delta x)}.$$

• We want the minimum of the right hand side (RHS), since $d$ has to be less than or equal to it. The RHS is minimized for $$\max(1-\underbrace{\cos(k_{\eta}\Delta x)}_{-1\,\mbox{to}\,1}),$$
  • The maximum here is $2$.
  • Hence, $$d\le\frac{1}{2}.$$
• This method applies to other PDEs too.
  • If the PDE is nonlinear, then linearize it (e.g., Taylor Series), then apply the analysis.
  • If the coefficients are not constant, then pretend they are and apply locally.
    • Take the most conservative $\Delta t$.
    • For example, for $\alpha\Delta t/\Delta x^2\le 1/2$, we have $\Delta t\le (1/2)(\Delta x^2/\alpha)$. Evaluate at $\min(\Delta x^2/\alpha)_i$ on the domain.