0
$\begingroup$

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?

$\endgroup$
  • 1
    $\begingroup$ 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] $\endgroup$ – yohbs Feb 21 '17 at 15:37
  • $\begingroup$ I think now it is ok. $\endgroup$ – Agnieszka Feb 21 '17 at 15:40
  • 3
    $\begingroup$ 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. $\endgroup$ – Marius Ladegård Meyer Feb 21 '17 at 15:44
  • $\begingroup$ I have approximated DiracDelta ( by function appro - see above), I still get zero everywhere... $\endgroup$ – Agnieszka Feb 22 '17 at 9:07
  • 2
    $\begingroup$ 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?) $\endgroup$ – xzczd Feb 23 '17 at 11:07

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Browse other questions tagged or ask your own question.