# Fokker-Planck equation in Mathematica

    ClearAll["Global*"]
α1 = 1;
σ1 = 0.1;
α2 = 1.2;
β = 0.4;
tf = 5000;
V = 0.00001;
proc = ItoProcess[{\[DifferentialD]R[
t] == (-1/R[t]^2*(α1 - α2*(Cos[ψ[t]])^2) +
1/R[t]^3*(2*α1*σ1*
Cos[χ[t] + ψ[t]] - α2*σ1*
Cos[ψ[
t]]*(3*Cos[χ[t] + ψ[t]]*Cos[ψ[t]] -
Cos[χ[t]])) -
V*Cos[χ[t] + ψ[
t]]) \[DifferentialD]t + \[DifferentialD]w1[
t], \[DifferentialD]ψ[
t] == (β/R[t]^2*Sin[ψ[t]] - β*σ1/
R[t]^3*(3*Cos[χ[t] + ψ[t]]*Sin[ψ[t]] -
Sin[χ[t]]) - α2/
R[t]^3*(Sin[ψ[t]]*Cos[ψ[t]]) -
1/R[t]^4*(-α1*σ1*
Sin[χ[t] + ψ[t]] - α2*σ1*
Sin[ψ[
t]]*(3 Cos[χ[t] + ψ[t]]*Cos[ψ[t]] -
Cos[χ[t] + ψ[t]])) -
V/R[t]*Sin[χ[
t]]) \[DifferentialD]t + \[DifferentialD]w2[
t], \[DifferentialD] χ[
t] ==  (-β/R[t]^2*Sin[ψ[t]] - β*σ1/
R[t]^3*(3*Cos[χ[t] + ψ[t]]*Sin[ψ[t]] -
Sin[χ[t]])) \[DifferentialD]t + \[DifferentialD]w3[
t]}, {R[t], ψ[t]}, {{R, ψ, χ}, {1, 1, 0}}, {t,
0}, {w1 \[Distributed] WienerProcess[],
w2 \[Distributed] WienerProcess[],
w3 \[Distributed] WienerProcess[]}];

K1 = 1;
K2 = 1;
K3 = 1;
pde = D[P[R, ψ, χ, t],
t] == -1/R*
D[(R*(-1/R^2*(α1 - α2*(Cos[ψ])^2) +
1/R^3*(2*α1*σ1*
Cos[χ + ψ] - α2*σ1*
Cos[ψ]*(3*Cos[χ + ψ]*Cos[ψ] -
Cos[χ])) - V*Cos[χ + ψ])*
P[R, ψ, χ, t]), R] -
D[(β/R^2*Sin[ψ] - β*σ1/
R^3*(3*Cos[χ + ψ]*Sin[ψ] -
Sin[χ]) - α2/R^3*(Sin[ψ]*Cos[ψ]) -
1/R^4*(-α1*σ1*
Sin[χ + ψ] - α2*σ1*
Sin[ψ]*(3 Cos[χ + ψ]*Cos[ψ] -
Cos[χ + ψ])) - V/R*Sin[χ])*
P[R, ψ, χ, t], ψ] -
D[ (-β/R^2*Sin[ψ] - β*σ1/
R^3*(3*Cos[χ + ψ]*Sin[ψ] - Sin[χ]))*
P[R, ψ, χ, t], χ] +
K1*1/R*D[D[R*P[R, ψ, χ, t], R], R] +
K2*D[D[P[R, ψ, χ, t], ψ], ψ] +
K3*D[D[P[R, ψ, χ, t], χ], χ]


I have a problem here where I've three coupled stochastic equations being modeled using Ito Process under the variable proc. I'm able to solve the equation. Now, we use Fokker Planck to convert a stochastic equation into a deterministic probability evolution equation. However, I'm unable to solve the PDE part and get a PDF for it.

• If you already have the trajectories from the langevin equations you can reconstruct the PDF at each timestep (maybe with some statistical fluctuatitions) Why would you need to solve the associated fokker plank (which is harder to solve)? – SSC Napoli Jun 17 '15 at 14:37
• How did you try solving the PDE? NDSolve should do it, with friendly initial conditions (I don't think it will accept delta functions). – Ian Jun 18 '15 at 2:15
• Yes, but perhaps my execution is not correct. Here are the conditions i'd like the code to work with: P(\chi)=P(\chi+2\pi) P(\psi)=P(\psi+2\pi) for all R and P(R going to infinity for all \chi and \psi)=0. – Skidmore_TCIS Jun 18 '15 at 6:33

This isn't a working answer, but maybe it will help.

I've guessed some about your auxiliary conditions too, and maybe that's the reason I'm not getting a solution. The error the following code gives me is:

NDSolve::ndnum: Encountered non-numerical value for a derivative at t == 0..


Naturally, this occurs at R=0 in your equation. Perhaps the support of P is not supposed to include R=0? The output of RandomFunction for the Ito SDE you give takes wild jumps if R gets near zero too.

Anyway, here's how I attempted a solution with NDSolve:

T = 1;
inf = 10;
\[Epsilon] = 10^-3;
P0 = MultinormalDistribution[{1, 1, 0}, \[Epsilon] IdentityMatrix[3]];
NDSolve[{pde,
P[R, \[Psi], \[Chi], 0] == PDF[P0, {R, \[Psi], \[Chi]}],
P[R, -\[Pi], \[Chi], t] == P[R, \[Pi], \[Chi], t],
P[R, \[Psi], -\[Pi], t] == P[R, \[Psi], \[Pi], t],
P[inf, \[Psi], \[Chi], t] == 0, P[-inf, \[Psi], \[Chi], t] == 0
}, P, {R, -inf, inf}, {\[Psi], -\[Pi], \[Pi]}, {\[Chi], -\[Pi], \[Pi]}, {t, 0, T}]

• Thank you for your help! I added this piece to the existing code and it's giving a bunch of errors besides the one you mentioned above, mostly that 1/0^2 is encountered. Still working through it! :) – Skidmore_TCIS Jun 19 '15 at 12:09
• @Skidmore If you get a working solution, it would be great if you could share it here as a self-answer. – Ian Jun 19 '15 at 12:13