0

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?

Agnieszka
  • 677
  • 4
  • 13
  • 1
    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] – yohbs Feb 21 '17 at 15:37
  • I think now it is ok. – Agnieszka Feb 21 '17 at 15:40
  • 3
    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. – Marius Ladegård Meyer Feb 21 '17 at 15:44
  • I have approximated DiracDelta ( by function appro - see above), I still get zero everywhere... – Agnieszka Feb 22 '17 at 09:07
  • 2
  • 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?)
    • – xzczd Feb 23 '17 at 11:07
  • I have changed the positon of MaxStep, chosen smaller k, included boundary conditions and chosen MaxStepSize->0.1. My Mathematica quits the kernel all the time when I start the calculation... – Agnieszka Feb 23 '17 at 14:01
  • You need to add @xzczd in your comment or I won't get a reminder. I suspect the RAM of your PC is eaten up. k = 5 still results in a quite sharp approximation of DiracDelta, to accurately represent it, a dense spatial grid is needed. (Just try With[{n = 25}, Plot[appro'[t], {t, -10, 10}, PlotRange -> All]~Show~ Plot[Interpolation[Table[{t, appro'[t]}, {t, -10, 10, 20/(n - 1)}], InterpolationOrder -> 4][t] // Evaluate, {t, -10, 10}, PlotRange -> All, PlotStyle -> Red]], you can adjust the n and observe. ) For a 5 dimentional problem, k = 5 may be too expensive. – xzczd Feb 24 '17 at 07:58
  • Also, notice MaxSteps and MaxStepSize is for adjustion of time step (the t in your case), while for your problem, the crucial part is the adjustion of spatial grid size. (p00, p11, a, b. ) – xzczd Feb 24 '17 at 08:06
  • ……And, your translation for $\frac{1}{2}\sum\limits_{i=1}^4 \sum\limits_{j=1}^4\partial^2_{x_i,x_j}[D_{ij}p]$ is wrong. You missed 2 +. A simple and clear way to deduce it with Mathematica is: var = {p00, p11, a, b}; sum2 = 1/2 Total@ Flatten@Table[ D[p[p00, p11, a, b, t] Dd[[i, j]], var[[i]], var[[j]]], {i, 4}, {j, 4}]; – xzczd Feb 24 '17 at 08:48
  • @xzczd , I have set k=1, added missing 2 +, set MaxSteps -> 10^3,MaxStepSize -> {1, 0.01, 0.01, 0.01, 0.01}, Mathematica still quits the kernel, also when I set smaller interval for calculating, p00,...,b – Agnieszka Feb 24 '17 at 09:23
  • By setting MaxStepSize -> {1, 0.1, 0.1, 0.1, 0.1} it seems that the calculation doesn't abort... I still wait for the output. – Agnieszka Feb 24 '17 at 10:17
  • Well, as mentioned above, MaxStepSize is an option for the adjustion of time step. To adjust the spatial grid size, use e.g. Method -> {"MethodOfLines", "SpatialDiscretization" -> {"TensorProductGrid", "MaxPoints" -> 25, "MinPoints" -> 25, "DifferenceOrder" -> 2}} (This defines a 2nd order spatial grid, with 25 points for each dimension. BTW suppose NDSolve uses 10^4 steps for time integration, you'll need 25^4 8/1024^3. 10^4 "GB" memory to finish the calculation.) BTW, I've define a mol function for this task in the link above. – xzczd Feb 24 '17 at 11:03
  • @xzczd , Although I have enough space on my computer (over 50GB), the computation is always aborted after some time and Mathematica quits the kernel. – Agnieszka Feb 28 '17 at 08:37
  • I think it would be good to first make sure the equation itself is correct. Can you talk a bit more about the definition of sigma and mu? – xzczd Feb 28 '17 at 09:02