# NDSolve to solve a PDF

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. • Let's ignore condition 2 for a moment. Since you specify one boundary condition for$t$and two for$x$, isn't the solution fully determined...? Is there any other solution? – Marius Ladegård Meyer Apr 15 '15 at 16:34 • In this case the normalization constraint is important, for the problem is singular at the end points for$x$. – Alex Apr 15 '15 at 17:05 • So what solution are you expecting? – Marius Ladegård Meyer Apr 15 '15 at 20:40 • It's actually a well-known equation and its closed-form analytic solution is also known. The solution involves a$[0,\infty)\$-integral of the product of two Whittaker W-functions (integration is with respect to the first index of the W-functions). – Alex Apr 15 '15 at 23:45
• Calculate that analytic solution for these parameters with Mathematica and compare. As I said, the system is fully determined without constraint 2, so including it makes the system overdetermined. Btw, you should immediately replace Subscript[x,0] with x0, because Subscript is not atomic in Mathematica, so NDSolve will "see" the x inside it and treat it like the independent variable of p. – Marius Ladegård Meyer Apr 16 '15 at 4:59

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. >> • One can also use the initial condition p[x, 0] == PDF[NormalDistribution[0, width], x - x0]. At first, it would finish running. But restarting the kernel fixed that. – Michael E2 Apr 16 '15 at 16:25
• Many thanks. Let me play with this now to see how this solution compares to the exact solution. – Alex Apr 16 '15 at 19:07
• @Alex You're welcome. I just realized I left a "not" out of my comment -- at first, the NDSolve would not finish with a normal approximation to delta. That seemed to be one-time glitch. – Michael E2 Apr 16 '15 at 20:06
• I actually didn't pay attention to that. What do you mean by "it wouldn't finish running"? – Alex Apr 16 '15 at 22:05
• "Numerically, using the method of lines (the default in this case), condition 1 is not integrated." Yeah, actually this has been (vaguely) mentioned in the possible issue of DiracDelta: Numerical routines will typically miss the contributions from measures at single points. – xzczd Apr 18 '15 at 5:34