# Solution of Fokker-Planck equation with delta initiall condition - everywhere zero

I try to solve a Fokker-Planck equation (which I derived from corresponding Ito equation)

$\partial_tp=-\sum\limits_{i=1}^4 \partial_{x_i}(\mu_i p)+\frac{1}{2}\sum\limits_{i=1}^4 \sum\limits_{j=1}^4\partial^2_{x_i,x_j}[D_{ij}p], D_{ij}=\sum\limits_{k=1}^2\sigma_{ik}\sigma_{jk}$,

my code is:

appro = With[{k = 1}, Erf[k #]/2 + 1/2 &]

sigma := 1/
Sqrt[2]*{{Sqrt[gamma/n]*(n + 1)*2*
a, (Sqrt[gamma/n]*(n + 1)*2*b)}, {-Sqrt[gamma*n]*2*
a, -(Sqrt[gamma*n]*2*b)},
{-Sqrt[gamma*n]*p00 + Sqrt[gamma/n]*(n + 1)*p11,
I (-Sqrt[gamma*n]*p00 + Sqrt[gamma/n]*(n + 1)*p11)},
{-Sqrt[gamma*n]*p00 +
Sqrt[gamma/n]*(n + 1)*
p11, -I (-Sqrt[gamma*n]*p00 + Sqrt[gamma/n]*(n + 1)*p11)}};

Dd = Table[0, {i, 4}, {j, 4}];

Do[Dd[[i, j]] = Sum[sigma[[i, k]]*sigma[[j, k]], {k, 1, 2}], {i, 1,
4}, {j, 1, 4}]

mu = {-gamma*n*p00 + gamma*(n + 1)*p11,
gamma*n*p00 - gamma*(n + 1)*p11, (2*dn + d)*b -
gamma*(2*n + 1)*a, -(2*dn + d)*a - gamma*(2*n + 1) b}

gamma = 1;
d = 1;
dn = 1;
n = 1

psol = First[p /. NDSolve[{D[p[p00, p11, a, b, t], t] +
D[mu[[1]]*p[p00, p11, a, b, t], p00] +
D[mu[[2]]*p[p00, p11, a, b, t], p11] +
D[mu[[3]]*p[p00, p11, a, b, t], a] +
D[mu[[4]]*p[p00, p11, a, b, t], b] -
0.5*(D[D[Dd[[1, 1]]*p[p00, p11, a, b, t], p00] +
D[Dd[[2, 1]]*p[p00, p11, a, b, t], p11] +
D[Dd[[3, 1]]*p[p00, p11, a, b, t], a] +
D[Dd[[4, 1]]*p[p00, p11, a, b, t], b], p00] +
D[D[Dd[[1, 2]]*p[p00, p11, a, b, t], p00] +
D[Dd[[2, 2]]*p[p00, p11, a, b, t], p11] +
D[Dd[[3, 2]]*p[p00, p11, a, b, t], a] +
D[Dd[[4, 2]]*p[p00, p11, a, b, t], b], p11] +D[
D[Dd[[1, 3]]*p[p00, p11, a, b, t], p00] +
D[Dd[[2, 3]]*p[p00, p11, a, b, t], p11] +
D[Dd[[3, 3]]*p[p00, p11, a, b, t], a] +
D[Dd[[4, 3]]*p[p00, p11, a, b, t], b], a]+ D[
D[Dd[[1, 4]]*p[p00, p11, a, b, t], p00] +
D[Dd[[2, 4]]*p[p00, p11, a, b, t], p11] +
D[Dd[[3, 4]]*p[p00, p11, a, b, t], a] +
D[Dd[[4, 4]]*p[p00, p11, a, b, t], b], b]) == 0,
p[p00, p11, a, b, 0] ==
appro'[p00 - 0.3]*appro'[p11 - 0.7]*appro'[a - 0.2]*
appro'[b - 0.1],
p[10, p11, a, b, t] == p[-10, p11, a, b, t] ==
p[p00, 10, a, b, t] == p[p00, -10, a, b, t] ==
p[p00, p11, 10, b, t] == p[p00, p11, -10, b, t] ==
p[p00, p11, a, 10, t] == p[p00, p11, a, -10, t] == 0},
p, {t, 0, 10}, {p00, -10, 10}, {p11, -10, 10}, {a, -10,
10}, {b, -10, 10}, Method -> {"MethodOfLines",


"SpatialDiscretization" -> {"TensorProductGrid", "MaxPoints" -> 25, "MinPoints" -> 25, "DifferenceOrder" -> 2}}]]

My Mathemtica quits the kernel, when I start the calculation.

What I do wrong?

• Your code returns an error. The best thing you can do in this situation is to write a simpler version of your equation just to see that your syntax is right (a minimal working example). It is currently not. For example, your Dd variable definition does not work. I think a better way to achieve what you want is to write Dd = sigma.Transpose[sigma] Commented Feb 21, 2017 at 15:37
• I think now it is ok. Commented Feb 21, 2017 at 15:40
• I don't think a numeric solver like NDSolve will be able cope with a DiracDelta, it will just see it as something that is zero everywhere. The DiracDelta is for symbolic calculations. Commented Feb 21, 2017 at 15:44
• I have approximated DiracDelta ( by function appro - see above), I still get zero everywhere... Commented Feb 22, 2017 at 9:07
• 1. The position of MaxSteps is wrong. 2. The bcart warning is a more serious problem, check this post for more information. You need to add proper boundary condition. 3. The approximate DiracDelta is probably too slender, try a smaller k, a denser spatial grid may be also necessary. (I think you've already read this post?) Commented Feb 23, 2017 at 11:07