I am currently working on a Masters project in which I will be looking at numerical methods for a specific PDE. The PDE that I will be considering is a Fokker-Planck equation on a unit hyper-cube.

More specifically I am studying the following PDE (in non-divergence form) $$ \frac{\partial u}{\partial t} = \frac{x_{d}(1-x_{d})}{2}\frac{\partial^2u}{\partial x_{d}^2}+\left(1-2x_{d}-M_{d}(x)\right)\frac{\partial u}{\partial x_{d}}-(1-\omega_{d})u, $$ for all $x \in [0,1]^D$ and $t\in[0,\infty)$. The boundary conditions will be zero flux except at the $(0,...,0)$ and $(1,...,1)$ corners.

where summation over $d$ is assumed on the right side of the equation. We also define $$ M_d(x)=\sum_{k=1}^{D} M_{dk}(x_k-x_d), \ \ \ \text{and}\ \ \ \omega_d = M_d(x)=\sum_{\substack{k=1\\k\neq d}}^{D} M_{dk}. $$ The flux vector of this equation is given by $$ F_d(u,x) = \frac{1}{2}\frac{\partial}{\partial x_d}\left[x_\alpha(1-x_d)u\right]+\sum_{k=1}^{D} M_{dk}(x_k-x_d)u. $$

I will be using the finite volume method (FVM) to approximate a discrete solution to this PDE. Had this been a constant coefficient problem, or even a positive coefficient problem, I would be fine. However, I am not sure how to deal with the fact that the leading order coefficients go to zero on the boundaries of our domain.

Based on my intuition, I have a hunch that we should be able to show that the solution is a smooth function in the interior of the domain with delta distributions on the boundaries. Sadly, I am not sure how to show this or if this is correct.

Any references or suggestions that you think may be helpful would be appreciated.

Thanks in advance for the help!