Consider the following PDE \begin{align} V(x,y) = F(x,y) + \frac{\partial V}{\partial x} f(x) + \frac{\partial V}{\partial y} g(y) \end{align}
where $x$ and $y$ are states, $V$ is the value function, $F$ the payoff function and $f$ and $g$ are the law of motions.
By $v^n$ we denote some approximation of $V$, where $n = 0,1,\ldots$ is the index of iterations. U use value function iteration and update via an implicit scheme \begin{align} \frac{v^{n+1} - v^n}{\Delta} + v^{n+1} = F + \frac{\partial v^{n+1}}{\partial x} f + \frac{\partial v^{n+1}}{\partial y} g \end{align} where $\Delta$ is some scaling parameter (CFL condition). I try to approximate $V$ over some finite and equidistant grid $(x_i,y_j)\in \{x_1,\ldots,x_I\}\times\{y_j, \ldots,y_J\}$. I have no boundary conditions, thus I use at $i = 1$ and $j = 1$ the respective forward difference and at $i=I$ and $j = J$ the backward difference. In between I use the central difference. That specification gives me 7 different equations for calculating $v^{n+1}$ dependending on $i$ and $j$. In particular
At $i=j=1$ \begin{align} \frac{v^{n+1}_{1,1} - v^n_{1,1}}{\Delta} + v^{n+1}_{1,1} = F_{1,1} + \frac{v^{n+1}_{2,1} - v^{n+1}_{1,1}}{\Delta x} f_1 + \frac{v^{n+1}_{1,2} - v^{n+1}_{1,1}}{\Delta y} g_1 \end{align}
At $i=1$, $J > j > 1$ \begin{align} \frac{v^{n+1}_{1,j} - v^n_{1,j}}{\Delta} + v^{n+1}_{1,j} = F_{1,j} + \frac{v^{n+1}_{2,j} - v^{n+1}_{1,j}}{\Delta x} f_1 + \frac{v^{n+1}_{1,j+1} - v^{n+1}_{1,j-1}}{2\Delta y} g_j \end{align}
At $I > i > 1$, $j = 1$ \begin{align} \frac{v^{n+1}_{i,1} - v^n_{i,1}}{\Delta} + v^{n+1}_{i,1} = F_{i,1} + \frac{v^{n+1}_{i+1,1} - v^{n+1}_{i-1,1}}{2\Delta x} f_i + \frac{v^{n+1}_{i,2} - v^{n+1}_{i,1}}{\Delta y} g_1 \end{align}
At $I > i>1$, $J > j > 1$ \begin{align} \frac{v^{n+1}_{i,j} - v^n_{i,j}}{\Delta} + v^{n+1}_{i,j} = F_{i,j} + \frac{v^{n+1}_{i+1,j} - v^{n+1}_{i-1,j}}{2\Delta x} f_i + \frac{v^{n+1}_{i,j+1} - v^{n+1}_{i,j-1}}{2\Delta y} g_j \end{align}
At $i=I$, $J > j > 1$ \begin{align} \frac{v^{n+1}_{I,j} - v^n_{I,j}}{\Delta} + v^{n+1}_{I,j} = F_{I,j} + \frac{v^{n+1}_{I,j} - v^{n+1}_{I-1,j}}{\Delta x} f_I + \frac{v^{n+1}_{I,j+1} - v^{n+1}_{I,j-1}}{2\Delta y} g_j \end{align}
At $I > i > 1$, $j = J$ \begin{align} \frac{v^{n+1}_{i,J} - v^n_{i,J}}{\Delta} + v^{n+1}_{i,J} = F_{i,J} + \frac{v^{n+1}_{i+1,J} - v^{n+1}_{i-1,J}}{2\Delta x} f_i + \frac{v^{n+1}_{i,J} - v^{n+1}_{i,J-1}}{\Delta y} g_J \end{align}
At $i = I$, $j = J$ \begin{align} \frac{v^{n+1}_{I,J} - v^n_{I,J}}{\Delta} + v^{n+1}_{I,J} = F_{I,J} + \frac{v^{n+1}_{I,J} - v^{n+1}_{I-1,J}}{\Delta x} f_I + \frac{v^{n+1}_{I,J} - v^{n+1}_{I,J-1}}{\Delta y} g_J \end{align} where $\Delta x := x_i - x_{i-1}$ and $\Delta y := y_j - y_{j-1}$.
I know that for the one dimensional case we would have three euqations ($i=1$, $1 < i <I$, $i = I$) and could write this more compact in matrix notation and solve $v^{n+1}$ by matrix inversion. However, I don't know how to write down this system compactly.
- Can you point me to literature where those kind of problems are discussed?
- Do have an idea how to compute $v^{n+1}$ efficiently?
I was thinking about to extend the one dimensional method. Firstly we may collect terms with $v^{n+1}_{(\cdot)}$ and have something like
\begin{align} a v^{n+1}_{i-1,j-1} + b v^{n+1}_{i,j-1} + c v^{n+1}_{i-1,j} + d v^{n+1}_{i,j} + e v^{n+1}_{i+1,j} + h v^{n+1}_{i,j+1} + l v^{n+1}_{i+1,j+1} = F_{i,j} + \frac{v^n_{i,j}}{\Delta} \end{align} where greek symbols denote coeffcients.
Where for instance $a$ is given by \begin{align} a = \begin{cases} 0, &i=j=1\\ 0, &i=1,J>j>1\\ 0, &I>i>1,j=1\\ \frac{f_i}{2\Delta x} + \frac{g_j}{2\Delta y},\quad & I>i>1,J>j>1\\ \frac{f_i}{\Delta x} + \frac{g_j}{2\Delta y},\quad & i=I,J>j>1\\ \frac{f_i}{2\Delta x} + \frac{g_j}{\Delta y},\quad & I>i>1,j=J\\ \frac{f_i}{\Delta x} + \frac{g_j}{\Delta y},\quad & i=I,j=J\\ \end{cases} \end{align}
Now I'd like to get a linear system of the form $Av^{n+1}=B$ where $A$ is a $IJ \times IJ$ matrix of coefficients $\{a,b,c,d,e,h,l\}$ and $v^{n+1}$ and $B$ are $IJ \times 1$ vectors, such that we can solve for $v^{n+1}=A^{-1}B$.
Example: $I = J = 3$
I worked through a minimal example. Note that $v'_{(\cdot)}:=v^{n+1}_{(\cdot)}$. \begin{align} \underbrace{ \begin{bmatrix} d & h & 0 & e & l & 0 & 0 & 0 & 0\\ b & d & h & 0 & e & l & 0 & 0 & 0\\ 0 & b & d & 0 & 0 & e & 0 & 0 & 0\\ c & 0 & 0 & d & h & 0 & e & l & 0\\ a & c & 0 & b & d & h & 0 & e & l\\ 0 & a & c & 0 & b & d & 0 & 0 & e\\ 0 & 0 & 0 & c & 0 & 0 & d & h & 0\\ 0 & 0 & 0 & a & c & 0 & b & d & h\\ 0 & 0 & 0 & 0 & a & c & 0 & b & d \end{bmatrix}}_{=:A} \times \underbrace{ \begin{bmatrix} v'_{11}\\ v'_{12}\\ v'_{13}\\ v'_{21}\\ v'_{22}\\ v'_{23}\\ v'_{31}\\ v'_{32}\\ v'_{33} \end{bmatrix}}_{=:v'} = \underbrace{ \begin{bmatrix} B_{11}\\ B_{12}\\ B_{13}\\ B_{21}\\ B_{22}\\ B_{23}\\ B_{31}\\ B_{32}\\ B_{33} \end{bmatrix}}_{=:B} \end{align}
There is some pattern in filling $A$ (look at the nine $3 \times 3$ matrices). Any idea how to automatize the filling by a double loop over $i$ and $j$?
THIS IS A COMNENT, but much too long to be put into the comments section.
The calculus (such as done by clueless) might be simplified thanks to preliminary change of variables :
$$V(x,y) = F(x,y) + \frac{\partial V}{\partial x} f(x) + \frac{\partial V}{\partial y} g(y)$$
Let : $\begin{cases} X=\int \frac{dx}{f(x)} \\ Y=\int \frac{dy}{g(x)} \end{cases}\quad $ which allows the computation of the functions $x=x(X)$ and $y=y(Y)$
So, $V(x,y)=V\left(x(X),y(Y)\right)= V(X,Y)$ and $F(x,y)=F\left(x(X),y(Y)\right)= F(X,Y)$
The PDE is transformed into : $$\frac{\partial V}{\partial X} + \frac{\partial V}{\partial Y}-V(X,Y)=-F(X,Y)$$ This is a first order linear non-homogeneous PDE.
The associated homogeneous PDE : $$\frac{\partial V}{\partial X} + \frac{\partial V}{\partial Y}-V(X,Y)=0$$ can be solved thanks to the method of characteristics. The set of characteristic equations is : $$\frac{dX}{1}=\frac{dY}{1}=\frac{dV}{V}$$ from which two equations of characteristic curves can be derived :
First : $\quad X-Y=c_1$
Second : $\quad \frac{V}{X}=c_2$
The general solution of the homogeneous ODE on implicit form is : $$f\left( X-Y \:,\:\frac{V}{X} \right)=0$$ where $f$ is any differentiable function of two variables.
Or, expressed on explicit form : $$V=X \: \Phi(X-Y)$$ where $\Phi$ is any differentiable function of one variable.
In order to find the general solution of the inhomogeneous PDE, one have to compute a particular solution $V_p$ of $$\frac{\partial V_p}{\partial X} + \frac{\partial V_p}{\partial Y}-V_p(X,Y)=-F(X,Y)$$ The general solution of the inhomogeneous PDE is : $$V(X,Y)=V_p+X \: \Phi(X-Y)$$
The finite differences method of computation can be used to compute an approximate of $V_p$
As pointed out at the begining, this is not an answer to the question raised by clueless, but a comment suggesting a preliminary change of variables in order to simplify the form of the PDE. I don't know if this can be useful in the context of discretized calculus.