Solution to the Smoluchowski equation for a harmonic potential.

PHYS - 436 STATISTICAL PHYSICS IV
Jeremy Stevens, Tobias Kippenberg (LPQM)

Clear memory

ClearAll["Global`*"]

Define potential

V[x_] := m w^2 x^2/2

Define stationary Smoluchowski equation

pde := -(1/(m g))D[-D[V[x],x]p[x],x]+d D[D[p[x],x],x] == 0

Find full solution

soln:= Evaluate[DSolve[{pde},p[x],x]]

soln

psol:= Evaluate[Simplify[p[x]/.soln[[1]]]]

psol

Try to solve for C[1] and C[2] so the constraints on p(x) are satisfied

constsol = Abort[]Solve[{Integrate[psol,{x,-Infinity,Infinity}]==1,psol>0},{C[1],C[2]},Reals]

(* ABORT THIS CALCULATION: Too complicated! *)

Try different approach by satisfying the constraints one at a time

constsol2 := Evaluate[Solve[Integrate[psol,{x,-Infinity,Infinity}]==1,{C[1],C[2]},Reals]]

constsol2

Abort[]FullSimplify[Solve[psol>0,{C[1],C[2]}],Reals]

(* ABORT THIS CALCULATION: Too complicated! *)

Plot[Erfi[x/Sqrt[2]],{x,-10,10}]

(* We see that no non-zero value of C[1] will satisfy the constraint p(x) > 0  *)

Define physical solution

truesol := psol/.constsol2[[1]]/.C[1]->0

truesol

Give values to variables and plot solution

m = 1;  (* Mass*)
w =1; (* Frequency *)
g = 10; (* Gravity *) 
d =1 ; (* Diffusion constant *)

truesol2 [y_] = Evaluate[truesol /. x-> y]

Plot[{truesol2[x],V[x]},{x,-10,10},PlotRange->{{-5,5},{0,0.5}}]