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$
14
  • 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
    Commented Feb 21, 2017 at 15:37
  • $\begingroup$ I think now it is ok. $\endgroup$
    – Agnieszka
    Commented Feb 21, 2017 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$ Commented Feb 21, 2017 at 15:44
  • $\begingroup$ I have approximated DiracDelta ( by function appro - see above), I still get zero everywhere... $\endgroup$
    – Agnieszka
    Commented Feb 22, 2017 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
    Commented Feb 23, 2017 at 11:07

0

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Browse other questions tagged or ask your own question.