#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);