How to solve this stationary Kolmogorov Forward Equation with a Dirac Delta

dirac deltaordinary differential equationspartial differential equationsprobability distributionsstochastic-differential-equations

I have the following stationary Kolmogorov Forward Equation

$$ 0 = -\dfrac{\partial}{\partial_x}[axg(x)] – \lambda g(x) + \lambda \delta(x-\underline{x}),$$

where $a>0,\underline{x}>0$ and $\lambda>0$. Basically we have an impulse at a point $\underline{x}$ and then the density flows to the right at a rate $ax$. At each instant mass is pulled back into $\underline{x}$ at a rate $\lambda$. If I ignore the Dirac delta I can easily solve for the density $g(x)$ by integrating the Kolomogorov Forward equation.

$$g(x) = g(\underline{x})\left(\dfrac{\underline{x}}{x}\right)^{1+\frac{\lambda}{a}} $$

I then impose that all the probability mass must sum to 1 and that pins down what $g(\underline{x})$ is.

$$g(\underline{x}) = \dfrac{\lambda}{a}\underline{x}^{-1} $$

Solving this problem numerically gives me the same result as above. I have three questions.

  1. The way I solved for the stationary density completely ignores the Dirac delta term. How is that possible?

  2. A colleague mentioned that as I have an impulse at $\underline{x}$, my density should sum to one with the mass split between $g(x)$ and a mass $m$, where $m$ captures the Dirac impulse: $1=\int_{\underline{x}}^{\infty} g(x)dx + m$ . Is this correct, and how could the solution be obtained for the stationary density?

  3. Is it possible to have a closed form expression linking an initial condition $g(\underline{x},t=0)=\delta$ to the stationary solution $g(x,t=\infty)=\dfrac{\lambda}{\underline{x}a}\left(\dfrac{\underline{x}}{x}\right)^{1+\frac{\lambda}{a}}$?

How the density looks:
How the density looks

Best Answer

The heuristic way to look at problems involving diracs is to first solve the problem below and above the ciritical point because there the dirac does not play a role and then glue the solutions together appropriatly. (The correct way to handle distributions is of course a weak formulation of the problem).

Indeed, in doing so we get $$ g(x)= \begin{cases} c_0x^{-\frac{a+\lambda}{a}}&\mbox{ for } x<\underline{x}\\ c_1x^{-\frac{a+\lambda}{a}}&\mbox{ for } x>\underline{x} \end{cases} $$ Heuristically speaking the derivative of a piecewise differentiable function which has a discontinuity where the function jumps one unit up yields a direac distribution at the point of the jump. Hence in order to achieve the correct jump height we need to ensure that $$ \lim_{x\to\underline{x}^+}axg(x)- \lim_{x\to\underline{x}^-}axg(x) = \lambda\\ \Leftrightarrow a\underline{x}^{-\frac{\lambda}{a}}(c_1-c_0)=\lambda\\ \Leftrightarrow c_1=c_0+\frac{\lambda}{a}\underline{x}^{\frac{\lambda}{a}} $$

which gives the solution $$ g(x)= \begin{cases} c_0x^{-\frac{a+\lambda}{a}}&\mbox{ for } x<\underline{x}\\ (c_0+\frac{\lambda}{a}\underline{x}^{\frac{\lambda}{a}})x^{-\frac{a+\lambda}{a}}&\mbox{ for } x>\underline{x} \end{cases} $$

Finally $c_0$ can be adapted to the boundary condition. Since you did not specify a domain and boundary condition I am going to assume that

  1. the domain is $\mathbb{R}$ (otherwise a dirac at point $0$ makes little sense)
  2. the solution should satisfy $\lim_{x\to\pm\infty}g(x)=0$.

These assumptions would imply that $c_0=0$ and the solution simplifies to $$ g(x)= \begin{cases} 0&\mbox{ for } x<\underline{x}\\ \frac{\lambda}{a}\left(\frac{\underline{x}}{x}\right)^{\frac{\lambda}{a}}x^{-1}&\mbox{ for } x>\underline{x} \end{cases} $$

which is (apart from the restriction of the support, which in essence captures the dirac) exactly your solution!

I am not 100% sure about the reasoning for your second question however it is possible to evaluate the integral and one can show that

$$ \int_\mathbb{R}g(x)=\int_\underline{x}^\infty g(x)=1 $$

Regarding the third question I am honestly unsure whether or not a solution exists for the Dirac delta distribution. What I can say is that for every inital data given by a e.g. $L^\infty$-function a solution exists and can be explicitely expressed.


The formal way is to look at the weak formulation of the equation namely $$ \int_\Omega axg(x)\phi(x)'-\lambda g(x)\phi(x)dx+\lambda\phi(\underline{x})=0 $$ for a suitable set of testfunctions $\phi$. This is basically the very definition of the deriviative in the distributional sense. This weak formulation also yields the formal derivation of the previously mentioned "jump condition". Ill show you the derivations but if you are interested in understanding the concept of weak solutions and distributions I strongly recommend a book on the topic (e.g. Evans) since this is way too much information to put it here.

Nevertheless inserting a piecewise differentiable function with discontinuity at $\underline{x}$ into the weak formulation gives. \begin{align} 0&=\int_A^B axg(x)\phi(x)'-\lambda g(x)\phi(x)dx+\lambda\phi(\underline{x})\\ &=\int_A^\underline{x} axg(x)\phi(x)'-\lambda g(x)\phi(x)dx+\int_\underline{x}^B axg(x)\phi(x)'-\lambda g(x)\phi(x)dx+\lambda\phi(\underline{x})\\ &=axg(x)\phi(x)|_{x=A}^{\underline{x}^-} +\int_A^\underline{x} (\underbrace{-(axg(x))'-\lambda g(x)}_{=0})\phi(x)dx+ axg(x)|_{x=\underline{x}^+}^B+\int_\underline{x}^B (\underbrace{-(axg(x))'-\lambda g(x)}_{=0})\phi(x)dx+\lambda\phi(\underline{x})\\ \end{align} Here we have used integration by parts on the first part of the integrals and the fact that $g$ satisfies the differential equation below and above $\underline{x}$. By definition our testfunctions vanishes outside a compact intervall and we conclude: \begin{align} 0&=a\underline{x}g(\underline{x}^-)\phi(\underline{x})-aAg(A)\underbrace{\phi(A)}_{=0}+aBg(B)\underbrace{\phi(B)}_{=0}-a\underline{x}g(\underline{x}^+)\phi(\underline{x})+\lambda\phi(\underline{x})\\ &=\left(a\underline{x}g(\underline{x}^-)-a\underline{x}g(\underline{x}^+)+\lambda\right)\phi(\underline{x}) \end{align}

Related Question