- #1
rs12345
- 1
- 0
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
[tex]\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}[/tex]
and the associated full phase space Kramers equation
[tex]\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}[/tex],
where [tex]p(x,t)[/tex] is a probability density function, [tex]\gamma[/tex] is the damping coefficient, [tex]x[/tex] the spatial coorindate, [tex]v[/tex] the velocity, [tex]F(x,t)[/tex] the force upon the particle, [tex]m[/tex] the mass and [tex]k_BT[/tex] 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 [tex]v[/tex] and using an upwind first order discretisation for derivatives with respect to [tex]x[/tex] 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 [tex]p[/tex] 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
[tex]\frac{\partial p}{\partial x}=\frac{\partial p}{\partial v}=0[/tex]
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 [tex]p[N+1] = p[N-1][/tex] 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
[tex]\frac{\partial p(x,t)}{\partial t}=\frac{k_BT}{m\gamma}\frac{\partial^2p(x,t)}{\partial x^2}[/tex].
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
[tex]\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}[/tex].
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
[tex]\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}[/tex]
where [tex] k[/tex] 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
[tex]\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}[/tex]
and
[tex]\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}[/tex].
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.
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
[tex]\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}[/tex]
and the associated full phase space Kramers equation
[tex]\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}[/tex],
where [tex]p(x,t)[/tex] is a probability density function, [tex]\gamma[/tex] is the damping coefficient, [tex]x[/tex] the spatial coorindate, [tex]v[/tex] the velocity, [tex]F(x,t)[/tex] the force upon the particle, [tex]m[/tex] the mass and [tex]k_BT[/tex] 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 [tex]v[/tex] and using an upwind first order discretisation for derivatives with respect to [tex]x[/tex] 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 [tex]p[/tex] 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
[tex]\frac{\partial p}{\partial x}=\frac{\partial p}{\partial v}=0[/tex]
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 [tex]p[N+1] = p[N-1][/tex] 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
[tex]\frac{\partial p(x,t)}{\partial t}=\frac{k_BT}{m\gamma}\frac{\partial^2p(x,t)}{\partial x^2}[/tex].
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
[tex]\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}[/tex].
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
[tex]\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}[/tex]
where [tex] k[/tex] 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
[tex]\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}[/tex]
and
[tex]\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}[/tex].
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.