Numerical Boundary Conditions for Fokker Planck and Kramers equations

[rs12345]

Hi Everyone,

Apologies if this is posted in the wrong place.

I am trying to produce numerical solutions for the following initial value partial differential equations, namely the overdamped Schmoluchowski equation

$$\frac{\partial p(x,t)}{\partial t}=-\frac{\partial}{\partial x}\left(\frac{F(x,t)p(x,t)}{m\gamma}\right)+\frac{k_BT}{m\gamma}\frac{\partial^2p(x,t)}{\partial x^2}$$

and the associated full phase space Kramers equation

$$\frac{\partial p(x,v,t)}{\partial t}=-v\frac{\partial p(x,v,t)}{\partial x}+\frac{\partial}{\partial v}\left[p(x,v,t)\left(\gamma v -\frac{F(x,t)}{m}\right)\right]+\frac{k_BT(x)\gamma}{m}\frac{\partial^2p(x,v,t)}{\partial v^2}$$,

where $$p(x,t)$$ is a probability density function, $$\gamma$$ is the damping coefficient, $$x$$ the spatial coorindate, $$v$$ the velocity, $$F(x,t)$$ the force upon the particle, $$m$$ the mass and $$k_BT$$ Boltzmann's constant multiplied by temperature which we will consider homogenous in the overdamped case.

I wish, for the time being, just to produce qualitatively sensible answers from a simple technique and so wish to do so using simple explicit schemes eg FTCS using a real space discretisaton.

Now, I have observed the 'correct' behaviour for the overdamped case by implementing FTCS, that is a central space derivative approximation approach and in the full phase space case by using FTCS for the derivatives with respect to $$v$$ and using an upwind first order discretisation for derivatives with respect to $$x$$ to avoid the known instabilites in a pure advection term (albeit at the cost of introducing some additional spurious diffusion).

I do however, seem to be unable to treat the boundary conditions properly. Since $$p$$ is a probability density function it must be conserved. I believe I am right in saying that in order to affect this I require zero divergence and thus

$$\frac{\partial p}{\partial x}=\frac{\partial p}{\partial v}=0$$

on the boundaries. To do this on the boundaries I must set the relevant ghost points or boundary points to reflect this. Ie if there are N points on my grid ie would set $$p[N+1] = p[N-1]$$ for the central space approximations. This would presumably then need changing for the upwind discretisations to reflect the slightly different approximation.

However, when I do this, even for the simple overdamped case (with only FTCS approximations) my probability is not conserved and either 'leaks' out or explodes in magnitude. To see where it was going wrong I started with the simple one dimensional diffusion equation

$$\frac{\partial p(x,t)}{\partial t}=\frac{k_BT}{m\gamma}\frac{\partial^2p(x,t)}{\partial x^2}$$.

Doing so resulted in a conserved probability with just some small net probability flow at the start of the simulation. I observed similar behaviour when implementing the simplied form of the overdamped equation

$$\frac{\partial p(x,t)}{\partial t}=\alpha\frac{\partial p(x,t)}{\partial x}+\frac{k_BT}{m\gamma}\frac{\partial^2p(x,t)}{\partial x^2}$$.

However, when I utilise the full overdamped equation with a spatially dependent force, this doesn't appear to work. For example if considering a harmonic oscillator I have

$$\frac{\partial p(x,t)}{\partial t}=\frac{\partial}{\partial x}\left(\frac{kxp(x,t)}{m\gamma}\right)+\frac{k_BT}{m\gamma}\frac{\partial^2p(x,t)}{\partial x^2}$$

where $$k$$ is the spring constant. Why would this form not permit probability conservation using the technique I am using?

Could it be how I am discretising it? I have tried both

$$\frac{\partial}{\partial x}\left(xp(x,t)\right)\simeq \frac{x_{i+1}p(x_{i+1})-x_{i-1}p(x_{i-1})}{2\Delta x}$$

and

$$\frac{\partial}{\partial x}\left(xp(x,t)\right)\simeq p(x_i)+\frac{x_{i}(p(x_{i+1})-p(x_{i-1}))}{2\Delta x}$$.

Any help on this would be greatly appreciated along with any potential 'heads up' on potential issues when trying to make this all work in the full phase space situation with upwind discretisations or any other pointers in general!

Thankyou for reading!

Rich.

 

[jvicens]

Hi Rich,

Thank you for your post. It seems like you have a good understanding of the equations and their behavior, and are using appropriate numerical methods to solve them. The issue with probability conservation may be due to the discretization of the boundary conditions, as you have mentioned. It is important to carefully consider how the boundary conditions are being handled in order to ensure that they are accurately represented in the numerical solution. It may also be helpful to check for any errors in the implementation of the equations or the numerical method.

In terms of the full phase space Kramers equation, it is possible that the additional terms in the equation, such as the advection term and the spatially dependent force, may be contributing to the difficulty in conserving probability. It may be worth exploring different numerical methods or discretizations to see if that improves the conservation of probability.

In general, it is always a good idea to carefully check the boundary conditions and make sure they are accurately represented in the numerical solution. It may also be helpful to compare your results to analytical solutions or other numerical solutions to ensure that they are reasonable. I hope this helps and good luck with your research!