Multiplying Eq. (7.69) with the integrating factor eR nxdx leads to
−(xny′)′ =−s
Dxny. (7.74)
This is the canonical Sturm-Liouville form [101]
−(py′)′+qy =λwy, (7.75)
with p(x) = xn, the weighting function w(x) = xn, q(x) ≡ 0, and the spectral parameter λ=−s/D.
We now observe that u := y −1 transforms the homogeneous problem (7.69) with inhomogeneous boundary conditions (7.70) into an inhomogeneous problem with homogeneous boundary, conditions
−(pu′)′ =λwu+λw (7.76)
Au(0) +Bu(a) = 0, (7.77)
where u(x) = (u(x), u′(x))T, and the two possible choices of A and c correspond to the Dirichlet problem and the Dirichlet-Neumann problem respectively. This is easier to solve, since it determines a self-adjoint operatorL defined by
Lu= 1
w[−(pu′)′] (7.78)
in the weighted Hilbert spaceH =L2(J, w), where we have defined the open inter-valJ = (0, b). This operator is not to be confused with the Laplace transformation operator in Eq. (7.55), which can be recognized from the index indicating the trans-formed variable.
This can be seen as follows: Let u, v ∈ H. Then the inner product is given by hu, vi = Rb
0 uvw dx; taking into account that¯ u and v satisfy the homogeneous boundary conditions (7.77), we get after integrating twice by parts
hu,Lvi = − Z b
0
¯
u(pv′)′dx
= h
p(vu¯′−uv¯ ′)ib 0−
Z b 0
(p¯u′)′v dx
= hLu, vi. (7.79)
Using the definition from Eq. (7.78) the boundary value problem given by Eqs. (7.69–
7.70) can be simplified to
(L −λ1)u = λ, (7.80)
Au(0) +Bu(b) = 0. (7.81)
7.6 Formal solution of the boundary value
has nontrivial solutions uk with eigenvalues αk, k= 1,2, . . .
Luk=αkuk. (7.83)
Because of the self-adjointness of L, the eigenvalues αk are real and the eigen-functions uk form an orthonormal basis of H. Furthermore αk > 0 holds, since αk =huk,Luki.
Hence the solution u of the inhomogeneous problem given by Eqs. (7.80) and (7.81) can be expressed as an expansion in this basis,
u=
∞
X
k=1
ckuk, (7.84)
with ck =huk, ui. The coefficients ck can be derived from Eq. (7.80):
huk,Lui − huk, λui=huk, λi. (7.85) Again, making use of the definition of a self-adjoint operator, we can pullLinto the first component of the inner product. Employing Eq. (7.83) we get
ck = huk, λi
αk−λ. (7.86)
The solution of the inhomogeneous problem reads u=
∞
X
k=1
huk, λi
αk−λuk, (7.87)
and hence
y= 1 +
∞
X
k=1
huk, λi
αk−λuk. (7.88)
The eigenfunctions ukdo not depend on λ=−s/D, and making use of the linearity of the inverse Laplace transformation one gets
y(t, x) = L−1s [y(s, x)](t)
= L−1s [1] +
∞
X
k=1
huk,1iukL−1s
λ αk−λ
(7.89)
= δ(t)+
∞
X
k=1
huk,1iuk(αkDe−αkDt−δ(t)).
This can be further simplified. Using Eq. (7.84) we know that P∞
k=1huk,1iuk = 1 and hence the two delta functions cancel out. Then we have
y(t, x) =
∞
X
k=1
huk,1iukαkDe−αkDt. (7.90) This function is normalized because, knowing thatαk is positive,
Z ∞ 0
y(t, x)dt =
∞
X
k=1
huk,1iukαkD Z ∞
0
e−αkDtdt
=
∞
X
k=1
huk,1iuk= 1. (7.91)
We are now able to solve the specific boundary value problems for the three different kinds of boundaries at zero.
7.6.1 Entrance boundary
As we know from the classification of the origin, we have an entrance boundary for n≥1 orν ≤0. The general solution of the homogeneous differential equation (7.82) is
u(x) =AxνJν √ αx
+BxνYν √ αx
, (7.92)
where Jν and Yν are the Bessel functions of the first and second kind respectively.
The derivative of the general solution reads u′(x) =xν√
α
AJν−1 √ αx
+BYν−1 √ αx
. (7.93)
In the limit x → 0 and for ν ∈ Z the term containing the Bessel functions of the second kind diverges whereas the one involving the Bessel functions of the first kind is zero; the boundary condition limx→0u′(x) = 0 can only be fulfilled if B = 0. The second boundary condition u(b) = 0 gives us the eigenvalues of the problem by the requirement that Jν(√αkb) = 0. Calling the ith zero of the Bessel function of the first kindji, one can use the orthogonality relation for Bessel functions [96]
Z b 0
Jν
jm
x b
Jν
jn
x b
x dx= 1
2b2Jν+12 (jm)δmn (7.94) to compute the constantA. Identifyingji with√αib, the normalized eigenfunctions of the problem read
uk(x) =
√2bν−1 xbν
Jν jkx b
Jν+1(jk) . (7.95)
For non-integer ν the Bessel functions Jν and J−ν are two linear independent solutions of Eq. (7.82), and writing the general solution as the linear combination
u(x) = AxνJν √ αx
+BxνJ−ν √ αx
(7.96) instead of using Eq. (7.92) makes the analysis easier. The derivative reads
u′(x) =xν√ α
AJν−1(√
αx)−BJ1−ν(√ αx)
. (7.97)
In this case the function xνJν−1(√
αx) diverges as we approach the origin, and hence we have to setA= 0 in order to fulfill the left boundary condition. The right boundary condition determines the eigenvalues by a similar requirement as in the previous case, namely that jk = √αkb is a root of the Bessel function J−ν. Again the second constant is evaluated using Eq. (7.94), and the normalized functions are
uk(x) =
√2bν−1 xbν
J−ν jkx b
J1−ν(jk) . (7.98)
It is easy to prove that the eigenfunctions given by Eqs. (7.95) and (7.98) fulfill the normalization condition (7.91). Evaluating the scalar products we get after some tedious but straightforward algebra:
∞
X
k=1
huk,1iuk =
∞
X
k=1
212a1−ν jk
21−νjkν−1 Γ(ν)Jν+1(jk)+ 1
212aν−1
Jν+1(jk)zνJν(jkz)
=
∞
X
k=1
22−νzνjkν−1
jkΓ(ν)Jν+12 (jk)Jν(jkz) +
∞
X
k=1
2zν
jkJν+1(jk)Jν(jkz) (7.99)
= zν
∞
X
k=1
22−νjkν−1
jkΓ(ν)Jν+12 (jk)Jν(jkz) +zν
∞
X
k=1
2
jkJν+1(jk)Jν(jkz),
where we have set z =x/b. Making use of the Fourier-Bessel expansion zν =
∞
X
k=1
2
jkJν+1(jk)Jν(jkz), (7.100) we obtain
∞
X
k=1
huk,1iuk=zν
" ∞ X
k=1
2 jk
2−ν
1
Γ(ν)Jν+12 (jk)Jν(jkz) +zν
#
, (7.101)
and identifying the sum with the Fourier-Bessel expansion of z−ν −zν the desired property is shown.
Fig. 7.1 shows the comparison between the theoretical curves and simulations for the case of an entrance boundary at the origin.
0 5 10 15 20
Tb
0.05 0.10 0.15 0.20 0.25 fHTbL
n=1, x=1, b=4
0 2 4 6 8 10
Tb
0.1 0.2 0.3 0.4 fHTbL
n=2.5, x=1, b=4
0 1 2 3 4 5
Tb
0.2 0.4 0.6 0.8 1.0 1.2 1.4 fHTbL
n=1, x=1, b=2
Figure 7.1: First passage time probability density function for an entrance boundary at the origin. The parameters n, x and b are varied. The diffusion coefficient is alwaysD= 1. The continuous lines are the theoretical curves obtained by truncating the sum (7.90) after the first 200 terms, and the dots are the values obtained by Monte Carlo simulations.
In the case of an entrance boundary a further refinement of the simulation algo-rithm is possible. Knowing that the zero level can never be reached from the interior of the state space of the process, it is clear that negative values in the simulations must result from discretization errors. If this is the case we can reduce the time step until the propagated value of the process is positive.
7.6.2 Exit boundary
If the zero level is an exit boundary it is impossible to reach any interior point b, provided that the starting point of the process is “sufficiently close” to the boundary.
Hence the first passage time will not converge in general. However, the analysis of the previous sections together with an absorbing boundary condition at the origin delivers the probability density of the first exit time T0,b = min{T0, Tb} from the interval (0, b). For n≤ −1 the first terms in Eqs. (7.92) or (7.96) vanish in the limit x→0, whereas the second terms are finite. Therefore the constantB must be zero in order to fulfill the required boundary conditionu(b) = 0. Thus the eigenfunctions are given by Eq. (7.95), just as in the case of an entrance boundary for ν ∈Z.
In Fig. 7.2 the theoretical curves are again compared with the results obtained via Monte Carlo simulations. It is interesting to notice that the density is bimodal for some choices of the parameters.
0 1 2 3 4 5
T0,b
0.5 1.0 1.5 2.0
fHT0,bL
n=-1, x=1, b=4
0 2 4 6 8 10
T0,b
0.2 0.4 0.6 0.8
fHT0,bL
n=-1, x=4, b=5
0 2 4 6 8 10
T0,b
0.1 0.2 0.3 0.4 0.5 0.6 0.7
fHT0,bL
n=-2, x=4, b=5
0 2 4 6 8
T0,b
0.1 0.2 0.3 0.4 0.5 0.6 0.7
fHT0,bL
n=-2.5, x=4, b=5
Figure 7.2: First exit time probability density functions for an entrance boundary at the origin. The continuous lines represent again the theoretical curves and the dots the points obtained from the simulations.
7.6.3 Regular boundary
The case when the origin is a regular boundary is the most complicated one. In Ref. [99] the regular boundary is described as follows:
For a regular boundary a variety of boundary behaviour can be prescribed in a consistent way, including the contingencies of complete absorption or reflecting, elastic or sticky barrier phenomena, and even the possibility
of the particle (path), when attaining the boundary point, waiting there for an exponentially distributed duration followed by a jump into the interior of the state space according to a specified probability distribution function. In the latter event, the process only exhibits continuous sample paths over the interior of the state space.
However, we first observe that imposing an absorbing boundary condition at the origin gives us the probability densities of the first exit times. The eigenfunctions are computed in the same way as in the case of an exit boundary at the origin, and are again given by Eq. (7.95). Fig. 7.3 illustrates the results.
0 2 4 6 8
T0,b
0.2 0.4 0.6 0.8 1.0 1.2 1.4
fHT0,bL
n=-0.5, x=1, b=4
0 2 4 6 8 10
T0,b
0.2 0.4 0.6 0.8
fHT0,bL
n=-0.5, x=4, b=5
0 2 4 6 8 10 12 14
T0,b
0.1 0.2 0.3 0.4
fHT0,bL
n=0.5, x=1, b=4
0 2 4 6 8 10 12 14
T0,b
0.2 0.4 0.6 0.8 1.0
fHT0,bL
n=0.5, x=4, b=5
Figure 7.3: First exit time probability densities in the case when the origin is a regular boundary
If n > 0 the first passage time probability density is approximated well for small times imposing a reflecting boundary at the origin, but the maximum of the distribution obtained from the simulations is higher then the one in the theoretical curve, as Fig. 7.4 indicates.
0 5 10 15 20 25 30
Tb
0.05 0.10 0.15 0.20 fHTbL
n=0.5, x=1, b=4
0 5 10 15 20
Tb
0.05 0.10 0.15 0.20 0.25 fHTbL
n=0.8, x=1, b=4
Figure 7.4: First passage time probability densities for a regular boundary at the origin.
2 5 10 20 50 100 200
Tb 10-4
0.001 0.01 0.1 1 fHTbL
Figure 7.5: Logarithmic plot of first passage time densities for a regular boundary at the origin. The blue curve corresponds to the choice n = 0.5, x = 1, b = 4, and the red one to n= 0.8, x= 1, b= 4.
Fig. 7.5 shows a logarithmic plot of the first passage time densities, and one can see that the latter are heavy-tailed, i.e. they exhibit a power law decay for long times. A possible interpretation of the power law decay for a regular boundary at the origin as opposed to the exponential decay observed for the entrance boundary are the possible zero crossings of the path and the positive amount of time spent in the vicinity of the “sticky” boundary.
It is interesting to note that sticky boundaries may have applications e.g. for the simulation of the partial adsorption of polymer molecules to walls and for the modeling of solvent quality [102].