# Simulating Brownian motion in 3-D for first hitting time?

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 $$$$\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),$$$$ 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;
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))
break
end
end
end


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.

• 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? – Juan M. Bello-Rivas Jan 1 at 18:26
• @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. – Ankit Jan 2 at 5:44
• @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!!! – Ankit Jan 2 at 5:46

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

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


by

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

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