How Spatial Grid Boundaries Affect Crank–Nicolson PDE Solutions

# How does the diameter of the spatial grid affects the solution of a Crank Nicholson algorithm?


# How does the diameter of the spatial grid affects the solution of a Crank Nicholson algorithm?












this is my first question so I hope I express myself clearly.

I'm trying to implement an Implicit and a Crank Nicolson algorithm for the generic PDE $\partial_\tau u(\tau,x)+a \partial_x^2 u(\tau,x) +d \partial_x u(\tau,x)+ c u(\tau,x)+d(\tau,x)=0$, where a,b,c are constants, x is the log-spot of an asset with black scholes dyanmics, that is $x=\log(S/S_0)$ and $\tau=T-t$ is the time to expiry. I solve the PDE forward from time 0 to expiry T.

Actually implementing both algorithms is quite straight forward, but my solution heavily depends on the choice of my spatial grid, namely left and right end points $x_\min$ and $x_\max$ appearing in

$x_\min < x_1 < \dots < x_N < x_\max $.

That means if I vary both values the solution variies, whereas if I fix both values and vary the discretization parameters in time $\Delta \tau$ and space $\Delta x$ the solution remains nearly unaffected.

Does someone know how to interpret this bevahiour? May it be possible that this is a direct consequence of the choice of my boundary conditions? They are probably not the best choice, but unfortunately the best I have. The solution of the PDE is the difference of a risk free OTC-derivative to that where counterparty credit risk, DVA and a funding cost adjustment is included.

Best,

Simon

Edit:

Algorithm in full detail

Ok here are some details about the algorithm, but now in terms of Let us denote the values of the function at the grid points by $u_n^m= u(\tau_m,x_n)$. Then I approximate the spatial derivatives with central differences as usual and arrive at $ \partial_\tau u(\tau_m,x_n) = -a\frac{u_{n+1}^m-2u_n^m+u_{n-1}^m }{\delta x^2} - b \frac{u_{n+1}^m-u_{n-1}^m}{2\delta x}-cu_n^m -d_n^m\\[0.8em] =\left(-\frac{a}{\delta x^2}+\frac{b}{2\delta x}\right)u_{n-1}^m+\left(\frac{2a}{\delta x^2}-c\right)u_n^m+\left(-\frac{a}{\delta x^2}-\frac{b}{2 \delta x}\right) u_{n+1}^m-d_n^m $

or in more compact form

$\partial_\tau u(\tau_m,x_n)=\alpha_1 u_{n-1}^m+\alpha_2 u_n^m +\alpha_3 u_{n+1}^m -d_n^m :=g(u^m,d^m) $

The Crank Nicholson algorithm averages the forward time difference at point $(\tau_m,x_n)$ with the backward difference at $(\tau_{m+1},x_n)$

- $\partial_\tau u(\tau_m,x_n) \approx \frac{u_n^{m+1}-u_n^m}{\delta t}$ forward difference at $(\tau_m,x_n)$

- $\partial_t u(\tau_{m+1},x_n) \approx \frac{u_n^{m+1}-u_n^m}{\delta t}$ backward difference at $(\tau_{m+1},x_n)$

But values coincide and we average the spatial derivatives for both points

$\frac{u_n^{m+1}-u_n^m}{\delta t} = \frac{1}{2} g(u^m,d^m) +\frac{1}{2}g(u^{m+1},d^{m+1})$ or $u_n^{m+1} -\frac{\delta t }{2} g(u^{m+1},d^{m+1}) = u_n^m+\frac{\delta t }{2} g(u^m,d^m) $.

We can write the right hand side in matrix notation $\begin{pmatrix} u_1^m\\u_2^m \\ \vdots \\ \vdots \\ u^m_N \end{pmatrix} + \frac{\delta t }{2} \begin{pmatrix} \alpha_2 & \alpha_3 & 0 & \dots & 0 \\ \alpha_1 & \alpha_2 & \alpha_3 & & \vdots \\ 0 & \alpha_1 & \alpha_2 & \ddots & 0 \\ \vdots & 0 & \ddots & \ddots & \alpha_3 \\ 0 & 0 & & \alpha_1 & \alpha_2 \end{pmatrix} \begin{pmatrix} u_1^m\\u_2^m \\ \vdots \\ \vdots \\ u^m_N \end{pmatrix} + \frac{\delta t}{2} \begin{pmatrix} \alpha_1 u_0^m\\0 \\ \vdots \\ 0 \\ \alpha_3 u^m_{N+1} \end{pmatrix} - \frac{\delta t}{2} \begin{pmatrix} d_1^m\\d_2^m \\ \vdots \\ \vdots \\ d^m_N \end{pmatrix}$

Doing the same for the left hand side analogously and by defining By defining \begin{gather*} \boldsymbol{u}^m=\left(u_1^m,\dots,u_N^m\right)^T, \quad \boldsymbol{b}^m=\left(\alpha_1 u_0^m,0,\dots,0,\alpha_3 u_{N+1}^m\right)^T, \quad \boldsymbol{d}^m=\left(d_1^m,\dots,d_N^m\right)^T, \\ \eta= \frac{\delta t}{2} \end{gather*} we arrive at the Crank Nicholson algorithm. \begin{align*} \left(\mathbf{I}-\eta \mathbf{A}\right) \boldsymbol{u}^{m+1}- \eta \left(\boldsymbol{b}^{m+1}-\boldsymbol{d}^{m+1} \right)= \left(\boldsymbol{I}+\eta \mathbf{A}\right)\boldsymbol{u}^m + \lambda \left(\boldsymbol{b}^m-\boldsymbol{d}^m\right) \end{align*} Which can be written more compactly by \begin{align*} \mathbf{C}=&\mathbf{I}-\eta \mathbf{A} \qquad \text{and} \qquad \mathbf{D}=\mathbf{I}+\eta \mathbf{A},\\[0.6em] \boldsymbol{e}^m &=\eta` \left(\boldsymbol{b}^m+\boldsymbol{b}^{m+1}-\boldsymbol{d}^m-\boldsymbol{d}^{m+1}\right) \end{align*} the \textbf{Crank Nicolson Scheme} in matrix notation can be written in compact form \begin{align*} \mathbf{C} \boldsymbol{u}^{m+1} = \mathbf{D} \boldsymbol{u}^{m} +\boldsymbol{e}^m \end{align*}.

Boundary Conditions

The generic PDE is the result of a transformed PDE with coefficients depending on $S$.

I know by Feynman Kac theorem that the solution $U$ of the untransformed PDE satisfies

$U(t,S_t) =-c_1\int_t^T e^{-(r+\lambda_B+\lambda_C)(u-t)} \mathbb{E}_t \left[V^+(u,S(u))\right]du - c_2\int_t^T e^{-(r+\lambda_B+\lambda_C)(u-t)} \mathbb{E}_t \left[V^-(u,S(u))\right]du$

for some constants $\lambda_B,\lambda_C,c_1,c_2$. As I said $U$ is the difference from a risk free derivative $V$ and the same derivative but including counterparty credit risk. The derivative for which I want to determine the adjustment for the moment is a simple forward contract with value $V(t,S_t) = S_t - Ke^{-r(T-t)}$. I argue as follow:

We see that $V(t,S_t) \gg 0$ if $S_t \gg 0$. We assume then that \begin{align*} V(u,S_u) \gg 0 \quad \text{for all} \quad u\geq t \end{align*} Thus we can determine the upperboundary condition for $U$ by \begin{align*} U(t,S)&= -c_1\int_t^T D_{r+\lambda_B+\lambda_C}(t,u) \mathbb{E}_t \left[V(u,S(u)\right]du \\[0.8em] &= -c_1V(t,S_t)\int_t^T e^{-(\lambda_B+\lambda_C)(u-t)} du \\[0.8em] &= c_1\frac{1}{\lambda_B+\lambda_C} \left(1-e^{-(\lambda_B+\lambda_C)(T-t)}\right) \left(S_t - Ke^{-r(T-t)}\right). \end{align*} Analogously we can determine the lower boundary condition. We assume that if $S_t \to 0$ then $S_u \approx 0$ for all $u\geq t$. Therefore \begin{align*} V(u,S_u) = -Ke^{-r(T-t)} \quad \text{for all} \quad u\geq t \quad \text{as} \quad S_t \to 0 \end{align*} holds and therefore we can deduce \begin{align*} U(t,S)&= -c_2\int_t^T D_{r+\lambda_B+\lambda_C}(t,u) \mathbb{E}_t \left[V(u,S(u)\right]du \\[0.8em] &=- c_2V(t,S_t)\int_t^T e^{-(\lambda_B+\lambda_C)(u-t)} du \\[0.8em] &= c_2\frac{1 }{\lambda_B+\lambda_C} \left(1-e^{-(\lambda_B+\lambda_C)(T-t)}\right) Ke^{-r(T-t)} \end{align*} These considerations result in the boundary conditions. \begin{align*} u(\tau_m,x_{\text{max}} ) &=\frac{c_1}{\lambda_B+\lambda_C} \left(1-e^{-(\lambda_B+\lambda_C)\tau}\right) \left(x_0e^{x_{\text{max} } } - Ke^{-r\tau_m }\right) \\[0.8em] u(\tau_m,x_{\text{min}} )&= \frac{c_2}{\lambda_B+\lambda_C} \left(1-e^{-(\lambda_B+\lambda_C)\tau_m}\right) Ke^{-r\tau_m} \end{align*}