FNC Sections 11.1 to 11.3

Author

Andrew Pua

Published

January 31, 2024

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 ut=kuxx, 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 ut=uxx.
    • Define a vector u by u(t)=[u0(t)u1(t)un(t)].
    • Leads to a linear, constant-coefficient system of ordinary differential equations: du(t)dt=Dxxu(t)
    • Dxx 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 τ rather than h.
  • Example 1: Euler IVP integrator
    • Iteration given by ui+1=ui+hf(ti,ui),i=0,,n1.
    • Now becomes vectorized: uj+1=uj+τ(Dxxuj)=(I+τDxx)uj,i=0,,n.
    • But unstable, see Demo 11.2.3.
  • Example 2: Backward Euler or one-step Adams-Moulton (AM1) IVP integrator
    • Iteration given by ui+1=ui+hfi+1 with the Chapter 6 convention fi=f(ti,ui)
    • Now vectorized as: uj+1=uj+τ(Dxxuj+1)(IτDxx)uj+1=uj.
    • Require inversion of sparse matrix IτDxx, 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 ui+1=ui+h(12fi+1+12fi)
    • Now becomes vectorized: uj+1=uj+τ(12Dxxuj+1+12Dxxuj)(I0.5τDxx)uj+1=(I+0.5τDxx)uj.
    • 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 τ0.
  • But now explore case where τ stays fixed while the number of nodes n goes to infinity.
  • The choice of τ is crucial to control the growth of the solutions to differential equations. Want boundedness as n.
  • 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