# Errors on Finite Differences + Implicit Scheme + Black & Scholes

I'm solving the classical Black & Scholes (BS) PDE for a European option using finite difference and the implicit scheme. In other words, I'm trying to solve

$$\displaystyle\frac{\partial V}{\partial t} + \frac 12\sigma^2\frac{\partial^2V}{\partial x^2} + \mu\frac{\partial V}{\partial x} - rV = 0$$

with boundary condition $$V(T,x) = \Phi(e^x) = \max\{e^x - K, 0\}$$. The drift parameter is as usual $$\mu = r - y - \frac 12 \sigma^2$$.

The problem is when my $$x$$ grid is too narrow, my option prices diverge from the BS. Specifically, I'm solving the ATM case with S0 = K = 2775 and if I use [np.log(2675), np.log(2875)] as the boundary points, the solution has a huge error (~400). However, when I use a wider grid [np.log(1775), np.log(3775)], it works like a charm.

I thought the implicit scheme would be numerically stable enough to not exhibit this kind of behavior, so should I assume there is something wrong with my Python code?

My other parameters are r, y, T, sigma = (0.0278, 0.0189, 1, 0.15).

def buildDs(nx):
D1 = np.diag(np.ones(nx), 1) - np.diag(np.ones(nx), -1)
D2 = np.diag(np.ones(nx), 1) + np.diag(np.ones(nx), -1) - 2 * np.eye(nx + 1)

D1[0, 2] = -1
D1[0, 1] = 4
D1[0, 0] = -3

D1[-1, -1] = 3
D1[-1, -2] = -4
D1[-1, -3] = 1

D2[0, 2] = 1
D2[0, 1] = -2
D2[0, 0] = 1

D2[-1, -1] = 1
D2[-1, -2] = -2
D2[-1, -3] = 1

return csc_matrix(D1), csc_matrix(D2)

def pxFD(Phi, r, mu, sigma,
xs, ts,
nx = 2000, nt = 10000):

dx = (xs[1] - xs[0]) / nx
dt = (ts[1] - ts[0]) / nt

xs = np.linspace(xs[0], xs[1], num = nx + 1, endpoint = True)
ts = np.linspace(ts[0], ts[1], num = nt + 1, endpoint = True)

V = np.zeros((nt + 1, nx + 1))
I = identity(nx + 1)

D1, D2 = buildDs(nx)
L = 1 / 2 * (sigma / dx) ** 2 * D2 + mu / (2 * dx) * D1 - r * I
P = I - dt * L

V[-1] = Phi(xs)
for j in reversed(range(nt)):
V[j] = spsolve(P, V[j + 1])
return V, xs, ts

def pxFDGBM(Phi, r, y, sigma,
Ss, ts,
nx = 2000, nt = 10000):

mu = r - y - 1 / 2 * sigma ** 2
xs = np.log(Ss)
Philog = lambda x: Phi(np.exp(x))

Vs, xs, ts = pxFD(Phi = Philog,
r = r,
mu = mu,
sigma = sigma,
xs = xs,
ts = ts,
nx = nx,
nt = nt)
return Vs, np.exp(xs), ts


The PDE is defined for $$x \in ]-\infty, +\infty[$$ but the finite difference scheme requires a truncated domain $$[x_{\min}, x_{\max}]$$, and the choice of $$x_{\min}$$ and $$x_{\max}$$ will affect the quality of the result, regardless of the scheme being explicit, implicit, or mixed.
A good rule of thumb is to choose the truncation $$[x_{\min}, x_{\max}]$$ such that the probability of $$X_t < x_{\min}$$ or $$X_t > x_{\max}$$ is very low under the original SDE, because these are states that won't have much impact on the option value.
In the BS case a truncation at $$\pm 5$$ standard deviation is usually a good choice. In you example that would be $$x_{\min} = \ln(S_0) - 5 \sigma \sqrt{T} \approx \ln(1311)$$ and $$x_{\max} = \ln(S_0) + 5 \sigma \sqrt{T} \approx \ln(5875)$$.
Another thing that can improve accuracy is the choice of the boundary condition at $$x_{\min}$$ and $$x_{\max}$$.