How can one simulate diffusion of a particle which has a natural axis and diffusion coefficient of the particle in the direction of the axis is D_1 and in perpendicular to the axis is D_2? Could anyone provide a code?
-
1$\begingroup$ Diffusion equation where the diffusion constant is a tensor instead of a scalar should work. $\endgroup$ – alarge Jun 23 '17 at 4:11
If I am reading you correctly, recall that the diffusion operator is just
$$ D_1 u_{yy} + D_2u_{xx} = D_1 \frac{u_{i+1,j} - 2u_{i,j} + u_{i-1,j}}{\Delta y^2} + D_2 \frac{u_{i,j+1} - 2u_{i,j} + u_{i,j-1}}{\Delta x^2} $$
to second order. That should give you an idea for how to write the loop.
But there's another way to do this. If $A$ is the diffusion operator which is a Strang matrix (-2 on the diagonal and 1 on the off diagonals), then notice that if you make $u$ a matrix where $u_{i,j}(t) = u(x_j,y_i,t)$, i.e. the each column is a separate $x$ and each row is a separate $y$ (notice the indexing swap there), then this is very easy to code up in a vectorized form. Notice that $Au$ uses the columns of $u$, so this is diffusing along $y$. Meanwhile, $uA$ is diffusing along the $x$. Thus, multiplying by $\frac{D_z}{\Delta z^2}$ for $z=x,y$ respectively,
$$ A_yu + uA_x $$
is an easy vectorized way to perform the diffusion. What's nice about this is that BLAS is multithreaded, so this will auto-parallelize in many languages like Julia, MATLAB, and Python.
-
$\begingroup$ This assumes that the axis is always aligned along x and y, which is not necessarily reasonable. e.g., an anisotropic molecule will also diffuse rotationally. $\endgroup$ – AJK Nov 20 '17 at 6:32
-
$\begingroup$ If that's the case just use polar coordinates and discretize the LaPlacian there? There's no issue... $\endgroup$ – Chris Rackauckas Jan 19 '18 at 12:57
If you want to simulate the (stochastic) Brownian motion of an object with a diffusion coefficient that depends on the axis, this is discussed in detail in the paper:
Brownian Motion of a Particle of General Shape in Newtonian Fluid Makino and Doi, J. Phys. Soc. Jpn. 73, pp. 2739-2745 (2004) https://doi.org/10.1143/JPSJ.73.2739
The critical issue is that rotational diffusion will alter what the axes are, leading to very different behavior.
It should also be possible to write and solve a Fokker-Planck equation for the probability distribution of the orientation variables and object position. This is (I think) discussed in Doi + Edwards's book The Theory of Polymer Dynamics.
The problem here amounts to apply the random force in the body reference frame of the particle. There, you can set the friction and/or diffusion coefficients separately. Once the forces are computed, you can rotate them back to the laboratory reference frame. A book such as "Dynamics of colloids" will provide you with all the necessary background. https://www.amazon.fr/Introduction-Dynamics-Colloids-J-K-G-Dhont-ebook/dp/B00DTSEAY2/ref=sr_1_1?ie=UTF8&qid=1495569759&sr=8-1&keywords=dynamics+of+colloids