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

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.

• Seems that DSolve can solve your equation in v11.2, what's your problem? (BTW, the equation in your code is different from the one in the linked post. ) Commented Nov 30, 2017 at 10:58
• Thank you for looking at this. I've added the output I see. I'm using MMA 11.1. Do you mean you make no changes to the inputs I have and get the expected result from 11.2? Commented Nov 30, 2017 at 19:27
• The output in v11.1 is the same as in v11.2, but as mentioned in my last comment and in Bill Watts' answer below, the equation in the link is different from the one in your post, so we don't know what's the expected result at the moment. Commented Dec 1, 2017 at 5:55
• So, you mean the answer given by DSolve is incorrect? Do you know what the correct answer is? Commented Dec 1, 2017 at 8:20
• I have updated the question to now include what the answer should be. Commented Dec 16, 2017 at 12:40

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.

• @Hedgehog, I have changed my answer to address your pde of interest Commented Dec 4, 2017 at 22:10
• It seems that people do not use MMA for this domain of problems. I'm starting to see why - its not MMA but the PDE approach requires you to already know so much about the solution (e.g. the appropriate transformations) that you rpobably already know the solution. That is I've only seen people use the PDE approach when the solution is known and someone is demonstrating that solving the PDE can get you the same answer. I hope that helps anyone considering embarking along this route. Commented Dec 31, 2017 at 3:51
• I think MMA can solve many problems where we don't know the solutions, but yes, you do need some insight about what type solution your looking for. One pde can have many form of solution. In this case, the problem I had was with the DiracDelta IC, where p is 0 everywhere initially except for x0. Also, you sort of have a diffusivity equation, but not quite. Commented Dec 31, 2017 at 4:30