# Numerically solved PDE of Ornstein–Uhlenbeck process on 2-Simplex violates conservation of probability

I'm working to create a solution of an Ornstein-Uhlenbeck process with a force that takes mass towards the centre of a Simplex. I'm assuming absorbing boundaries.

The Mathematica code below quickly provides a solution. However, the probability mass within the domain grows significantly early, but it should only ever diminish, due to mass being absorbed.

I don't think the error lies in my formulation of the forward Kolmogorov (Fokker-Plank) equation.

If it doesn't lie there, I suppose it could lie in the numerical approximation, with errors growing? I greatly appreciate any insight into this problem.

Update: This is my first foray into numerically solving PDEs, I am learning as I'm troubleshooting. I have tried drastically refining the mesh, accuracy, and step size. This hasn't lessened the violation of conservation, perhaps indicating a problem in formulation. But I've been fastidious with that as well, and can't find a problem. Indeed, my formulation of the differential equation appears to matched that provided automatically by Mathematica (see below).

ClearAll["Global*"]

\[Eta] = 5.; (*side length*)
xopt = {\[Eta]/2, \[Eta]/(2  Sqrt[3])}; (*centroid*)
\[Kappa] = .75; (*rate of reversion to centroid,diffusion constant=1*)
Tmax = 5.; (*length of time*)

\[CapitalOmega] =
Polygon[Rationalize[{{0, 0}, {\[Eta],
0}, {\[Eta]/2, (\[Eta]   Sqrt[3])/2}},
0]];  (*domain is equilateral triangle*)

bC = Rationalize[DirichletCondition[P[x1, x2, t] == 0, True],
0]; (*absobing boundary condition*)
iC = Rationalize[
P[x1, x2, 0] ==
Piecewise[{{1/((Sqrt[3]  \[Eta]^2)/4),
RegionMember[\[CapitalOmega], {x1, x2}]}}, 0],
0]; (*uniform initial condition*)

(*forward Kolmogorov equation*)
fwrdKol =
Rationalize[
D[P[x1, x2, t],
t] == -D[\[Kappa]  (xopt[[1]] - x1)*P[x1, x2, t], {x1, 1}] -
D[\[Kappa]  (xopt[[2]] - x2)*P[x1, x2, t], {x2, 1}] +
1/2  D[P[x1, x2, t], {x1, 2}] + 1/2  D[P[x1, x2, t], {x2, 2}], 0];

(*numerical solution*)
Psol = NDSolveValue[{fwrdKol, iC, bC},
P, {x1, x2} \[Element] \[CapitalOmega], {t, 0, Tmax}];

(*visualise solution at a t=Tmax/2*)
ContourPlot[Psol[x1, x2, Tmax/2], {x1, x2} \[Element] \[CapitalOmega]]

(*probability mass within domain*)
domP[t_] :=
NIntegrate[
Rationalize[Psol[x1, x2, t],
0], {x1, x2} \[Element] \[CapitalOmega], AccuracyGoal -> 4]

(*visualise*)
Plot[domP[t], {t, 0, 5}, PlotTheme -> "Scientific", PlotRange -> All,
FrameLabel -> {"t", "Prob. Mass Domain"}]


Automatically generated PDE:

eqn = ItoProcess[{\[DifferentialD]x1[
t] == \[Kappa] (x1opt -
x1[t])  \[DifferentialD]t + \[DifferentialD]w1[
t], \[DifferentialD]x2[
t] == \[Kappa] (x2opt -
x2[t])   \[DifferentialD]t + \[DifferentialD]w2[t]}, {x1[t],
x2[t]}, {{x1, x2}, {x10, x20}},
t, {w1 \[Distributed] WienerProcess[],
w2 \[Distributed] WienerProcess[]}]



$$\ \partial_t P =-\left(\frac{3}{2}\ - \ y \right) \ \partial_y \ P \ -\left(\frac{5}{2}-x \right) \partial_x \ P +\frac{1}{2}\Delta \ P \ + \frac{3}{2} \ P$$ by the fact that the directional field multiplies P inside the first order differential $$\vec \nabla (\vec x - \vec x_0) P$$ For the given constant start distribution, for some time, all terms except the linear $$P$$ term are zero, so there is an exponential increase in the beginning. From the wikipedia article, the Kolmogorov PDE has the vector field outside $$\nabla$$.

• Thanks @Roland F for you input. If I understand you point correctly, I think you may be referencing the Backward Kolmogorov PDE, which tells you from an initial point, how likely are you to end up at another. I'm instead using the Forward Kolmogorov PDE. $$\frac{\partial P}{\partial t} + \nabla \cdot \mathbf{J} = 0$$ where $$J_i = \kappa (\hat{x}_i - x_i) P - \frac{1}{2} \frac{\partial}{\partial x_i} P$$ So that the linear term should show up, if I'm not mistaken. Commented Jul 12 at 16:13

The initial condition was inconsistent and needed to exclude the boundary:

iC = Rationalize[
P[x1, x2,0] ==
Piecewise[{{1/
Area[\[CapitalOmega]], (RegionMember[
RegionBoundary[\[CapitalOmega]], {x1, x2}] == False \[And]
RegionMember[\[CapitalOmega], {x1, x2}])}}, 0], 0];
`