Let us discretize space as
z_i = i \; \Delta z
and time as
t_n = n \; \Delta t.
Let \mu be the fraction of density transfered from site i to site i+1 during one time step, so that the corresponding current would be
J_{i,i\!+\!1} = \mu \; \rho_i
and let us denote the density in site i at time n as
\rho^{n}_i = \rho(z_i,t_n).
Except for the most left (i=0) and the most right (i=i_mx) sites (where boundary conditions have to applied properly), the discretized diffusion equation reads
\rho^{n+1}_i = (1-2\mu) \rho^{n}_i + \mu (\rho^{n}_{i-1} + \rho^{n}_{i+1}).
Note that the physical diffusion constant is related to the discretization parameters as
D = \frac{\mu \Delta z^2}{\Delta t}.
To avoid numerical instabilities, one can set \mu=0.1 and choose a proper time step.
We shall later apply the above method to the DDD model. In this case, \mu is locally replaced by
\mu \rightarrow \mu_i = \mu(\rho^{n}_i).
As a simple test of the method, we start with ordinary diffusion in a box of width 100 with hard walls at z=0 and z=100. The initial density at t=0 is a narrow Gaussian, centered in the middle of the box. The diffusion constant was set to 10.
In a semilog plot, one gets the expected evolution of density:
The variance of z-distributions follows initially the behaviour of free diffusion. Later it saturates because of the confinement in the box:
Next we turn to the model with density dependent diffusion constant (DDD model). The initial density is now localized in the z=0 layer. The evolution of the particle densities, as resulting from numerical integration of the FP equation, looks as follows:
In the semilog. plot one can see the characteristic bi-phasic behaviour and the almost exponential regime at intermediate z.
Here the corresponding cummulative probabilities (= invasion profiles):
The variance as a function of time ...
... starts sub-diffusively and later appears to become slightly super-diffusive.
No comments:
Post a Comment