Solve the probability distribution using a Kolmogorov backward equation: A common case

Question

I'd like to use MMA to tackle some advanced problems. First step is to ensure that simple, known problems are solved correctly. It seems solving PDE's is still quite an art. see Laplace Transforms vs Change of Variables

pde = D[p[x, τ], τ] == -D[(μ*x)*p[x, τ], x] + D[((σ^2*x^2)/2)*p[x, τ], {x, 2}]

soln = DSolve[pde, p[x, τ], {x,τ}]

The results of this is:

DSolve[Derivative[0, 1][p][x, τ] == (-μ)*p[x, τ] + σ^2*p[x, τ] 
   - x*μ*Derivative[1, 0][p][x, τ] + 2*x*σ^2*Derivative[1, 0][p][x, τ]
   + (1/2)*x^2*σ^2*Derivative[2, 0][p][x, τ], p[x, τ], {x, τ}]

Trying with an initial condition:

sol = DSolve[{pde, p[x, 0] == DiracDelta[x - x0]}, p[x, τ], {x, τ}]

The result is:

{{p[x, τ] ->
  Integrate[DiracDelta[-1 + x0], {K[1], -∞ ,∞}]/(Sqrt[2*Pi]*σ*Sqrt[τ])}}

Boundary and Initial Conditions

The initial condition is p[x, 0] == DiracDelta[x - x0] The function p[x,t] is a probability density function, in this case two boundary conditions are p[-∞, τ] == 0 and p[∞, τ] == 0.

Expected result

This is geometric Brownian motion so the expected result should be the lognormal density.

The problem I am interested in is as posted in this question here.

I'd greatly appreciate someone showing the canonical set up in MMA for this type of problem. Specifically how to get the expected density function returned.

Answer 1

So let's try your pde.

pde = D[p[x, τ], τ] == -D[(μ*x)*p[x, τ], x] + D[((σ^2*x^2)/2)*p[x, τ], {x, 2}]

Instead of choosing a specific initial condition, let's choose a generic f[x].

sol = DSolve[{pde, p[x, 0] == f[x]}, p[x, τ], {x, τ}] //Flatten
(* {p[x,τ] -> Integrate[f[1],{K[1],-Infinity,Infinity}]/(Sqrt[2 π] σ Sqrt[τ])} *)

This answer isn't all that useful. f[1] is a constant that can be taken outside the integral, and we are left with just an integral of infinities. The denominator indicates singularity at τ = 0, so a DiracDelta as an initial condition may not be far off. The denominator also gives us a clue to try a solution:

p[x_, τ_] = f[x]/Sqrt[τ]
DSolve[pde, f[x], x] // Flatten // Simplify

and

p[x_, τ_] = (f[x] /. %)/Sqrt[τ]

The output is a rather lengthy and is not that easy to visualize in this format, but MMA found the solution, and you can see it in MMA

Update

We can solve this pde with general separation of variables.

Clear["Global`*"]

pde = D[p[x, τ], τ] == -D[(μ*x)*p[x, τ], x] + D[((σ^2*x^2)/2)*p[x, τ], {x, 2}]

Separate p in the form:

p[x_, τ_] = X[x] T[τ]

pdenew = pde/p[x, τ] // Apart // Collect[#, x D[X[x], x]] &

(*D[T[τ],τ]/T[τ]\[Equal]-μ+σ^2+(σ^2 x^2 D[X[x],x,x])/(2 X[x])+x D[X[x],x] ((2 σ^2)/X[x]-μ/X[x])*)

The lhs is dependent on τ only and the rhs is dependent on x only, so each side must be equal to a constant. That constant can be either positive, negative, or zero, and the form of the solution depends on the sign of the constant. Since I don't know the boundary conditions, I will consider all three cases here.

A zero constant

t0eq = pdenew[[1]] == 0;

DSolve[t0eq, T[τ], τ] // Flatten;

T0 = T[τ] /. (DSolve[t0eq, T[τ], τ][[1]] /. C[1] -> 1);

x0eq = pdenew[[2]] == 0;

DSolve[x0eq, X[x], x] // Flatten // Simplify;

X0 = X[x] /. % /. {C[1] -> c1, C[2] -> c2}

I won't show the outputs here because they are lengthy and a mess in this format. Next positive constant.

tpeq = pdenew[[1]] == α^2

DSolve[tpeq, T[τ], τ] // Flatten

Tp = T[τ] /. % /. C[1] -> 1

xpeq = pdenew[[2]] == α^2

DSolve[xpeq, X[x], x] // Flatten // Simplify

Xp = X[x] /. % /. {C[1] -> c3, C[2] -> c4}

Now the negative constant

tmeq = pdenew[[1]] == -β^2

DSolve[tmeq, T[τ], τ] // Flatten

Tm = T[τ] /. % /. C[1] -> 1

xmeq = pdenew[[2]] == -β^2

DSolve[xmeq, X[x], x] // Flatten // Simplify

Xm = X[x] /. % /. {C[1] -> c5, C[2] -> c6}

and the total solution:

p[x_, τ_] = T0 X0 + Tp Xp + Tm Xm

In general the boundary conditions will allow many pieces of the solution to be quickly thrown out.