# Black Scholes PDE explicit Scheme

I am currently working on the implementation of classic schemes to solve the BS PDE and it seems that I make a mistake in my code because the result looks far from the result of the BS formula.

Here is the maths of the implicit scheme:

Substituting $$\tau = T - t$$ into the Black-Scholes equation and using the chain rule, we get:

$$-\frac{\partial V}{\partial \tau} + \frac{1}{2}\sigma^2S^2\frac{\partial^2 V}{\partial S^2} + rS\frac{\partial V}{\partial S} - rV = 0$$

To solve this equation using an explicit scheme, we first need to discretize it in both time and space. Let $$\Delta t$$ be the time step and $$\Delta S$$ be the space step. We define:

$$\tau_n = n\Delta \tau \qquad\text{and}\qquad S_j = j\Delta S$$

for integers $$n$$ and $$j$$. We also define the numerical solution at the grid points $$(S_j, \tau_n)$$ as $$V_j^n \approx V(S_j, \tau_n)$$.

We approximate the derivative with finite differences as follows:

$$\frac{\partial V}{\partial \tau}(S_j,\tau_n) \approx \frac{V_j^{n+1} - V_j^n}{\Delta \tau}$$

$$\frac{\partial V}{\partial S}(S_j,\tau_n) \approx \frac{V_{j+1}^n - V_{j-1}^n}{2\Delta S}$$

$$\frac{\partial^2 V}{\partial S^2}(S_j,\tau_n) \approx \frac{V_{j+1}^n - 2V_j^n + V_{j-1}^n}{\Delta S^2}$$

Substituting these approximations into the Black-Scholes equation, we get:

$$-\frac{V_j^{n+1} - V_j^n}{\Delta \tau} + \frac{1}{2}\sigma^2 S_j^2 \frac{V_{j+1}^n - 2V_j^n + V_{j-1}^n}{\Delta S^2} + rS_j\frac{V_{j+1}^n - V_{j-1}^n}{2\Delta S} - rV_j^n = 0$$

Rearranging this equation, we obtain an explicit scheme for $$V_j^{n+1}$$:

$$V_j^{n+1} = \frac{1}{2}\sigma^2 j^2\Delta S^2 \Delta \tau \frac{V_{j+1}^n - 2V_j^n + V_{j-1}^n}{\Delta S^2} + rj\Delta S\Delta \tau\frac{V_{j+1}^n - V_{j-1}^n}{2\Delta S} + (1-r\Delta \tau)V_j^n$$

So : $$V_j^{n+1} = V_{j+1}^n \left(\frac{\Delta \tau j}{2} \left(r+ \sigma^2 j\right)\right) + V_{j}^n \left(1 - \Delta \tau \left(\sigma^2 j^2 + r\right)\right) + V_{j-1}^n \left(\frac{\Delta \tau j}{2} \left(-r + \sigma^2 j\right)\right)$$

Here is the code :

from mpl_toolkits import mplot3d
import numpy as np
import matplotlib.pyplot as plt
from scipy.stats import norm

# Parameters
S0 = 100  # Initial stock price
K = 100   # Strike price
T = 1     # Time to maturity
r = 0.05  # Risk-free rate
sigma = 0.2  # Volatility
N = 2000   # Number of time steps
M = 100  # Number of stock price steps

# Grid setup
dt = T / N
ds = S0*2 / M
s = np.arange(S0*2, 0, -ds)
t = np.arange(T, 0,-dt)
T, S = np.meshgrid(t, s)
j = np.arange(M-1,1,-1)

#Courant-Friedrichs-Lewy (CFL) condition,

if(dt<= (ds**2/(2*(sigma*S0)**2))):
print('Courant-Friedrichs-Lewy (CFL) condition is verified')
else:
print('Courant-Friedrichs-Lewy (CFL) condition IS NOT verified')

# Initial conditions
V = np.maximum(S-K, 0)

# Boundary conditions
V[-1,:] = 0
V[0,:] = S0-K*np.exp(-r*T[0,:])

# Explicit scheme
for i in range(N-2,-1,-1):
V[1:-1,i] =  V[2:,i+1] * (0.5*dt*j*( r + sigma**2*j )) + V[1:-1,i+1]*(1- dt*( (sigma**2)*j**2 + r)) + V[0:-2,i+1]*(0.5*dt*j*( -r + (sigma**2)*j ))

#Plot the result
fig = plt.figure(figsize=(20,12))
ax = plt.axes(projection='3d')
ax.plot_surface(S, T, V, cmap='coolwarm')
ax.view_init(10, 250)
ax.set_xlabel('Stock Price')
ax.set_ylabel('Time')
ax.set_zlabel('Option Value')
ax.set_title('Black-Scholes Option Value')
plt.show()


If someone is able to spot the mistake it would be a great help! Thanks a lot :)

• when i tried running your code, it complained about t array being true or false. maybe try set T = 1.0 (float) instead of 1. Not sure if this is the cause. Mar 27 at 5:18
• Solving the PDE in its original variables will inevitably cause problems for not even that large volatilities. It is better to apply the standard transformation first. For a detailed post showing the problems see here. Mar 27 at 10:55

Edit: based on comment.

In your question, you didn't specify the boundary conditions being considered.

V[0,:] = S0-K*np.exp(-r*T[0,:])


As this condition from the Put-Call parity is to approximate large S, it should be S here instead of the initial price, S0. In this case, it should be

V[0,:] = S0*2-K*np.exp(-r*T[0,:])


as per your maximum defined S.

• Indeed, when S is large compared to the strike price, the call option price is : $$C = S - K e^{-r\tau}$$ This can be understood using the call put parity. We see that if S is very large compared to K, the put price is zero, therefore: $$C = C - P = S - K e^{-r \tau}$$ Apr 10 at 14:00
• Yes, but you are using $S_0$ instead of $S$? Apr 12 at 2:10