4
$\begingroup$

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: enter image description here

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

$\endgroup$
  • $\begingroup$ 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? $\endgroup$ – Marius Ladegård Meyer Apr 15 '15 at 16:34
  • $\begingroup$ In this case the normalization constraint is important, for the problem is singular at the end points for $x$. $\endgroup$ – Alex Apr 15 '15 at 17:05
  • $\begingroup$ So what solution are you expecting? $\endgroup$ – Marius Ladegård Meyer Apr 15 '15 at 20:40
  • $\begingroup$ 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). $\endgroup$ – Alex Apr 15 '15 at 23:45
  • $\begingroup$ 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. $\endgroup$ – Marius Ladegård Meyer Apr 16 '15 at 4:59
4
$\begingroup$

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. >>

Mathematica graphics

$\endgroup$
  • $\begingroup$ 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. $\endgroup$ – Michael E2 Apr 16 '15 at 16:25
  • $\begingroup$ Many thanks. Let me play with this now to see how this solution compares to the exact solution. $\endgroup$ – Alex Apr 16 '15 at 19:07
  • $\begingroup$ @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. $\endgroup$ – Michael E2 Apr 16 '15 at 20:06
  • $\begingroup$ I actually didn't pay attention to that. What do you mean by "it wouldn't finish running"? $\endgroup$ – Alex Apr 16 '15 at 22:05
  • 1
    $\begingroup$ "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. $\endgroup$ – xzczd Apr 18 '15 at 5:34

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.