Lecture 12: Numerical Methods for Partial Differential Equations

Topics

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

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

Numerical Examples

Introduction

In finance, most partial differential equations (PDE) we encounter are of the form

$\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}} \large{ 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 }$

Mostly referred to as the Wave equation

Parabolic

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

$\;$
The canonical form:

$\large{ \PDut = d \;\PDuxx }$

Mostly referred to as the Heat equation

Elliptic

$b^2(x,t) - a(x,t)c(x,t) < 0$

$\;$
The canonical form:

$\large{ \PDutt + \PDuxx = 0 }$

Mostly referred to as the Laplace or Poisson equation

Initial and Boundary Conditions

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

$\large{ x(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 $c = 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}$

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.

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$

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

Approximate the Black-Scholes Eqn by

$\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} 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

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 }$

Which is the CFL condition for the explicit method. This says for the explicit method to be stable, the time step needs to be pretty small.
$r\neq 0$ case?

The Implicit Method

Approximate the Black-Scholes Eqn by

$\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} 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} \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)$