Solving a Fokker-Planck equation NDSolveValue

Question

I want to solve the following Fokker-Planck equation:

$$\partial_t p = \frac{1}{2}\partial_{xx}\left[\varepsilon^2p\right] + \frac{1}{2}\partial_{yy}\left[\varepsilon^2p\right] - \partial_x\left[(9x-x^3)p\right] - \partial_y\left[-yp \right]$$

where $\varepsilon \in (0,1)$. The initial condition is a Dirac delta function at the origin.

The code I have so far is

delt[x_, y_] := 
  PDF[MultinormalDistribution[{0, 0}, {{.0015, 0}, {0, .0015}}], {x, 
    y}];
MvFP = \!\(
\*SubscriptBox[\(\[PartialD]\), \(t\)]\(p[x, y, 
     t]\)\) == .5*\[Epsilon]^2*\!\(
\*SubscriptBox[\(\[PartialD]\), \(x, x\)]\(p[x, y, 
       t]\)\) + .5*\[Epsilon]^2*\!\(
\*SubscriptBox[\(\[PartialD]\), \(y, y\)]\(p[x, y, t]\)\) - \!\(
\*SubscriptBox[\(\[PartialD]\), \(x\)]\((\((9*x - 
\*SuperscriptBox[\(x\), \(3\)])\)*p[x, y, t])\)\) + \!\(
\*SubscriptBox[\(\[PartialD]\), \(y\)]\((y*p[x, y, t])\)\);
\[Epsilon] = .2;
sol = NDSolveValue[
   {MvFP,
    p[x, y, 0] == delt[x, y],
    p[-5, y, t] == 0,
    p[5, y, t] == 0,
    p[x, -5, t] == 0,
    p[x, 5, t] == 0
    },
   p, {t, 0, 1}, {y, -5, 5}, {x, -5, 5}, 
   Method -> {"MethodOfLines", "TemporalVariable" -> t, 
     "SpatialDiscretization" -> "FiniteElement"}
   ];
Manipulate[Plot3D[
  sol[x, y, t], {x, -5, 5}, {y, -5, 5},
  PlotRange -> All,
  PlotPoints -> 50,
  ColorFunction -> "Rainbow"
  ],
 {t, 0, 1}
 ]

The resulting plot makes little sense. If I put the boundary at $x\in (-1,1)$ and $y \in (-1,1)$, the plot of the solution makes sense, but I would like to be able to visualize the solution over a larger region.

I suspect the issues have something to do with the options for the "Method" of NDSolveValue, but after seeing some other posts, I don't have any good ideas how to pick appropriate options.

Answer 1

Rule of thumb: if NDSolve doesn't work well on solving PDE, and you're sure you haven't made any simple mistake, then usually the problem lies in spatial discretization. When TensorProductGrid is chosen for spatial discretization, usually we adjust MinPoints and MaxPoints sub-option, here is an example. When FiniteElement is chosen for spatial discretization, usually we adjust MaxCellMeasure sub-option.

molfem[measure_: Automatic] := {"MethodOfLines", 
   "SpatialDiscretization" -> {"FiniteElement", 
     "MeshOptions" -> MaxCellMeasure -> measure}};
sol = NDSolveValue[
       {MvFP,
         p[x, y, 0] == delt[x, y], DirichletCondition[p[x, y, t] == 0, True]
         },
       p, {t, 0, 1}, {y, -5, 5}, {x, -5, 5}, Method -> molfem[0.01]
       ]; // AbsoluteTiming
(* {4.31252, Null} *)

molfem is a function I keep in my SystemOpen@"init.m" file because adjustion of MaxCellMeasure option is so frequently needed when playing with NDSolve. DirichletCondition[p[x, y, t] == 0, True] is equivalent to your original b.c.s but conciser.