NDSolve to solve a PDF

Question

I would like to solve numerically the following PDE: $$ \frac{\partial}{\partial t} p(x,t)=-\frac{\partial}{\partial x}p(x,t)+\frac{1}{2}\frac{\partial^2}{\partial x^2}[x^2\,p(x,t)],\;x\ge0,\;t\ge0, $$ where $p(x,t)$ denotes a probability density function. The PDF is subject to the following conditions:

1. $p(x,0+)=\delta(x-x_0)$, where $x_0\ge0$ is a given constant and $\delta(x)$ denotes the Delta function;

2. $\int_{0}^\infty p(x,t)\,dx=1$ for all $t\ge0$;

3. Two limit-type boundary conditions: $$ \lim_{x\to0+}\left\{p(x,t)-\frac{1}{2}\frac{\partial}{\partial x}[x^2\,p(x,t)]\right\}=0, $$ and $$ \lim_{x\to+\infty}\left\{p(x,t)-\frac{1}{2}\frac{\partial}{\partial x}[x^2\,p(x,t)]\right\}=0. $$

I am a novice in Mathematica, but having searched stackexchange.com for similar questions, I was able to come up with the following code:

ClearAll["Global`*"];
$Assumptions = x >= 0 && Subscript[x, 0] >= 0 && t >= 0;

tmax = 15;
xmax = 15;
Subscript[x, 0] = 2;
sol = NDSolve[{D[p[x, t], {t, 1}] == -D[p[x, t], {x, 1}] + 
      1/2*D[x^2*p[x, t], {x, 2}],
        p[x, 0] == D[HeavisideTheta[x - Subscript[x, 0]], {x, 1}],
        ((p[x, t] - 1/2*D[x^2*p[x, t], {x, 1}]) /. x -> xmax) == 0,
    ((p[x, t] - 1/2*D[x^2*p[x, t], {x, 1}]) /. x -> 0) == 0},
   p, {t, 0, tmax}, {x, 0, xmax}];

pxt = First[p /. sol];

Plot3D[pxt[x, t], {x, 0, xmax}, {t, 0, tmax}, PlotRange -> All]

It basically yields the trivial solution $p(x,t)\equiv 0$, which is because the normalization constraint 2 is not taken into account.

My question is whether it's possible to feed a constraint of this type of NDSolve?

UPDATE This equation can be solved analytically. Its exact solution (evaluated using Mathematica) is illustrated in the following figure:

I would like to compare the exact solution and the numerical solution. What I need help with is getting the latter using Mathematica.

Answer 1

This is a long comment that answers the question, What is the problem? I don't have an answer to the question, What is the fix?, at this point.

The problem is, I believe, condition 1. I'm afraid my background in PDEs is weak, but numerically, using the method of lines (the default in this case), condition 1 is not integrated. Basically, the initial value for p[x, 0] is set equal to the value of D[HeavisideTheta[x - x0], {x, 1}], which is the value of DiracDelta[x - 2], at each point on the spatial grid. The value is used as an initial condition to begin integration along the line. The value is either 0 or undefined. Depending on xmax, the spatial grid is unlikely to include the point x == 2., so the integrator sees initial conditions of all zeros. Given those initial conditions, one should expect to get the trivial solution.

Here is a way to test it: Find a value of xmax such that the spatial grid lands on x == 2. The default grid size in this case is 49 points or 48 intervals. So choosing xmax to be a divisor of 48 that is at least two ends up evaluating DiracDelta[x - 2] at x == 2..

tmax = 15;
xmax = 3;
x0 = 2;
dh[x_?NumericQ] := (Sow[x]; DiracDelta[x - x0]);
{sol, {grid}} = Reap@NDSolve[{
     D[p[x, t], {t, 1}] == -D[p[x, t], {x, 1}] + 1/2*D[x^2*p[x, t], {x, 2}], 
     p[x, 0] == dh[x],
     ((p[x, t] - 1/2*D[x^2*p[x, t], {x, 1}]) /. x -> xmax) == 0,
     ((p[x, t] - 1/2*D[x^2*p[x, t], {x, 1}]) /. x -> 0) == 0}, 
    p, {t, 0, tmax}, {x, 0, xmax}];

NDSolve::ndnum: Encountered non-numerical value for a derivative at t == 0.`. >>
NDSolve::ndnum: ....

grid
(*
  {0., 0.0625, 0.125, 0.1875, 0.25, 0.3125, 0.375, 0.4375, 0.5, 0.5625, \
   0.625, 0.6875, 0.75, 0.8125, 0.875, 0.9375, 1., 1.0625, 1.125, \
   1.1875, 1.25, 1.3125, 1.375, 1.4375, 1.5, 1.5625, 1.625, 1.6875, \
   1.75, 1.8125, 1.875, 1.9375, 2., 2.0625, 2.125, 2.1875, 2.25, 2.3125, \
   2.375, 2.4375, 2.5, 2.5625, 2.625, 2.6875, 2.75, 2.8125, 2.875, \
   2.9375, 3.}
*)

Or one can set the spatial grid for xmax = 15 to contain 2.:

Reap@NDSolve[{
   D[p[x, t], {t, 1}] == -D[p[x, t], {x, 1}] + 1/2*D[x^2*p[x, t], {x, 2}], 
   p[x, 0] == D[HeavisideTheta[x - x0], {x, 1}],
   ((p[x, t] - 1/2*D[x^2*p[x, t], {x, 1}]) /. x -> xmax) == 0,
   ((p[x, t] - 1/2*D[x^2*p[x, t], {x, 1}]) /. x -> 0) == 0}, 
  p, {t, 0, tmax}, {x, 0, xmax},
  Method -> {"MethodOfLines", 
    "SpatialDiscretization" -> {"TensorProductGrid", 
      "MinPoints" -> 4*15 + 1, "MaxPoints" -> 4*15 + 1, 
      "DifferenceOrder" -> Automatic}}]

NDSolve::ndinpd: The initial conditions did not evaluate to an array of numbers of depth 1 on the spatial grid. Initial conditions for partial differential equations should be specified as scalar functions of the spatial variables. >>
....

---

Update: Possible workaround

One might approximate the delta function by any of the standard formulas. For instance a unit-area box:

Block[{tmax = 15,
  xmax = 15,
  x0 = 2,
  width = 1*^-2},
 sol = NDSolve[{
    D[p[x, t], {t, 1}] == -D[p[x, t], {x, 1}] + 1/2*D[x^2*p[x, t], {x, 2}], 
    p[x, 0] == width UnitBox[(x - x0) / width],
    ((p[x, t] - 1/2*D[x^2*p[x, t], {x, 1}]) /. x -> xmax) == 0,
    ((p[x, t] - 1/2*D[x^2*p[x, t], {x, 1}]) /. x -> 0) == 0}, 
   p, {t, 0, tmax}, {x, 0, xmax},
   Method -> {"MethodOfLines", 
     "SpatialDiscretization" -> {"TensorProductGrid", 
       "MinPoints" -> Round[xmax/width], 
       "MaxPoints" -> 1 + Round[10 xmax/width],
       "DifferenceOrder" -> Automatic}}]; pxt = First[p /. sol];
 Plot3D[pxt[x, t], {x, 0, xmax}, {t, 0, tmax}, 
  PlotPoints -> 100, PlotRange -> All, AxesLabel -> Automatic]
 ]

NDSolve::mxsst: Using maximum number of grid points 15001 allowed by the MaxPoints or MinStepSize options for independent variable x. >>