Featured Post

#75 The application of Fokker-Planck equation-II 福克-普朗克方程的應用-II

Today we are going to talk about the first passage time problem. It is important in physical, biomedical science and economy.

In first passage time problem, we are interested in the quantity T(x; x0), the time it takes a particle to reach a point x from x0. However, due to the stochastic nature of the system, T(x; x0) is associated with a probability f(x,t), which denotes the probability that x is reached at time t. We need to solve the Fokker-Planck equation for P(x,t;x0,t0). We consider only 1D system for simplicity.
The boundary condition of the system looks like this:
A poor drunk (so he can't recognize the cliff) random walker is walking on a 1D line with a dangerous cliff at x = x_{end} and reflecting boundary condition at x=0.


We denote the probability that the dangerous cliff is not reached at time t as S(x0, t), the survival probability. From basic probability theory, the following 3 equation is obtainable:
The function f is difficult to obtain. However, the moment of arrival time T is much more easier. From integration by parts, we could write:


Now assume that our drunk random walker is move under a constant biased force f toward the negative side (perhaps we hope the drunk random walker to live longer!), the corresponding Fokker-Planck equation is:
The quantity we are interested in is
Considering the integration domain, an easier way to obtain the answer is to transform our original Fokker-Planck equation into its Kolmogorov backward equation:

Integrate the Kolmogorov backward equation with respect to x domain and t domain give rise to the following:
This ODE is not difficult to solve. The solution becomes
 
A simple computer simulation program could verify the result:
/*------MATLAB code------*/


%% basic setting
r = 5E-9; eta = 8.90E-4; zeta = 6*pi*eta*r;
kB = 1.38E-23;
syms x;
dt = 1E-10;
n_sim = 5E6;

T =1800; A = sqrt(2*zeta*kB*T*dt);
D = kB*T/zeta;
f = 1E-11;
x_end = 9E-9;
%%
n_round = 1000;
total_data = zeros(1,8);
for x0 = 1E-9:1E-9:8E-9
    time_get = zeros(1, n_round);
    for j = 1:1:n_round
        x_total = zeros(1, n_sim);
        x_total(1)=x0;
        for i = 2:1:n_sim
            x_old = x_total(i-1);
            x_total(i) = x_old + 1/zeta*(-dt*f+A*randn(1,1));
            if x_total(i) <0
                x_total(i) = - x_total(i);
            end
            if x_total(i)>=x_end
                time_get(j) = i*dt;
                break
            end
        end
        display(j);
    end
    total_data(int8(x0/1E-9))=mean(time_get);
end
%theoretical
x0 = 1E-9:1E-9:8E-9
tau_theo = zeta^2*D/f^2*(exp(f*x_end/(zeta*D))-exp(f*x0/(zeta*D))+zeta/f*(x0 - x_end));
figure
scatter(x0, total_data);hold on;
plot(x0, tau_theo,'LineWidth',2);
xlabel('x_0','FontWeight','bold','FontSize',14);
ylabel('\tau','FontWeight','bold','FontSize',14);
set(gca,'FontSize',14);

Comments