Setting Up a Black-Scholes Diffusion PDE Simulation

# Simulating the Black Scholes PDE - as a Diffusion equation


# Simulating the Black Scholes PDE - as a Diffusion equation












i am trying to simulate the Black Scholes equation in it's diffusion equation form.

I found this stack exchange post: Transformation from the Black-Scholes differential equation to the diffusion equation - and back

Where the derivation on how to get from Black Scholes to it's diffusion form is done.

Essentially what they do is they introduce four transformations and substitute them into the equation: $$ S = e^y,\quad t = T - \tau, \quad u=e^{r\tau} C(y,\tau), \quad x = y + (r - \frac{1}{2} \sigma^2) \tau $$

resulting in $$ \frac{\partial u}{\partial \tau} = \frac{1}{2} \sigma^2 \frac{\partial^2 u}{\partial x^2}$$

From this one can then derive the boundary conditions under the new variables, from:

$$ C(S,T) = max(S-K, 0), \quad C(0,t) = 0, \quad \lim_{S\to \infty} C(S,t) \propto S$$

one gets

$$u(x,0) = u_0(x) = max(e^{\frac{1}{2}(a+1)x} - e^{\frac{1}{2}(a-1)x})$$, where $a=\frac{2r}{\sigma^2}$.

That all is fine to me - but now i am struggling to understand how one does the actual numerical simulation of this?

I know from my Computational Physics class, that we can use the second-order product formula. I thought of doing it that way:

We choose a spatial domain $x \in [x_{min}, x_{max}]$ and discretize it into N grid points with spacing $ \Delta x$. Let

$$u_i(\tau) = u(x_i,\tau), \quad i=1,...,N$$

That way we become a vector like:

$$ \begin{equation} \vec{u}(\tau) = \begin{pmatrix} u_1(\tau)\\ \vdots\\ u_N(\tau) \end{pmatrix} \label{eq:uvector} \end{equation} $$

Now we essentially choose $\vec{u}(\tau)$ according to some initial starting conditions and simulate the system by using the second order product formula. This is done by looking at $$ \begin{equation} \frac{\partial ^2 u(x_i)}{\partial x^2}|_{x_i} \approx \frac{u_{i+1} - 2u_i + u_{i-1}}{(\Delta x)^2} \label{eq:numeric1} \end{equation} $$ with $D=\frac{1}{2} \sigma^2$ and we re-write the whole second derivative in matrix form like $$\begin{equation} \frac{\partial \vec{u}(\tau)}{\partial \tau} = -\frac{D}{(\Delta x)^2}\hat{H} \vec{u}(\tau) \label{eq:diffusion_new} \end{equation}$$

The state essentially then goes through this trajectory: $$ \vec{u}(\tau + \delta_{\tau}) = e^{\alpha \hat{H} \delta_{\tau}} \vec{u}(\tau) $$

We can then decompose the matrix $\hat{H}$ into two matrices A and B and use the product formula defined as $$ \begin{equation} e^{\alpha \hat{H} \delta_\tau} \approx \left( e^{ \frac{\alpha \delta_\tau \hat{A}}{2m}} e^{ \frac{\alpha \delta_\tau \hat{B}}{m}} e^{ \frac{\alpha \delta_\tau \hat{A}}{2m}} \right)^m \label{eq:solve2} \end{equation} $$

with m being the amount of time steps we apply onto our state.

Now i am just super confused on how to setup the grid and the initial conditions properly of that system. Maybe i have some general logical flaw somewhere - i am not sure right now. But according to what i think, maybe we just set each grid point $u_i$ to be:

$$u_i(0) = max(e^{\frac{1}{2}(a+1)x} - e^{\frac{1}{2}(a-1)x})$$

But for some reason i am now fundamentally confused about everything i am doing here as my computational physics class is already 2 years in the past and i haven't been doing much physics since then. The variable x is directly related to S, as it is defined as $$x = log(S) + (r-\frac{1}{2} \sigma^2) \tau$$, meaning that we can not have a stock price S=0 as the log would explode. What is the smallest $x_i$ we choose here then and what is the maximum? Do we distribute S around the strike price and calculate the $x_i$ according to that? Because of course choosing the x_i 's is what will create our u array...

I have uploaded my calculations and process of thought in more detail here: https://smallpdf.com/file#s=17e12553-4145-42ad-a939-42be9109d0a7

I would be incredibly greatful if anyone could give me even the slightest tip or hint at what i am doing wrong or where i am confused.

Thanks!!