I want to solve the initial value problem of a particle diffusing in a well on a disk.

Paticle density$\rho(\vec{r},t)$satisfy the following equation: $$ \partial\rho/\partial t=\nabla\cdot(\nabla V(r)\rho)+D_t\nabla_{x,y} ^2\rho\\ \frac{\partial\rho}{\partial\hat{n}}=0,\text{on the boundary } (|\vec{r}|==1)\\ \rho(r,0)=\text{Gaussian}(r) $$ $V(r)$ is a potential well, $V(r)=r^2$.

I try to solve this problem in Mathematica with the Finite element method.

But the integral of the particle density does not conserve (there appears flux on the boundary).

It seems to be due to my false boundary condition.

How should I set the boundary condition?


1 Answer 1


The diffusion equation you are using can written as a balance law, a simple continuity equation $$ \partial_t \rho + \nabla \cdot \vec{j} = 0\,, $$ using a current density $$ \vec{j} = - \nabla V \rho - D_t \nabla \rho\,. $$ For the mass inside a volume you integrate the differential balance law to obtain an integral balance law $$ \partial_t \int_\Omega \rho \text{d}V = - \int_\Omega \nabla \cdot \vec{j} \text{d}V = - \int_{\partial\Omega} \vec{j} \cdot \vec{n} \text{d}S\,. $$ The very basics of continuum mechanics... ;) In words: The rate of change of the mass/probability in $\Omega$ (left) equals the outflux (right) that is computed by integrating the current density normal to the surface.

To ensure conservation in your problem, you thus need to make sure that the integral on the right vanishes, for example by making $\vec{j}$ vanish everywhere on the boundary.

  • $\begingroup$ Thanks for your help. I've set the boundary condition to $\hat{n}*\vec{j}=0$ (a robin boundary condition). But my code still doesn't work. Is there anything wrong related to the Finite Element Method? $\endgroup$
    – 江蛮子
    Jul 7, 2022 at 2:25
  • $\begingroup$ Not generally, or not that I could tell without further info, no. Be aware that it is important to make a right discretisation of a problem, otherwise you can have growing errors crashing your simulation. To get started on that see, e.g., W.A. Strauss, PDE, Chapter 8 (8.2 for diffusion). If it's still not working, maybe ask a new question, maybe on mathematica.stackexchange.com Or have a look here: reference.wolfram.com/language/FEMDocumentation/tutorial/… Good luck troubleshooting! The physics should be fine now, at least. $\endgroup$
    – kricheli
    Jul 7, 2022 at 5:35

Your Answer

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

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