Numerical scheme and boundary condition for $2$-D Fokker Planck equation $\newcommand{\P}{\mathbb{P}}$ I have a $2$-D stationary Fokker-Planck equation
$$\frac{\partial^2 \P(A,B)}{\partial A^2}+\frac{\partial^2 \P(A,B)}{\partial B^2}=\frac{\partial f_1(A,B) \P(A,B)}{\partial A}+\frac{\partial f_2(A,B) \P(A,B)}{\partial B},$$
where $f_1$ and $f_2$ are polynomials.
I would like to compute the solution with the normalization condition
$\int \P(A,B)\mathrm{d}A\mathrm{d}B=1$. I have two question:

*

*Is this possible without any boundary conditions?


*If I chose to use finite differences method, a very naive choice would be to do something like this:
$$\frac{\P_{(a+1,b)}-2\P_{(a,b)}+\P_{(a-1,b)}}{\Delta A^2}+\frac{\P_{(a,b+1)}-2\P_{(a,b)}+\P_{(a,b-1)}}{\Delta B^2}\\=\frac{\P_{(a+1,b)}f_{1(a+1,b)}-\P_{(a-1,b)}f_{1(a-1,b)}}{2\Delta A}+\frac{\P_{(a,b+1)}f_{2(a,b+1)}-\P_{(a,b-1)}f_{2(a,b-1)}}{2\Delta B},$$
where $(a,b)$ are indexes of 2D uniform mesh, with an extra equation given by the discretization of the normalization condition.
Is this a good way to solve the problem ?
 A: Your equation may be rewritten in somewhat more standard terms by viewing it
as a problem on $\mathbb{R}^2$. For this, just compute the derivatives
on the left-hand side:
$$
\Delta p = \partial_A (f_1 p) + \partial_B (f_2 p) \\
= (\partial_A f_1) p + f_1 \partial_A p +
(\partial_B f_2) p + f_2 \partial_B p \\
= \left[f_1 \partial_A p + f_2 \partial_B p\right] + 
\left[(\partial_A f_1 + \partial_B f_2) p\right].
$$
Introducing $F = (f_1, f_2)^T$, you may rewrite the equation as
$$
-\Delta p + F \cdot (\nabla p) + (\nabla \cdot F) p = 0
$$
This is a stationary convection-diffusion-reaction equation. In theory,
you don't have boundary conditions for $p$, however the normalisation condition implies that $p$ will decay towards $0$ for
$A, B \rightarrow \pm \infty$. So an approach that could work is to
compute the solution on a "sufficiently large" box $$B = (-L, L) \times (-L, L) \subset \mathbb{R}^2, \quad L > 0,$$ where you impose homogeneous Dirichlet boundary conditions, i.e. $p|_{\partial B} = 0$.
Proving existence and uniqueness in the above setting depends on the
polynomials $f_1$ and $f_2$, but I assume in general, the solution will
only be unique modulo scaling by a constant. Although the normalisation constraint could be incorporated in a Finite Element framework (using Mixed Finite Elements more specifically), I think
that it might be worth a try to simply apply the Finite Difference Method
to the equation with Dirichlet boundary conditions and afterwards compute
the integral
$$C = \int\limits_B \! \mathrm{d}x~p$$
numerically and then scale your obtained solution by $C^{-1}$.
For specific choices of polynomials $f_1$ and $f_2$ there might even be
simpler solutions to this problem, e.g. in the case that $\nabla \cdot F = 0$.
