2
$\begingroup$

enter image description hereI'm solving a Fokker-Planck equation with Python by solving AB=C where A is coefficient matrix, B is the vector of unknowns, and C is the right-hand side. The Fokker-Planck is below:

$$ \frac{∂ρ(x,y,t)}{∂t} ​= -\beta \delta \left[\left(\frac{∂}{∂x​}F_x(x,y)ρ(x,y,t) + \frac{∂}{∂y​}F_y(x,y)ρ(x,y,t) \right) \\ + \delta \left(\frac{∂^2ρ(x,y,t)}{∂x^2} + \frac{∂^2ρ(x,y,t)}{∂y^2} \right)\right] \quad (1) $$

The challenge is that the solution B seems unstable since it appears 'ripples' during the simulation and this is an unexpected behaviour. I've tried to reduce the time step or increase the grid size but none of these helps. However, if I just make the $F$ smaller by scaling it by a small factor then the 'ripples' disappear. What would have caused this issue and how to sort this out please?

c = scipy.sparse.linalg.spsolve(A, b)
p[:,:] = c.reshape(Ny+1,Nx+1).T
# Normalize p
p /= np.sum(p)
$\endgroup$
6
  • 2
    $\begingroup$ This question would be more appropriate for stackoverflow $\endgroup$ – spiridon_the_sun_rotator Apr 29 at 7:04
  • $\begingroup$ I posted this question on stackoverflow but they asked me to post here $\endgroup$ – An Nguyen Apr 29 at 7:09
  • $\begingroup$ What scheme do you use for the first derivatives in space? $\endgroup$ – Pavlo. B. Apr 29 at 7:26
  • $\begingroup$ I used Forward Euler scheme $\endgroup$ – An Nguyen Apr 29 at 7:30
  • 1
    $\begingroup$ I’m voting to close this question because it is about debugging code, not physics. $\endgroup$ – Jon Custer May 4 at 14:33
1
$\begingroup$

This is a common and known problem - an immediate artefact of discretization and the straightforward approach to solving that you use.

  • If you want to do it yourself, you have to go for more sophisticated algorithms, notably for more sophisticated approximations for the derivatives. I suggest consulting the Numerical recipies, although there might be also something in the Reif's book on FPE.
  • Alternative approach is to have it solved for you by a numerical pde solver (which are available for python).
  • Finally, one can do 50/50, e.g., by separating variables and solving the coordinate part using an existing solver - they are more readily available and stable for the elliptics (Laplace/Poisson-like) equations than for the parabolic equations equations (such as FPE or diffusion equation).
$\endgroup$
2
  • $\begingroup$ Thanks, I would need to do more research on this then. I guess it would probably need more sophisticated grid mesh for the integration to be stable $\endgroup$ – An Nguyen Apr 29 at 9:05
  • $\begingroup$ Representing derivatives by first order differences is usually a bad idea - the error is too big. $\endgroup$ – Roger Vadim Apr 29 at 9:07

Not the answer you're looking for? Browse other questions tagged or ask your own question.