Lecture 12: Numerical Methods for Partial Differential Equations

Topics

• I. Introduction
  • Second Order PDE Classification
  • Initial and Boundary Conditions
  • Black-Scholes Equation

• II. Finite Difference Methods for Solving PDEs

  • Discretization
  • Finite Difference Methods
  • The Explicit Methods
  • Consistency, Stability and Convergence
  • The Implicit Methods
  • The Crank-Nicholson Method

• III. Numerical Examples

I. Introduction

• In finance, most partial differential equations (PDE) we encounter are of the second order (the highest derivative)

$$ \renewcommand{PDut}{\frac{\partial u}{\partial t}} \renewcommand{PDux}{\frac{\partial u}{\partial x}} \renewcommand{PDutt}{\frac{\partial ^2u}{\partial t^2}} \renewcommand{PDuxx}{\frac{\partial ^2u}{\partial x^2}} \renewcommand{PDuyy}{\frac{\partial ^2u}{\partial y^2}} a(x,t)\PDutt + 2b(x,t)\frac{\partial ^2u}{\partial t\partial x} + c(x,t)\PDuxx + d(x,t)\PDut + e(x,t)\PDux = f(t,x,u) $$

• It can be classified into three categories: Hyperbolic, Parabolic and Elliptic

• In general, the different types of equations make a big difference in how the solutions behave and

• on how we can solve them more effectively

Classification

Hyperbolic

• $b^2(x,t) - a(x,t)c(x,t) > 0 $ $\;$

• The canonical form: $$ \large{ \PDutt = c^2 \PDuxx } $$

• No steady state and no diffusion, mostly referred to as the Wave equation or Advection equation

Parabolic

• $b^2(x,t) - a(x,t)c(x,t) = 0 $ $\;$

• The canonical form: $$ \large{ \PDut = d \;\PDuxx } $$

• Evolving to a steady state due to diffusion, referred to as the Heat equation

Elliptic

• $b^2(x,y) - a(x,y)c(x,y) < 0 $ $\;$

• The canonical form: $$ \large{ \PDuxx + \PDuyy = f(x,y) } $$

• Steady state (time independent), mostly referred to as the Laplace ($f(x,y) = 0$) or Poisson equation.

Initial and Boundary Conditions

• Depending on the problem, we usually have initial conditions (IC)

$$ \large{ u(x, t=0) = u_0(x) } $$

• and/or boundary conditions (BC)

$$ \large{ a \; u + b\; \PDux = c, \;\forall \; x \in \partial Q, \;\forall \; t } $$

• and the BC is called Dirichlet if $b = 0$, Neumann if $a = 0$, or Robin if $a\cdot b \neq 0$ .

Black-Scholes Equation

$$ \renewcommand{PDuS}{\frac{\partial u}{\partial S}} \renewcommand{PDuSS}{\frac{\partial ^2u}{\partial S^2}} \PDut + \frac{1}{2}\sigma^2S^2\PDuSS + rS\PDuS - ru = 0 $$

• It belongs to the parabolic type, and can be converted into the standard Heat equation form by a change-of-variables

• For European Call option, the "initial" (the payoff function) and boundary conditions are

\begin{aligned} &u(S,T) = \max(S-K, 0); \\ &u(0,t) = 0; \\ &u(S,t) = S - Ke^{-r(T-t)}, \;\mbox{as } \; S \rightarrow \infty \end{aligned}

• For European Put option, the "initial" (the payoff function) and boundary conditions are

\begin{aligned} \renewcommand{FDut}{\frac{u_{i,k+1}-u_{i,k}}{\triangle t}} \renewcommand{FDutb}{\frac{u_{i,k}-u_{i,k-1}}{\triangle t}} \renewcommand{FDutc}{\frac{u_{i,k+1}-u_{i,k-1}}{2\triangle t}} \renewcommand{FDutt}{\frac{u_{i,k+1}-2u_{i,k}+u_{i,k-1}}{\triangle t^2}} \renewcommand{FDux}{\frac{u_{i+1,k}-u_{i,k}}{\triangle x}} \renewcommand{FDuxb}{\frac{u_{i,k}-u_{i-1,k}}{\triangle x}} \renewcommand{FDuxc}{\frac{u_{i+1,k}-u_{i-1,k}}{2\triangle x}} \renewcommand{FDuxx}{\frac{u_{i+1,k}-2u_{i,k}+u_{i-1,k}}{\triangle x^2}} &u(S,T) = \max(K-S, 0); \\ &u(0,t) = Ke^{-r(T-t)}; \\ &u(S,t) = 0, \;\mbox{as } \; S \rightarrow \infty \end{aligned}

II. Finite Difference Methods for Solving PDEs

• Very rare that the PDEs have closed form solutions, in general, can only be solved numerically.

• Will focus on the finite difference method (FDM): other numerical methods exist and can be more appropriate---very much depends on the actual problems at hand.

Discretization

• First step in solving PDEs using FDM is to represent the solution $u(x,t)$ as a discrete collection of values at a well distributed grid points in space and time in the proper domain.

• For example, for the rectangular domain in space and time with $t \in [0, T]$ and $x\in [x_{min}, x_{max}]$, we can represent the points simply as

$$t_0 = 0, t_1, t_2, \cdots, t_M = T$$

in time and

$$x_0 = x_{min}, x_1, x_2, \cdots, x_N = x_{max}$$

in space. The grid is uniform, i.e. $t_{k+1} = t_k + \triangle t, \;\triangle t = \frac{T}{M}$ and $x_{i+1} = x_i + \triangle x, \;\triangle x = \frac{x_{max}-x_{min}}{N}$.

• The discretized version of the solution will be represented by the values of

$$ u_{i,k} = u(x_i, t_k), \;\mbox{for } i = 0,1,\cdots,N, \mbox{ and } k = 0,1,\cdots,M $$

• The values of $u(x,t)$ at any other desired points can be approximated via interpolation.

• A uniform mesh may not necessary be the most efficient form to work with, in fact, it rarely is.

• A greatly simplified rule-of-thumb: the mesh needs to be refined around the region where the function varies a lot

• On the other hand, the mesh can be relatively coarse if the function is smooth and changing slowly

• In finance, the time grid is well advised to keep the important dates as grid points: such as cash flow dates, contractual schedule dates, etc. The spacing of these dates is most likely not uniform.

Finite Difference Methods(FDM)

• Toolkit for finite difference approximation

Partial Derivative	Finite Difference	Type	Order
$\PDux$	$\FDux$	forward	1st in $x$
$\PDux$	$\FDuxb$	backward	1st in $x$
$\PDux$	$\FDuxc$	central	2nd in $x$
$\PDuxx$	$\FDuxx$	symmetric	2nd in $x$
$\PDut$	$\FDut$	forward	1st in $t$
$\PDut$	$\FDutb$	backward	1st in $t$
$\PDut$	$\FDutc$	central	2nd in $t$
$\PDutt$	$\FDutt$	symmetric	2nd in $t$

The Explicit Method

• For convenience, reverse the time for Black-Scholes notation: $ t \leftarrow T - t$

$$ \small \PDut - \frac{1}{2}\sigma^2S^2\PDuSS - rS\PDuS + ru = 0 $$

• Approximate the Black-Scholes Eqn by

$$ \small \FDut - \frac{1}{2} \sigma^2 x_i^2 \FDuxx - r x_i \FDuxc + r u_{i,k} = 0 $$

• This can be rewritten simply as

\begin{aligned} \small u_{i,k+1} & = \left\{ \frac{1}{2}\sigma^2 x_i^2 \frac{\triangle t}{\triangle x^2} - r x_i \frac{\triangle t}{2\triangle x} \right\} u_{i-1,k} \\ & + \left\{ 1 - \sigma^2 x_i^2 \frac{\triangle t}{\triangle x^2} -r \triangle t \right\} u_{i,k} \\ & + \left\{ \frac{1}{2}\sigma^2 x_i^2 \frac{\triangle t}{\triangle x^2} + r x_i \frac{\triangle t}{2\triangle x} \right\} u_{i+1,k} \end{aligned}

Consistency, Stability and Convergence

• For any FDM that are used to solve practical problems, we should ask

  • 1) How acurate is the method?

  • 2) Does it converge?

  • 3) What is the best choice of step sizes?

Consistency

• For the explicit method above, the local truncation error is

$$ T(x,t) = O(\triangle x^2 + \triangle t) $$

• Taylor series expansion can verify this (or one can simply read this off of the Toolkit Table).

• So the explicit method above is consistent of order 1 in time and order 2 in spatial variable.

Von Neumann Stability Analysis

• The Von Neumann method of analyzing the stability of FDM boils down to substitute $u_{j,k}$ with

$$ \epsilon_{j,k} = e^{a t_k}e^{il_mx_j} $$

• The logic behind this: assuming $v_{j,k}$ is the solution of the FDM that did not suffer from finite precision arithmetic and is without the round off error, then the error

$$ \epsilon_{j,k} = v_{j,k} - u_{j,k} $$

$\hspace{0.2in}$ again satifies the FDM due to linearity.

• Looking for $\epsilon_{j,k} $ in the Fourier expansion representation,

$$ \epsilon(x,t) = \sum_{m= 1}^{M} e^{at}e^{il_mx} $$

where $e^{at}$ is a special form of the amplitude, and $l_m$ is the wavelength: $l_m = \frac{\pi m}{L}, m = 1, \cdots, M, \mbox{ and } M = \frac{L}{\triangle x}$

• Due to linearity again, the error analysis can be performed (without any loss of generality) by looking at a single term of the expansion

$$ e^{at_k}e^{il_mx_j} $$

• In fact, we are really interested in the amplification of the error,

$$ \frac{\epsilon_{j,k+1}}{\epsilon_{j,k}} = e^{a\triangle t} $$

• The stability analysis in the Von Neumann method is simply to analyze the form of $e^{a\triangle t}$ so that

$$ ||e^{a\triangle t}|| < 1 $$

• For the explicit method,

\begin{aligned} e^{a (t_k + \triangle t)}e^{il_mx_j} & = \left\{ \frac{1}{2}\sigma^2 x_j^2 \frac{\triangle t}{\triangle x^2} - r x_j \frac{\triangle t}{2\triangle x} \right\} e^{a t_k}e^{il_mx_{j-1}} \\ & + \left\{ 1 - \sigma^2 x_j^2 \frac{\triangle t}{\triangle x^2} -r \triangle t \right\} e^{a t_k}e^{il_mx_j} \\ & + \left\{ \frac{1}{2}\sigma^2 x_j^2 \frac{\triangle t}{\triangle x^2} + r x_j \frac{\triangle t}{2\triangle x} \right\} e^{a t_k}e^{il_mx_{j+1}} \end{aligned}

• Simplify it into,

\begin{aligned} e^{a\triangle t} & = \left\{ \frac{1}{2}\sigma^2 x_j^2 \frac{\triangle t}{\triangle x^2} - r x_j \frac{\triangle t}{2\triangle x} \right\} e^{-il_m\triangle x} \\ & + \left\{ 1 - \sigma^2 x_j^2 \frac{\triangle t}{\triangle x^2} -r \triangle t \right\} \\ & + \left\{ \frac{1}{2}\sigma^2 x_j^2 \frac{\triangle t}{\triangle x^2} + r x_j \frac{\triangle t}{2\triangle x} \right\} e^{il_m\triangle x} \end{aligned}

• Looking at the simple case that $r = 0$,

\begin{aligned} e^{a\triangle t} & = \frac{1}{2}\sigma^2 x_j^2 \frac{\triangle t}{\triangle x^2} e^{-il_m\triangle x} \\ & + 1 - \sigma^2 x_j^2 \frac{\triangle t}{\triangle x^2} \\ & + \frac{1}{2}\sigma^2 x_j^2 \frac{\triangle t}{\triangle x^2} e^{il_m\triangle x} \end{aligned}

• which is

\begin{aligned} e^{a\triangle t} &= \sigma^2 x_j^2 \frac{\triangle t}{\triangle x^2} \cos(l_m\triangle x) + 1 - \sigma^2 x_j^2 \frac{\triangle t}{\triangle x^2} \\ & = 1 - \sigma^2 x_j^2 \frac{\triangle t}{\triangle x^2} \cdot 2\sin^2(l_m\triangle x/2) \end{aligned}

• The stability condition will be saisfied if

$$ || 1 - \sigma^2 x_j^2 \frac{\triangle t}{\triangle x^2} \cdot 2\sin^2(l_m\triangle x/2) || < 1 $$

• Or if

$$ \triangle t < \frac{\triangle x^2}{ \sigma^2 x_{max}^2 } $$

• This is the CFL condition for the explicit method, which says, in order for the explicit method to be stable, the time step size needs to be really small.

• $r\neq 0$ case?

Lax Theorem: Consistency + Stability = Convergence

• The Lax equivalence theorem holds for PDE as well.

• Takeaway: when implementing FDM for PDEs, study the CFL (Courant–Friedrichs–Lewy) stability condition before building your discretization grid.

The Implicit Method

• Approximate the Black-Scholes Eqn by

$$ \small \FDutb - \frac{1}{2} \sigma^2 x_i^2 \FDuxx - r x_i \FDuxc = - r u_{i,k} $$

• Which can be rewritten as

\begin{aligned} \small u_{i,k-1} & = \left\{ -\frac{1}{2}\sigma^2 x_i^2 \frac{\triangle t}{\triangle x^2} + r x_i \frac{\triangle t}{2\triangle x} \right\} u_{i-1,k} \\ & + \left\{ 1 + \sigma^2 x_i^2 \frac{\triangle t}{\triangle x^2} + r \triangle t \right\} u_{i,k} \\ & + \left\{ -\frac{1}{2}\sigma^2 x_i^2 \frac{\triangle t}{\triangle x^2} - r x_i \frac{\triangle t}{2\triangle x} \right\} u_{i+1,k} \end{aligned}

• The local truncation error

$$ T(x,t) = O(\triangle x^2 + \triangle t) $$

• And the amplification factor for $r=0$

$$ e^{a\triangle t} = \frac{1}{1 + \sigma^2 x_j^2 \frac{\triangle t}{\triangle x^2} \cdot 2\sin^2(l_m\triangle x/2)} $$

• So the implicit method is unconditionally stable.

The Crank-Nicholson Method

• Crank-Nicholson is an average of the explicit and implicit scheme:

\begin{aligned} \small \FDut & = \frac{1}{4} \sigma^2 x_i^2 \left\{ \FDuxx + \frac{u_{i+1,k+1}-2u_{i,k+1}+u_{i-1,k+1}}{\triangle x^2} \right\} \\ & + \frac{1}{2} r x_i \left\{ \FDuxc + \frac{u_{i+1,k+1}-u_{i-1,k+1}}{2\triangle x} \right\} \\ & - \frac{1}{2} r \left\{ u_{i,k} + u_{i,k+1} \right\} \end{aligned}

• The local truncation error is (Note: expanding the Taylor series at $(i, k+\frac{1}{2})$)

$$ T(x,t) = O(\triangle x^2 + \triangle t^2) $$

III. Numerical Examples

• Will go through the Black Scholes example in next lecture.

Homework

• Derive the CFL condition for the Crank-Nicholson method for the Black-Scholes PDE, assuming $r=0$.

Reference

• L. Andersen and V. Piterbarg, Chapter 4 of Interest Rate Modeling, Volumn I, Atlantic Financial Press, 2010.