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 to the heat equation where is the constant diffusion coefficient.
- 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 and time
- 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 .
- Define a vector by
- Leads to a linear, constant-coefficient system of ordinary differential equations:
- seen before in Section 10.3, but adapted to the special case of periodic end conditions (
FNC.diffper()
)
- is not time step! Rather it is , where 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 rather than .
- Example 1: Euler IVP integrator
- Iteration given by
- Now becomes vectorized:
- But unstable, see Demo 11.2.3.
- Example 2: Backward Euler or one-step Adams-Moulton (AM1) IVP integrator
- Iteration given by with the Chapter 6 convention
- Now vectorized as:
- Require inversion of sparse matrix , 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
- Now becomes vectorized:
- 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: is the step size, yet it is used as the truncation error in Chapter 6
- All time-stepping methods are zero-stable as .
- But now explore case where stays fixed while the number of nodes goes to infinity.
- The choice of is crucial to control the growth of the solutions to differential equations. Want boundedness as .
- Theoretical analysis in this section:
- Decouple system of ODEs.
- Look at complex eigenvalues ’s produced by the decoupling.
- Key to formulating absolute stability: is the focus
- Prefer methods with no time step restrictions