Stoch Proc - Lab/Tutorial 10 Solutions
Finite Differences solution of PDEs

The theoretical section is based on course notes for “Numerical Methods in Finance” by Eckhard Platen.

Theory

This code solves PDEs of the following form:

((∂u(τ,x))/(∂τ))  =  a(τ,x)((∂u(τ,x))/(∂x)) + b(τ,x)(((∂^2)u(τ,x))/(∂(x^2))) + c(τ,x) u(τ,x)        s.t.        u(0,x) = h(x).

Defining mesh

Discretisation of τ
(Subscript[τ, j]) = j Δτ,           j = 0,1,2,...,m,    Δτ=((Subscript[τ, max])/m)

Discretisation of x
(Subscript[x, k]) = (Subscript[x, min]) + k Δx,    k=0,1,2,...,n,    Δx=(((Subscript[x, max]) - (Subscript[x, min]))/n)

Discretisation of u
u((Subscript[τ, j]),(Subscript[x, k]))  ≡  (Subscript[u, k]^j)
On the mesh define     u((Subscript[τ, j])) := (({{u((Subscript[τ, j]),(Subscript[x, 0]))}, {u((Subscript[τ, j]),(Subscript[x, 1]))}, {u((Subscript[τ, j]),(Subscript[x, 2]))}, {⋮}, {u((Subscript[τ, j]),(Subscript[x, n-2]))}, {u((Subscript[τ, j]),(Subscript[x, n-1]))}, {u((Subscript[τ, j]),(Subscript[x, n]))}})) ≡ (({{(Subscript[u, 0]^j)}, {(Subscript[u, 1]^j)}, {(Subscript[u, 2]^j)}, {⋮}, {(Subscript[u, n-2]^j)}, {(Subscript[u, n-1]^j)}, {(Subscript[u, n]^j)}})) =: (u^j).

We seek solution of PDE as
(({{ (u^0),  (u^1),  (u^2), ⋯,  (u^(m-1)),  (u^m)}})) = (({{(Subscript[u, 0]^0), (Subscript[u, 0]^1), (Subscript[u, 0]^2), ⋯, (Subscript[u, 0]^(m-1)), (Subscript[u, 0]^m)}, {(Subscript[u, 1]^0), (Subscript[u, 1]^1), (Subscript[u, 1]^2), ⋯, (Subscript[u, 1]^(m-1)), (Subscript[u, 1]^m)}, {(Subscript[u, 2]^0), (Subscript[u, 2]^1), (Subscript[u, 2]^2), ⋯, (Subscript[u, 2]^(m-1)), (Subscript[u, 2]^m)}, {⋮, ⋮, ⋮, ⋱, ⋮, ⋮}, {(Subscript[u, n-2]^0), (Subscript[u, n-2]^1), (Subscript[u, n-2]^2), ⋯, (Subscript[u, n-2]^(m-1)), (Subscript[u, n-2]^m)}, {(Subscript[u, n-1]^0), (Subscript[u, n-1]^1), (Subscript[u, n-1]^2), ⋯, (Subscript[u, n-1]^(m-1)), (Subscript[u, n-1]^m)}, {(Subscript[u, n]^0), (Subscript[u, n]^1), (Subscript[u, n]^2), ⋯, (Subscript[u, n]^(m-1)), (Subscript[u, n]^m)}})).

Discretisation of a
a((Subscript[τ, j]),(Subscript[x, k]))  ≡  (Subscript[a, k]^j)

Discretisation of b
b((Subscript[τ, j]),(Subscript[x, k]))  ≡  (Subscript[b, k]^j)

Discretisation of c
c((Subscript[τ, j]),(Subscript[x, k]))  ≡  (Subscript[c, k]^j)

Discretisation of h
h((Subscript[x, k]))  ≡  (Subscript[h, k])
Note that    (u^0) = (({{(Subscript[u, 0]^0)}, {(Subscript[u, 1]^0)}, {(Subscript[u, 2]^0)}, {⋮}, {(Subscript[u, n-2]^0)}, {(Subscript[u, n-1]^0)}, {(Subscript[u, n]^0)}})) = (({{(Subscript[h, 0])}, {(Subscript[h, 1])}, {(Subscript[h, 2])}, {⋮}, {(Subscript[h, n-2])}, {(Subscript[h, n-1])}, {(Subscript[h, n])}}))
so we have
(({{ (u^0),  (u^1),  (u^2), ⋯,  (u^(m-1)),  (u^m)}})) = (({{(Subscript[h, 0]), (Subscript[u, 0]^1), (Subscript[u, 0]^2), ⋯, (Subscript[u, 0]^(m-1)), (Subscript[u, 0]^m)}, {(Subscript[h, 1]), (Subscript[u, 1]^1), (Subscript[u, 1]^2), ⋯, (Subscript[u, 1]^(m-1)), (Subscript[u, 1]^m)}, {(Subscript[h, 2]), (Subscript[u, 2]^1), (Subscript[u, 2]^2), ⋯, (Subscript[u, 2]^(m-1)), (Subscript[u, 2]^m)}, {⋮, ⋮, ⋮, ⋱, ⋮, ⋮}, {(Subscript[h, n-2]), (Subscript[u, n-2]^1), (Subscript[u, n-2]^2), ⋯, (Subscript[u, n-2]^(m-1)), (Subscript[u, n-2]^m)}, {(Subscript[h, n-1]), (Subscript[u, n-1]^1), (Subscript[u, n-1]^2), ⋯, (Subscript[u, n-1]^(m-1)), (Subscript[u, n-1]^m)}, {(Subscript[h, n]), (Subscript[u, n]^1), (Subscript[u, n]^2), ⋯, (Subscript[u, n]^(m-1)), (Subscript[u, n]^m)}})).

Time derivative

We use the backward finite difference
((∂u((Subscript[τ, j]),(Subscript[x, k])))/(∂τ)) ≡ ((∂u(τ,x))/(∂τ))(Subscript[|, τ=(Subscript[τ, j]),x=(Subscript[x, k])]) ≈ (((Subscript[u, k]^j) - (Subscript[u, k]^(j-1)))/(△τ)),    j = 1,2,...,m.
On the mesh we have    ((∂u((Subscript[τ, j])))/(∂τ)) = (({{((∂u((Subscript[τ, j]),(Subscript[x, 0])))/(∂τ))}, {((∂u((Subscript[τ, j]),(Subscript[x, 1])))/(∂τ))}, {((∂u((Subscript[τ, j]),(Subscript[x, 2])))/(∂τ))}, {⋮}, {((∂u((Subscript[τ, j]),(Subscript[x, n-2])))/(∂τ))}, {((∂u((Subscript[τ, j]),(Subscript[x, n-1])))/(∂τ))}, {((∂u((Subscript[τ, j]),(Subscript[x, n])))/(∂τ))}})) ≈ (({{(((Subscript[u, 0]^j) - (Subscript[u, 0]^(j-1)))/(△τ))}, {(((Subscript[u, 1]^j) - (Subscript[u, 1]^(j-1)))/(△τ))}, {(((Subscript[u, 2]^j) - (Subscript[u, 2]^(j-1)))/(△τ))}, {⋮}, {(((Subscript[u, n-2]^j) - (Subscript[u, n-2]^(j-1)))/(△τ))}, {(((Subscript[u, n-1]^j) - (Subscript[u, n-1]^(j-1)))/(△τ))}, {(((Subscript[u, n]^j) - (Subscript[u, n]^(j-1)))/(△τ))}})) =: (((u^j) - (u^(j-1)))/(△τ)).

Spatial derivatives

For the first order derivative we use the central finite difference
((∂u((Subscript[τ, j]),(Subscript[x, k])))/(∂x)) ≡ ((∂u(τ,x))/(∂x))(Subscript[|, τ=(Subscript[τ, j]),x=(Subscript[x, k])]) ≈ (((Subscript[u, k+1]^j) - (Subscript[u, k-1]^j))/(2△x)),    k=1,2,...,n-1.

For the second order derivative we use the second order central finite difference
(((∂^2)u((Subscript[τ, j]),(Subscript[x, k])))/(∂(x^2))) ≡ (((∂^2)u(τ,x))/(∂(x^2)))(Subscript[|, τ=(Subscript[τ, j]),x=(Subscript[x, k])]) ≈  (((((Subscript[u, k+1]^j) - (Subscript[u, k]^j))/(△x))-(((Subscript[u, k]^j) - (Subscript[u, k-1]^j))/(△x)))/(△x)) = (((Subscript[u, k+1]^j)-2(Subscript[u, k]^j)+(Subscript[u, k-1]^j))/((((△x)^2)))),    k=1,2,...,n-1.

Then for τ=(Subscript[τ, j]) and x=(Subscript[x, k])
a((Subscript[τ, j]),(Subscript[x, k]))((∂u((Subscript[τ, j]),(Subscript[x, k])))/(∂x)) + b((Subscript[τ, j]),(Subscript[x, k]))(((∂^2)u((Subscript[τ, j]),(Subscript[x, k])))/(∂(x^2))) + c((Subscript[τ, j]),(Subscript[x, k])) u((Subscript[τ, j]),(Subscript[x, k]))
  ≈ (Subscript[a, k]^j)(((Subscript[u, k+1]^j)-(Subscript[u, k-1]^j))/(2△x))  + (Subscript[b, k]^j)(((Subscript[u, k+1]^j)-2(Subscript[u, k]^j)+(Subscript[u, k-1]^j))/((((△x)^2))))  + (Subscript[c, k]^j)(Subscript[u, k]^j)
  = (((Subscript[b, k]^j)/((((△x)^2))))-((Subscript[a, k]^j)/(2△x)))(Subscript[u, k-1]^j) + ((Subscript[c, k]^j)-((2(Subscript[b, k]^j))/((((△x)^2)))))(Subscript[u, k]^j) + (((Subscript[b, k]^j)/((((△x)^2))))+((Subscript[a, k]^j)/(2△x)))(Subscript[u, k+1]^j)
  = (Subscript[A, k,k-1]^j)(Subscript[u, k-1]^j) + (Subscript[A, k,k]^j)(Subscript[u, k]^j) +(Subscript[A, k,k+1]^j)(Subscript[u, k+1]^j)  =  (({{(Subscript[A, k,k-1]^j), (Subscript[A, k,k]^j), (Subscript[A, k,k+1]^j)}}))(({{(Subscript[u, k-1]^j)}, {(Subscript[u, k]^j)}, {(Subscript[u, k+1]^j)}}))
where    (Subscript[A, k,k-1]^j) = ((Subscript[b, k]^j)/((((△x)^2))))-((Subscript[a, k]^j)/(2△x)),    (Subscript[A, k,k]^j) = (Subscript[c, k]^j)- ((2(Subscript[b, k]^j))/((((△x)^2))))     and    (Subscript[A, k,k+1]^j) = ((Subscript[b, k]^j)/((((△x)^2))))+ ((Subscript[a, k]^j)/(2△x)).

Therefore on the mesh

A((Subscript[τ, j]))u((Subscript[τ, j])):= (({{a((Subscript[τ, j]),(Subscript[x, 0]))((∂u((Subscript[τ, j]),(Subscript[x, 0])))/(∂x)) + b((Subscript[τ, j]),(Subscript[x, 0]))((∂u((Subscript[τ, j]),(Subscript[x, 0])))/(∂(x^2))) + c((Subscript[τ, j]),(Subscript[x, 0])) u((Subscript[τ, j]),(Subscript[x, 0]))}, {a((Subscript[τ, j]),(Subscript[x, 1]))((∂u((Subscript[τ, j]),(Subscript[x, 1])))/(∂x)) + b((Subscript[τ, j]),(Subscript[x, 1]))((∂u((Subscript[τ, j]),(Subscript[x, 1])))/(∂(x^2))) + c((Subscript[τ, j]),(Subscript[x, 1])) u((Subscript[τ, j]),(Subscript[x, 1]))}, {a((Subscript[τ, j]),(Subscript[x, 2]))((∂u((Subscript[τ, j]),(Subscript[x, 2])))/(∂x)) + b((Subscript[τ, j]),(Subscript[x, 2]))((∂u((Subscript[τ, j]),(Subscript[x, 2])))/(∂(x^2))) + c((Subscript[τ, j]),(Subscript[x, 2])) u((Subscript[τ, j]),(Subscript[x, 2]))}, {⋮}, {a((Subscript[τ, j]),(Subscript[x, n-2]))((∂u((Subscript[τ, j]),(Subscript[x, n-2])))/(∂x)) + b((Subscript[τ, j]),(Subscript[x, n-2]))((∂u((Subscript[τ, j]),(Subscript[x, n-2])))/(∂(x^2))) + c((Subscript[τ, j]),(Subscript[x, n-2])) u((Subscript[τ, j]),(Subscript[x, n-2]))}, {a((Subscript[τ, j]),(Subscript[x, n-1]))((∂u((Subscript[τ, j]),(Subscript[x, n-1])))/(∂x)) + b((Subscript[τ, j]),(Subscript[x, n-1]))((∂u((Subscript[τ, j]),(Subscript[x, n-1])))/(∂(x^2))) + c((Subscript[τ, j]),(Subscript[x, n-1])) u((Subscript[τ, j]),(Subscript[x, n-1]))}, {a((Subscript[τ, j]),(Subscript[x, n]))((∂u((Subscript[τ, j]),(Subscript[x, n])))/(∂x)) + b((Subscript[τ, j]),(Subscript[x, n]))((∂u((Subscript[τ, j]),(Subscript[x, n])))/(∂(x^2))) + c((Subscript[τ, j]),(Subscript[x, n])) u((Subscript[τ, j]),(Subscript[x, n]))}})) ≈ (({{(Subscript[u, 0]^j)-2(Subscript[u, 1]^j)+(Subscript[u, 2]^j)}, {(Subscript[A, 1,0]^j)(Subscript[u, 0]^j)+(Subscript[A, 1,1]^j)(Subscript[u, 1]^j)+(Subscript[A, 1,2]^j)(Subscript[u, 2]^j)}, {(Subscript[A, 2,1]^j)(Subscript[u, 1]^j)+(Subscript[A, 2,2]^j)(Subscript[u, 2]^j)+(Subscript[A, 2,3]^j)(Subscript[u, 3]^j)}, {⋮}, {(Subscript[A, n-2,n-3]^j)(Subscript[u, n-3]^j) + (Subscript[A, n-2,n-2]^j)(Subscript[u, n-2]^j) +(Subscript[A, n-2,n-1]^j)(Subscript[u, n-1]^j)}, {(Subscript[A, n-1,n-2]^j)(Subscript[u, n-2]^j) + (Subscript[A, n-1,n-1]^j)(Subscript[u, n-1]^j) +(Subscript[A, n-1,n]^j)(Subscript[u, n]^j)}, {(Subscript[u, n-2]^j)-2(Subscript[u, n-1]^j)+(Subscript[u, n]^j)}})) = (A^j)(u^j)
where    (A^j) := (({{1, -2, 1, 0, ..., 0, 0, 0, 0}, {(Subscript[A, 1,0]^j), (Subscript[A, 1,1]^j), (Subscript[A, 1,2]^j), 0, ..., 0, 0, 0, 0}, {0, (Subscript[A, 2,1]^j), (Subscript[A, 2,2]^j), (Subscript[A, 2,3]^j), ..., 0, 0, 0, 0}, {..., ..., ..., ..., ..., ..., ..., ..., ...}, {0, 0, 0, 0, ..., (Subscript[A, n-2,n-3]^j), (Subscript[A, n-2,n-2]^j), (Subscript[A, n-2,n-1]^j), 0}, {0, 0, 0, 0, ..., 0, (Subscript[A, n-1,n-2]^j), (Subscript[A, n-1,n-1]^j), (Subscript[A, n-1,n]^j)}, {0, 0, 0, 0, ..., 0, 1, -2, 1}}))    and    (u^j) = ( ({{(Subscript[u, 0]^j)}, {(Subscript[u, 1]^j)}, {(Subscript[u, 2]^j)}, {⋮}, {(Subscript[u, n-2]^j)}, {(Subscript[u, n-1]^j)}, {(Subscript[u, n]^j)}})).

Note the treatment on the boundaries of the mesh where we use linear extrapolation.

PDE approximation (implicit scheme setup)

On the mesh we have the PDE problem
((∂u((Subscript[τ, j])))/(∂τ))  =  A((Subscript[τ, j]))u((Subscript[τ, j]))
which is
(({{((∂u((Subscript[τ, j]),(Subscript[x, 0])))/(∂τ))}, {((∂u((Subscript[τ, j]),(Subscript[x, 1])))/(∂τ))}, {((∂u((Subscript[τ, j]),(Subscript[x, 2])))/(∂τ))}, {⋮}, {((∂u((Subscript[τ, j]),(Subscript[x, n-2])))/(∂τ))}, {((∂u((Subscript[τ, j]),(Subscript[x, n-1])))/(∂τ))}, {((∂u((Subscript[τ, j]),(Subscript[x, n])))/(∂τ))}})) = (({{a((Subscript[τ, j]),(Subscript[x, 0]))((∂u((Subscript[τ, j]),(Subscript[x, 0])))/(∂x)) + b((Subscript[τ, j]),(Subscript[x, 0]))((∂u((Subscript[τ, j]),(Subscript[x, 0])))/(∂(x^2))) + c((Subscript[τ, j]),(Subscript[x, 0])) u((Subscript[τ, j]),(Subscript[x, 0]))}, {a((Subscript[τ, j]),(Subscript[x, 1]))((∂u((Subscript[τ, j]),(Subscript[x, 1])))/(∂x)) + b((Subscript[τ, j]),(Subscript[x, 1]))((∂u((Subscript[τ, j]),(Subscript[x, 1])))/(∂(x^2))) + c((Subscript[τ, j]),(Subscript[x, 1])) u((Subscript[τ, j]),(Subscript[x, 1]))}, {a((Subscript[τ, j]),(Subscript[x, 2]))((∂u((Subscript[τ, j]),(Subscript[x, 2])))/(∂x)) + b((Subscript[τ, j]),(Subscript[x, 2]))((∂u((Subscript[τ, j]),(Subscript[x, 2])))/(∂(x^2))) + c((Subscript[τ, j]),(Subscript[x, 2])) u((Subscript[τ, j]),(Subscript[x, 2]))}, {⋮}, {a((Subscript[τ, j]),(Subscript[x, n-2]))((∂u((Subscript[τ, j]),(Subscript[x, n-2])))/(∂x)) + b((Subscript[τ, j]),(Subscript[x, n-2]))((∂u((Subscript[τ, j]),(Subscript[x, n-2])))/(∂(x^2))) + c((Subscript[τ, j]),(Subscript[x, n-2])) u((Subscript[τ, j]),(Subscript[x, n-2]))}, {a((Subscript[τ, j]),(Subscript[x, n-1]))((∂u((Subscript[τ, j]),(Subscript[x, n-1])))/(∂x)) + b((Subscript[τ, j]),(Subscript[x, n-1]))((∂u((Subscript[τ, j]),(Subscript[x, n-1])))/(∂(x^2))) + c((Subscript[τ, j]),(Subscript[x, n-1])) u((Subscript[τ, j]),(Subscript[x, n-1]))}, {a((Subscript[τ, j]),(Subscript[x, n]))((∂u((Subscript[τ, j]),(Subscript[x, n])))/(∂x)) + b((Subscript[τ, j]),(Subscript[x, n]))((∂u((Subscript[τ, j]),(Subscript[x, n])))/(∂(x^2))) + c((Subscript[τ, j]),(Subscript[x, n])) u((Subscript[τ, j]),(Subscript[x, n]))}})).

Using the spatial finite difference approximations, this can be approximated as the ODE problem
((du((Subscript[τ, j])))/(dτ)) ≈(A^j)(u^j)
which is
(({{((∂u((Subscript[τ, j]),(Subscript[x, 0])))/(∂τ))}, {((∂u((Subscript[τ, j]),(Subscript[x, 1])))/(∂τ))}, {((∂u((Subscript[τ, j]),(Subscript[x, 2])))/(∂τ))}, {⋮}, {((∂u((Subscript[τ, j]),(Subscript[x, n-2])))/(∂τ))}, {((∂u((Subscript[τ, j]),(Subscript[x, n-1])))/(∂τ))}, {((∂u((Subscript[τ, j]),(Subscript[x, n])))/(∂τ))}})) ≈ (({{(Subscript[u, 0]^j)-2(Subscript[u, 1]^j)+(Subscript[u, 2]^j)}, {(Subscript[A, 1,0]^j)(Subscript[u, 0]^j)+(Subscript[A, 1,1]^j)(Subscript[u, 1]^j)+(Subscript[A, 1,2]^j)(Subscript[u, 2]^j)}, {(Subscript[A, 2,1]^j)(Subscript[u, 1]^j)+(Subscript[A, 2,2]^j)(Subscript[u, 2]^j)+(Subscript[A, 2,3]^j)(Subscript[u, 3]^j)}, {⋮}, {(Subscript[A, n-2,n-3]^j)(Subscript[u, n-3]^j) + (Subscript[A, n-2,n-2]^j)(Subscript[u, n-2]^j) +(Subscript[A, n-2,n-1]^j)(Subscript[u, n-1]^j)}, {(Subscript[A, n-1,n-2]^j)(Subscript[u, n-2]^j) + (Subscript[A, n-1,n-1]^j)(Subscript[u, n-1]^j) +(Subscript[A, n-1,n]^j)(Subscript[u, n]^j)}, {(Subscript[u, n-2]^j)-2(Subscript[u, n-1]^j)+(Subscript[u, n]^j)}})).

Finally,replacing the time derative with finite difference approximation gives
(((u^j) - (u^(j-1)))/(△τ)) ≈ (A^j)(u^j)
which is
(({{(((Subscript[u, 0]^j) - (Subscript[u, 0]^(j-1)))/(△τ))}, {(((Subscript[u, 1]^j) - (Subscript[u, 1]^(j-1)))/(△τ))}, {(((Subscript[u, 2]^j) - (Subscript[u, 2]^(j-1)))/(△τ))}, {⋮}, {(((Subscript[u, n-2]^j) - (Subscript[u, n-3]^(j-1)))/(△τ))}, {(((Subscript[u, n-1]^j) - (Subscript[u, n-2]^(j-1)))/(△τ))}, {(((Subscript[u, n]^j) - (Subscript[u, n-1]^(j-1)))/(△τ))}}))  ≈ (({{(Subscript[u, 0]^j)-2(Subscript[u, 1]^j)+(Subscript[u, 2]^j)}, {(Subscript[A, 1,0]^j)(Subscript[u, 0]^j)+(Subscript[A, 1,1]^j)(Subscript[u, 1]^j)+(Subscript[A, 1,2]^j)(Subscript[u, 2]^j)}, {(Subscript[A, 2,1]^j)(Subscript[u, 1]^j)+(Subscript[A, 2,2]^j)(Subscript[u, 2]^j)+(Subscript[A, 2,3]^j)(Subscript[u, 3]^j)}, {⋮}, {(Subscript[A, n-2,n-3]^j)(Subscript[u, n-3]^j) + (Subscript[A, n-2,n-2]^j)(Subscript[u, n-2]^j) +(Subscript[A, n-2,n-1]^j)(Subscript[u, n-1]^j)}, {(Subscript[A, n-1,n-2]^j)(Subscript[u, n-2]^j) + (Subscript[A, n-1,n-1]^j)(Subscript[u, n-1]^j) +(Subscript[A, n-1,n]^j)(Subscript[u, n]^j)}, {(Subscript[u, n-2]^j)-2(Subscript[u, n-1]^j)+(Subscript[u, n]^j)}})).

So for the implicit scheme we need to solve
(((u^j) - (u^(j-1)))/(△τ)) = (A^j)(u^j),    j∈{1,2,…,m},
which is
(I - (A^j) Δτ)(u^j) = (u^(j-1)),    j∈{1,2,…,m}.

A system of linear equations needs to be solved at each time step (Subscript[τ, j]) for    j∈{1,2,…,m}.

Theta Method

The theta method is a generalisation of the above the encompases implicit, explicit and mixture setups.

Let    0≤θ≤1    so that
(((u^j) - (u^(j-1)))/(△τ))  =  (1-θ)(A^(j-1))(u^(j-1)) + θ (A^j)(u^j),     j∈{1,2,…,m},
which is
(I - θ (A^j) Δτ)(u^j) = (I + (1 - θ) (A^(j-1))Δτ)(u^(j-1)),     j∈{1,2,…,m},
or
(X^j)(u^j) = (Y^(j-1))(u^(j-1)),     j∈{1,2,…,m},
where
(X^j) = I - θ  (A^j)Δτ    and    (Y^(j-1)) = I + (1 - θ) (A^(j-1))Δτ.

The linear system must be solved at each  step (Subscript[τ, j]) for    j∈{1,2,…,m}.

Explicit method

For the "explicit" FD method take    θ=0.

This gives
(u^j) = (I + (A^(j-1))Δτ)(u^(j-1)),    j∈{1,2,…,m}.

Notice at each time step there is no linear system to solve.
However, this scheme is numerically unstable (as we shall see).

Implicit method

For the "implicit" FD method take    θ=1.

This gives
(I - (A^j) Δτ)(u^j) = (u^(j-1)),    j∈{1,2,…,m}.

In this case there is a linear systemn to solve at each step.
This scheme is numerically stable.

Crank-Nicholson scheme

For the "Crank-Nicholson" FD method take    θ=1/2.

This gives
(I - (1/2) (A^j) Δτ)(u^j) = (I + (1/2) (A^(j-1))Δτ)(u^(j-1)),    j∈{1,2,…,m}.

In this case there is a linear systemn to solve at each step.
This scheme is numerically stable and has good convergence properties.

Q1

European option call option

The Black-Scholes PDE problem is
((∂V(t,x))/(∂t)) + (r-q) x ((∂V(t,x))/(∂x)) + (1/2)(σ^2)(x^2)(((∂^2)V(t,x))/(∂(x^2))) = r V(t,x)
s.t.    V(T,x) = max(x-K,0).

Note this is consistent with the probabilitistic representation
 V(T,x) = (e^(-r (T-T)))𝔼[max((Subscript[S, T])-K,0)|(Subscript[S, T])=x] = 𝔼[max(x-K,0)|(Subscript[S, T])=x]
  = max(x-K,0).

Let    τ(t)=T-t    and    u(τ,x)≡u(τ(t),x)=V(t,x).

Then
((∂V(t,x))/(∂t)) = ((∂u(τ,x))/(∂t)) = ((∂u(τ,x))/(∂τ))((∂τ)/(∂t)) = -((∂u(τ,x))/(∂τ)),
((∂V(t,x))/(∂x)) = ((∂u(τ,x))/(∂x))
and
(((∂^2)V(t,x))/(∂(x^2))) = (((∂^2) u(τ,x))/(∂(x^2)))
with
V(T,x) = u(τ(T),x) = u(0,x).

Therefore transformed PDE problem is
-((∂u(τ,x))/(∂τ)) + (r-q)x ((∂u(τ,x))/(∂x)) + (1/2)(σ^2)(x^2)(((∂^2)u(τ,x))/(∂(x^2)))= r u(τ,x)
or
((∂u(τ,x))/(∂τ)) =(r-q)x ((∂u(τ,x))/(∂x)) + (1/2)(σ^2)(x^2)(((∂^2)u(t,x))/(∂(x^2)))- r u(τ,x)
s.t.
u(0,x) = (SuperPlus[(x - K)]).

So for our general PDE solver we identify
a(τ,x) = (r-q) x,    b(τ,x) = (1/2)(σ^2)(x^2),    c(τ,x) = -r    and    h(x) = max(x-K,0).

Part (a)

For Part (a) we use the explicit method
(u^j) = (I + (A^(j-1))Δτ)(u^(j-1)),    j∈{1,2,…,m}.

Note that the matrix (A^j) does not depend on j.

We will see below that this scheme is numerically unstable.

Code

r=0.01;    q=0.02;    σ=0.5;    KK=20; 
τMax=1;    m=250;
dτ=τMax/m;
xMin=0;    xMax=60;    n=240;
dx=(xMax-xMin)/n;

b[x_]:=(σ^2)(x^2)/2;
a[x_]:=(r-q) x;
c[x_]:=-r;
h[x_]:=Max[x-KK,0];

AA=SparseArray[{Band[{2,1},{n,n-1}]->Table[((b[xMin+k dx])/(dx^2))-((a[xMin+k dx])/(2dx)),{k,1,n-1}],
Band[{2,2},{n,n}]->Table[c[xMin+k dx]-((2b[xMin+k dx])/(dx^2)),{k,1,n-1}],
Band[{2,3},{n,n+1}]->Table[((b[xMin+k dx])/(dx^2))+((a[xMin+k dx])/(2dx)),{k,1,n-1}]},{n+1,n+1}];
AA[[1,1;;3]]={1,-2,1};
AA[[n+1,n-1;;n+1]]={1,-2,1};

x=Table[xMin+k dx,{k,0,n,1}];
u=Table[0,{n+1},{m+1}];
u[[All,1]]=Map[h,x];

For[j=1,j≤m,j++,
u[[All,j+1]]=u[[All,j]]+AA.u[[All,j]]dτ;
];

ListPlot3D[u,AxesLabel->{"τ","(Subscript[S, τ])","V(τ,(Subscript[S, τ]))"},DataRange->{{0,τMax},{xMin,xMax}},PlotLabel->"European vanilla call pricing surface",PlotRange->All]

Part (b)

For Part (b) we use the Crank-Nicholson method
(I - (1/2) (A^j) Δτ)(u^j) = (I + (1/2) (A^(j-1))Δτ)(u^(j-1)),    j∈{1,2,…,m}.

Note that the matrix (A^j) does not depend on j.

Code

r=0.01;    q=0.02;    σ=0.5;    KK=20;
τMax=1;    m=250;
dτ=τMax/m;
xMin=0;    xMax=60;    n=240;
dx=(xMax-xMin)/n;

b[x_]:=(σ^2)(x^2)/2;
a[x_]:=(r-q) x;
c[x_]:=-r;
h[x_]:=Max[x-KK,0];

AA=SparseArray[{Band[{2,1},{n,n-1}]->Table[((b[xMin+k dx])/(dx^2))-((a[xMin+k dx])/(2dx)),{k,1,n-1}],
Band[{2,2},{n,n}]->Table[c[xMin+k dx]-((2b[xMin+k dx])/(dx^2)),{k,1,n-1}],
Band[{2,3},{n,n+1}]->Table[((b[xMin+k dx])/(dx^2))+((a[xMin+k dx])/(2dx)),{k,1,n-1}]},{n+1,n+1}];
AA[[1,1;;3]]={1,-2,1};
AA[[n+1,n-1;;n+1]]={1,-2,1};

x=Table[xMin+k dx,{k,0,n,1}];
u=Table[0,{n+1},{m+1}];
u[[All,1]]=Map[h,x];

XX=IdentityMatrix[n+1]-AA dτ/2;
For[j=1,j≤m,j++,
u[[All,j+1]]=LinearSolve[XX,u[[All,j]]+AA .u[[All,j]]dτ/2];
];

ListPlot3D[u,AxesLabel->{"τ","(Subscript[S, τ])","V(τ,(Subscript[S, τ]))"},DataRange->{{0,τMax},{xMin,xMax}},PlotLabel->"European vanilla call pricing surface",PlotRange->All]

k[x_]:=((x-xMin)/dx);    (* return matrix index for (Subscript[S, t])=x  *)
u[[k[26]+1,m+1]]

Black-Scholes equation

price[t_,TT_,St_,KK_,r_,q_,σ_,type_]:=
Module[{d1,d2,Nd1,Nd2},
d1=(1/(σ(Sqrt[TT-t])))(Log[(St/KK)]+(r-q+(σ^2)/2)(TT-t));
d2=d1-σ(Sqrt[TT-t]);
Nd1=CDF[NormalDistribution[],type*d1];
Nd2=CDF[NormalDistribution[],type*d2];
type*((E^(-q(TT-t)))St*Nd1-(E^(-r(TT-t)))KK*Nd2)
];
price[0,1,26,20,0.01,0.02,0.5,1]//N

Part (c)

For Part (c) we use the Crank-Nicholson method
(I - (1/2) (A^j) Δτ)(u^j) = (I + (1/2) (A^(j-1))Δτ)(u^(j-1)),    j∈{1,2,…,m}.

Note that the matrix (A^j) does not depend on j.

Code

r=0.01;    q=0.02;    σ=0.5;    KK=20;
τMax=1;    m=250;
dτ=τMax/m;
xMin=0;    xMax=60;    n=240;
dx=(xMax-xMin)/n;

b[x_]:=(σ^2)(x^2)/2;
a[x_]:=(r-q) x;
c[x_]:=-r;
h[x_]:=Max[x-KK,0];

AA=SparseArray[{Band[{2,1},{n,n-1}]->Table[((b[xMin+k dx])/(dx^2))-((a[xMin+k dx])/(2dx)),{k,1,n-1}],
Band[{2,2},{n,n}]->Table[c[xMin+k dx]-((2b[xMin+k dx])/(dx^2)),{k,1,n-1}],
Band[{2,3},{n,n+1}]->Table[((b[xMin+k dx])/(dx^2))+((a[xMin+k dx])/(2dx)),{k,1,n-1}]},{n+1,n+1}];
AA[[1,1;;3]]={1,-2,1};
AA[[n+1,n-1;;n+1]]={1,-2,1};
XX=IdentityMatrix[n+1]-0.5*AA dτ;

x=Table[xMin+k dx,{k,0,n,1}];
u=Table[0,{n+1},{m+1}];
payoff=Map[h,x];
u[[All,1]]=payoff;

For[j=1,j≤m,j++,
u[[All,j+1]]=LinearSolve[XX,u[[All,j]]+AA .u[[All,j]]dτ/2];
u[[All,j+1]]=Map[Max,Transpose[{u[[All,j+1]],payoff}]];    (* check early exercise condition *)
];

ListPlot3D[u,AxesLabel->{"τ","(Subscript[S, τ])","V(τ,(Subscript[S, τ]))"},DataRange->{{0,τMax},{xMin,xMax}},PlotLabel->"American vanilla call pricing surface",PlotRange->All]

k[x_]:=((x-xMin)/dx);    (* return matrix index for (Subscript[S, t])=x  *)
u[[k[26]+1,m+1]]

Part (d)

For Part (d) we use the Crank-Nicholson method
(I - (1/2) (A^j) Δτ)(u^j) = (I + (1/2) (A^(j-1))Δτ)(u^(j-1)),    j∈{1,2,…,m}.

Note that the matrix (A^j) does not depend on j.

Code

r=0.01;    q=0.02;    σ=0.5;    KK=20;
τMax=1;    m=250;
dτ=τMax/m;
xMin=0;    xMax=30;    n=120;
dx=(xMax-xMin)/n;

b[x_]:=(σ^2)(x^2)/2;
a[x_]:=(r-q) x;
c[x_]:=-r;
h[x_]:=Max[x-KK,0];

AA=SparseArray[{Band[{2,1},{n,n-1}]->Table[((b[xMin+k dx])/(dx^2))-((a[xMin+k dx])/(2dx)),{k,1,n-1}],
Band[{2,2},{n,n}]->Table[c[xMin+k dx]-((2b[xMin+k dx])/(dx^2)),{k,1,n-1}],
Band[{2,3},{n,n+1}]->Table[((b[xMin+k dx])/(dx^2))+((a[xMin+k dx])/(2dx)),{k,1,n-1}]},{n+1,n+1}];
AA[[1,1;;3]]={1,-2,1};
AA[[n+1,n-1;;n+1]]={1,-2,1};
XX=IdentityMatrix[n+1]-0.5*AA dτ;

x=Table[xMin+k dx,{k,0,n,1}];
u=Table[0,{n+1},{m+1}];
u[[All,1]]=Map[h,x];
u[[n+1,1]]=0;    (* set to zero on boundary *)

For[j=1,j≤m,j++,
u[[All,j+1]]=LinearSolve[XX,u[[All,j]]+AA .u[[All,j]]dτ/2];
u[[n+1,j+1]]=0;    (* set to zero on boundary *)
];

ListPlot3D[u,AxesLabel->{"τ","(Subscript[S, τ])","V(τ,(Subscript[S, τ]))"},DataRange->{{0,τMax},{xMin,xMax}},PlotLabel->"European up-and-out barrier call pricing surface",PlotRange->All]

k[x_]:=((x-xMin)/dx);    (* return matrix index for (Subscript[S, t])=x  *)
u[[k[26]+1,m+1]]