FNC Sections 11.1 to 11.3

Learning objectives

• Document a numerical instability when solving diffusion equations
• Work out a solution to this instability using the method of lines
• Work out the notion of absolute stability

Diffusion equations

• Focus is on diffusion equations (seems to be also referred to as diffusive processes, parabolic PDEs, evolutionary PDEs)
• In particular, finding a solution $u(x,t)$ to the heat equation $$u_t=k{u}_{xx},$$ where $k$ is the constant diffusion coefficient.
• $k$ usually set to one in the exposition
• Black-Scholes equation is the motivating example

Source of instability

• Solution involves a function of two variables
• Natural to apply material in Chapter 6 and discretize over space $x$ and time $t$
  • Observe the replacement of the derivatives by finite differences.
  • Left hand side uses a forward finite difference, while right hand side uses a centered difference.
  • Section 11.1 Problem 3 asks the reader to explore why.
• Seemingly works as in Demo 11.1.4, but “trouble lurks just around the corner” in Demo 11.1.5.

New in Julia: Section 11.1

• reshape() used in the demos for “transposing” a vector of strings
• @animate used in creating a movie for the generated plots, mp4() to save as a video file

Dealing with instability

• Separate the discretizations rather than simultaneously discretize!
• Discretize space only for $u_t=u_{xx}$.
  • Define a vector $\mathbf{u}$ by $$\mathbf{u}(t)=\left[\begin{array}{c}{u}_0(t)\\ u_1(t)\\ ⋮\\ u_n(t)\end{array}\right].$$
  • Leads to a linear, constant-coefficient system of ordinary differential equations: $$\frac{d\mathbf{u}(t)}{dt}={\mathbf{D}}_{xx}\mathbf{u}(t)$$
  • ${\mathbf{D}}_{xx}$ seen before in Section 10.3, but adapted to the special case of periodic end conditions (FNC.diffper())
  • $h$ is not time step! Rather it is $h=1/m$, where $m$ is the number of nodes in space discretization.
• Afterwards, just a matter of applying methods in Chapter 6!
• But “trouble lurks just around the corner” because of a missing concept. See Demo 11.2.3.

How to apply Chapter 6 to semi-discretized PDE

• Time steps now in $\tau$ rather than $h$.
• Example 1: Euler IVP integrator
  • Iteration given by $$u_{i+1}=u_i+hf(t_i,u_i),i=0,\dots ,n-1.$$
  • Now becomes vectorized: $${\mathbf{u}}_{j+1}={\mathbf{u}}_j+\tau ({\mathbf{D}}_{xx}{\mathbf{u}}_j)=(\mathbf{I}+\tau {\mathbf{D}}_{xx}){\mathbf{u}}_j,i=0,\dots ,n.$$
  • But unstable, see Demo 11.2.3.
• Example 2: Backward Euler or one-step Adams-Moulton (AM1) IVP integrator
  • Iteration given by $$u_{i+1}=u_i+h{f}_{i+1}$$ with the Chapter 6 convention $f_i=f(t_i,u_i)$
  • Now vectorized as: $$\begin{array}{rl}{\mathbf{u}}_{j+1}& ={\mathbf{u}}_j+\tau ({\mathbf{D}}_{xx}{\mathbf{u}}_{j+1})\\ (\mathbf{I}-\tau {\mathbf{D}}_{xx}){\mathbf{u}}_{j+1}& ={\mathbf{u}}_j.\end{array}$$
  • Require inversion of sparse matrix $\mathbf{I}-\tau {\mathbf{D}}_{xx}$, but fast and stable algorithms available
  • Reinforces the finding in Chapter 6 that implicit time stepping methods are more stable
• Example 3 (Problems 3 and 4 of Section 11.2): two-step Adams-Moulton (AM2) or trapezoid IVP integrator
  • Iteration given by $$u_{i+1}=u_i+h(\frac{1}{2}{f}_{i+1}+\frac{1}{2}{f}_i)$$
  • Now becomes vectorized: $$\begin{array}{rl}{\mathbf{u}}_{j+1}& ={\mathbf{u}}_j+\tau (\frac{1}{2}{\mathbf{D}}_{xx}{\mathbf{u}}_{j+1}+\frac{1}{2}{\mathbf{D}}_{xx}{\mathbf{u}}_j)\\ (\mathbf{I}-0.5\tau {\mathbf{D}}_{xx}){\mathbf{u}}_{j+1}& =(\mathbf{I}+0.5\tau {\mathbf{D}}_{xx}){\mathbf{u}}_j.\end{array}$$
  • Called Crank–Nicolson method
  • Mixes Euler and backward Euler, still stable, refer to Problem 8 of Section 11.3

Absolute stability

• Notation might be unnerving: $\tau$ is the step size, yet it is used as the truncation error in Chapter 6
• All time-stepping methods are zero-stable as $\tau \to 0$.
• But now explore case where $\tau$ stays fixed while the number of nodes $n$ goes to infinity.
• The choice of $\tau$ is crucial to control the growth of the solutions to differential equations. Want boundedness as $n\to \mathrm{\infty }$.
• Theoretical analysis in this section:
  • Decouple system of ODEs.
  • Look at complex eigenvalues $\lambda$’s produced by the decoupling.
  • Key to formulating absolute stability: $\xi =\tau \ast \lambda$ is the focus
• Prefer methods with no time step restrictions