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=ku_{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) = \begin{bmatrix} u_0(t) \\ u_1(t) \\ \vdots \\ u_n(t) \end{bmatrix}.\]
- Leads to a linear, constant-coefficient system of ordinary differential equations: \[\frac{d \mathbf{u}(t)}{d t} = \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 + h f(t_i,u_i), \qquad i=0,\ldots,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, \qquad i=0,\ldots,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 + hf_{i+1}\] with the Chapter 6 convention \(f_i=f(t_i,u_i)\)
- Now vectorized as: \[\begin{split}
\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{split}\]
- 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\left(\frac{1}{2}f_{i+1}+\frac{1}{2}f_i\right)\]
- Now becomes vectorized: \[\begin{split}
\mathbf{u}_{j+1} &= \mathbf{u}_j + \tau \left(\frac{1}{2}\mathbf{D}_{xx} \mathbf{u}_{j+1}+\frac{1}{2}\mathbf{D}_{xx} \mathbf{u}_{j}\right)\\
(\mathbf{I} - 0.5\tau \mathbf{D}_{xx}) \mathbf{u}_{j+1} &= (\mathbf{I} + 0.5\tau \mathbf{D}_{xx})\mathbf{u}_j.
\end{split}\]
- 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\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 * \lambda\) is the focus
- Prefer methods with no time step restrictions