0
$\begingroup$

I want to simulate Brownian motion in 3-D for the following conditions: $$p(x=0,y=0,z=0,t=0)=1$$ $$p(x,y,z=c,t)=0$$ where $p$ is the probability of finding molecules in the 3-D environment. I want to find the number of particles getting absorbed on a certain region on the absorbing wall in 5 secs. As can be inferred the absorbing boundary is at $z=c$, and I want the probability of molecules hitting a circular area centered at $(0,0,2.5\times10^{-4})$ with a radius $5\times10^{-6}$m. The fluid environment has drift,i.e. $u,v,w$ in three directions. The particles in $x,y,z$ directions are moving as \begin{equation} \text{d}X_j(t)=u\text{d}t+\sqrt{2D}\text{d}B_j(t),\\ \text{d}Y_j(t)=v\text{d}t+\sqrt{2D}\text{d}B_j(t),\\ \text{d}Z_j(t)=w\text{d}t+\sqrt{2D}\text{d}B_j(t), \end{equation} where $\text{d}B$ is the derivative of Brownian process which I am simulating through normal random varaible.

I am using following MATLAB code to simulate particle diffusion:

close all; 
clear all;
clc;
TotalNumberReceivedUp=0;
Ts=5;
TotalSimulations1=100000;
flg=0;
for particlesNum=1:TotalSimulations1
flg=flg+1
D=4*10^-9; 
c=250*10^-6;u=5*10^-6;v=2*10^-6;w=3*10^-6; 
e=5*10^-6;
X=zeros;Y=zeros;Z=zeros;
X(1)=0;Y(1)=0;Z(1)=0;
j=1;
for i=5*10^-4:5*10^-4:Ts
t2=i;    
t1=t2-5*10^-4;
j=j+1;
r = normrnd(0,1);
q= normrnd(0,1);
p= normrnd(0,1);
X(j)=X(j-1)+(sqrt(2*D))*sqrt(t2-t1)*p+u*(t2-t1);
Y(j)=Y(j-1)+(sqrt(2*D))*sqrt(t2-t1)*q+v*(t2-t1);
Z(j)=Z(j-1)+(sqrt(2*D))*sqrt(t2-t1)*r+w*(t2-t1);
if ((Z(j)>=c) && ...
    (sqrt((X(j))^2+(Y(j))^2)<=e))
TotalNumberReceivedUp=TotalNumberReceivedUp+1;
break
end
end
end
ProbabilityUp=TotalNumberReceivedUp/TotalSimulations1

But my results are not coming good, the probable particles hitting are more than those I am getting through analysis. I think it may be because of the excess over the boundary problem which has been discussed in "A comparison of four methods for simulating the diffusion process" for 1-D movement with two absorbing barriers. I have tried to use the same factor in 3-D but it doesn't work well.

$\endgroup$
  • $\begingroup$ What happens if you decrease the time step length ($t2 - t1$)? Does it become closer to the analytical result? Also, your statement is ambiguous: what is the absorbing boundary? Is it the whole plane $z = c$ or just the circle of radius $\sqrt{e}$ within that plane? $\endgroup$ – Juan M. Bello-Rivas Jan 1 at 18:26
  • $\begingroup$ @Juan M. Bello-Rivas Actually the whole wall is absorbing but my problem is to find the number of molecules adsorbed on a circular region of radius $e$, which is on the wall. $\endgroup$ – Userhanu Jan 2 at 5:44
  • $\begingroup$ @Juan M. Bello-Rivas Earlier I was using step length $5\times10^{-2}$ and now $5\times10^{-4}$, but the answer is more or less same!!! $\endgroup$ – Userhanu Jan 2 at 5:46
1
$\begingroup$

If the whole wall is absorbing, then you should replace:

    if ((Z(j)>=c) && ...
        (sqrt((X(j))^2+(Y(j))^2)<=e))
        TotalNumberReceivedUp=TotalNumberReceivedUp+1;
        break
    end

by

    if Z(j) >= c
        if sqrt(X(j)^2 + Y(j)^2) <= e
            TotalNumberReceivedUp = TotalNumberReceivedUp + 1;
        end

        break
    end

Otherwise, the current absorbing boundary is the cylinder $\{ (x, y, z) \in \mathbb{R}^3 \, | \, x^2 + y^2 = e^2, z \ge c \}$.

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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