- Email: [email protected]
Green's function approach to the solution of the time dependent Fokker-Planck equation with an absorbing boundary
Physica 135A (1986) 63-79 North-Holland, Amsterdam
GREEN’S FUNCTION APPROACH TO THE SOLUTION OF THE TIME DEPENDENT FOKKER-PLANCK EQUATION WITH AN ABSORBING BOUNDARY S.V.G. MENON, Vinod KUMAR
and D.C. SAHNI
Theoretical Physics Division, Bhabha Atomic Research Centre, Bombay,
400 085, India
Received 27 May 1985 Revised 2 September 1985
The solution of the time dependent multivelocity Fokker-Planck equation in plane geometry with an absorbing boundary is formulated in terms of the infinite medium Green’s function and the boundary distribution. The boundary distribution is to be determined by solving an integral equation. It is shown that the velocity moments such as density, current and energy density can be expressed in terms of three reduced distributions which depend on the longitudinal velocity component only. Numerical results for monovelocity pulsed and steady sources are presented.
1. Introduction Recently there has been a lot of interest in the solution of the Fokker-Planck Equation (FPE) describing Brownian motion near an absorbing boundaryle4). In this connection, Harris’) has applied the method of half-range expansions to the one-dimensional (1-D) FPE and obtained an analytical expression for the density of Brownian particles. However, Razi Naqvi et al.‘) have shown that similar results can be deduced without resorting to half-range treatment. Burschka and Titulaer3) and Mayya and Sahni4) have used an eigenfunction expansion method (EEM) for solving the time independent 1-D FPE. The EEM essentially reduces the problem to solving a set of algebraic equations for the expansion coefficients. However, as pointed out by Burschka and Titulaer3) even a very large order expansion (-140) does not yield converged results for the integral parameters such as Milne extrapolation length, boundary density and energy density. Therefore, Mayya and Sahni4) attempted to sum up the infinite series using an asymptotic form of expansion coefficients after making some adhoc assumptions about the emergent distribution. They have obtained 0378-4371/86/$03.50 0 Elsevier Science Publishers B.V. (North-Holland Physics Publishing Division)
S.V.G.
64
an expression obtained
for density
expressions
ing the emergent exploited to solve
a relation
near
connecting
the 1-D FPE.
In their
et al
the boundary.
for expansion
distribution.
MENON
coefficients
In a recent the albedo work
Later
Titulaer’)
and Mayya’)
for large order without
work, kernel
Selinger
and the boundary
the albedo
prescrib-
and Titulaer’)
kernel,
which
have
distribution expresses
the
velocity distribution attained by a particle (of a given velocity) injected into the system (at the boundary) on its first return, is found via a numerical simulation of the Langevin equation equivalent to the FPE. The EEM has been applied by Waldenstrom et a1.8) for solving FPE in plane gemetry (but in three dimensional (3-D) velocity space) to calculate the Milne extrapolation length and particle density (L,P, approximation in their notation). also been treated using EEM’-‘I) f or calculating particles on a spherical droplet. It has been
FPE in spherical geometry has condensation rate of Brownian found that EEM fails to give
converged results”) for droplets of small radii. In this paper we develop the Green’s function method
(GFM)
for solving
the
time dependent FPE in plane geometry (but in 3-D velocity space) with an absorbing boundary. Our idea is to use an integral representation for the distribution function satisfying the boundary condition (BC) (see section 2) in terms of the infinite medium Green’s function (GF). This integral representation contains the distribution at the boundary which should be determined by solving an integral equation. We then introduce three reduced distribution functions (see section 3) which depend only on the longitudinal (along the space co-ordinate) velocity component. The velocity moments of general interest (density, current and completely determined
energy density) are then shown (see section in terms of the reduced distribution functions
forms one of the new results in this paper. density and Milne extrapolation length obtained
from the 1-D FPE. Therefore
4) to be and this
This reduction shows that the particle are exactly same as that would be
the treatment
of the full 3-D FPE (as has
been done in ref. 7) is needed only if some other parameters For the absorbing barrier problem of 1-D FPE, the integral the distribution function was proposed by Harris]‘) on
are to be studied. representation for physical grounds.
However, we have shown’“) that this integral representation can be derived via that adjoint GF of FPE. In ref. 13 we had also derived an expression for the Milne extrapolation length in terms of the boundary distribution (see section 5) which is explicitly determined in GFM. In section 6 we present, for the first time, a numerical scheme to solve the integral equation for the boundary distribution. In section 7 numerical results for the 1-D FPE are presented for two cases, viz. a pulsed source and a steady source. The steady source case for sufficiently large time closely approximates the solution of time independent FPE ‘). It is found that the results of our numerical procedure converge much faster than those obtained by EEM “). Finally section 8 gives our conclusions.
GREEN’S FUNCTION APPROACH TO FOKKER-PLANCK
65
EQUATION
2. Green’s function formalism for FPE We consider the general function W(r, u, t):
time dependent
FPE 14,1s) for the distribution
Vu *(VW + V”W) + Q(r, u, t) ,
aWlat + u*yW=
(1)
where the space (J), velocity (u) and time (t) variables are in appropriate units. Eq. (1) is to be solved with a specified source distribution Q(r, u, t). We assume that the medium extends over the semi-infinite half space, --cc < y, z < x), O
W(O,Y,
2;
v,,
v,, t) = KW(O,Y,
VY’
-v,,
2;
VY’
v,, t),
u, >o,
(2)
where K is a parameter. The limiting values, K = 0 and K = 1, correspond to perfectly absorbing and specularly reflecting boundaries respectively*). In order to cast eq. (1) and BC of eq. (2) in the integral form, we first consider the infinite medium (-m
+ u .V,G = V,, * (UC + V,G) + 6(r - r’)6(u - u’)8(t
- t’) .
(3)
The GF (solution of eq. (3)) has been given by Chandrasekhar14): G = (87rA))3’2 exp[- 1 {aVi - 2hV, * (R + V) + b(R + V)‘}] .
(4)
In eq. (4) and the following we use the abbreviations V=u-u’-(V,,l$,V,),
R=r-r’=(X,Y,Z),
V, = u - u’ exp(-7)
e, = 1 - exp(-r),
= (V,,, Vdy, V,,), e2 = 1 - exp(-2r)
a = 27/A,
A = 2re, - 4e:,
r=t-t’,
Following the procedure for eq. (1):
h = 2e,lA,
(5)
, b = e,lA .
outlined in ref. 13 we obtain an integral representation
W(r,u,t)=jdr’jdx’jdy’jdl./dufG(r,u,t;r’,u’,t’)Q(r’,u’,t’)
d; j
+ [ dtj 0
-cc
pm
d:‘, (m)
d:
G(r_v.
,I%, y’, t’, u’, t’)v;W(O,
y’, z’, u’, t;&
S.V.G. MENON
66
In eq.
(6),
the first term
distribution
Q(r,
particles
is the contribution
u, t) of particles.
(of velocity
u:) crossing
to W(r, u, t) from
The second
the source
term is the contribution
the boundary
crossing being uiW(O, y, z, v’, t’). We now specialise eq. (6) to the plane Q(x,
et al
due to
(X = 0) at time t’, the rate of
symmetric
problems,
i.e., Q(r,
u, t) =
u, r), W(r, u, t) = W(x, u, t) and the BC is
wo, u,, uy, u,, Therefore
t) = KW(O, -u,, uy, u*, t),
eq. (6) reduces
u,
>o
to
W(x, u, t) = / dt ‘/dl’j0 0
du’ G(x, u, t; x’, u’, c’)Q(x’,
u’, t’)
(-1 du’ G(x, u, t; 0, u’, t’)u;W(O,
0
(7)
(8)
u’, t') ,
(ml
where x
X
(9)
dz’ G(r, u, t; r’, u’, t') -s Integrations
-x in eq. (9) can be easily
C? = G,(2vAb)-’ where
G., is defined
exp[-(Vi,.
carried
+ V5,)/2bA],
(10)
as
G., = (271-A)-“* exp[ - : { uVf, - 2hV,,(X We observe
out to obtain
+ V,) + b(X + Vl)2}] .
that the BC of eq. (7) can be imposed
exactly
(11)
in the representation
of eq. (8) for W. Taking the limit x-+0 and using the BC of eq. (7) one can get an integral equation for the boundary distribution W(0, u, f) u, < 00). In the next section we shall consider certain (-~
3. Reduced
distribution
The velocity
moments
functions of distribution
function
which are of general
interest
GREEN’S
FUNCTION
APPROACH
TO FOKKER-PLANCK
EQUATION
67
are the density n(x, t), current J(x, t) and energy density E(x, t) which are defined as n(x, t) =
i (=)
E(x, t) =
I (=I
du W(x, u, t),
J(x, t) =
d u uM’(x, u, t), i Cm)
du u%V(x, u, t) .
(12)
In order to compute the density n(x, t) and current component JX(x, t), it is enough to consider the reduced distribution function f(x, u,, t) defined as Cc
m
Ax, u,, t) = I dUYI du, W(x, u, t) . -02 --m
(13)
For instance, the density is given by
6
I
4 = du, f(x, u,, t> , -m
and the current component
can be expressed as
(14) Integrating eq. (8) with reference to uYand u, we find that the reduced distribution function f(x, u,, t) is given by
f(x,Y,.‘)=jdr’jdx’jdu:G,(x,u,,r;x’,u:,Ir)a,(x~,u~,~‘) 0
-m
0 I
cr
du: G,(x, u,, t; 0, u:,q$f(o, 0
where
Q,
Dz
--m
u:,t’)
7
(15)
-m
m
-a
du, Qk
u, 4.
We note that starting from the 1-D FPE one can obtain eq. (15) in an analogous manner13). Thus, as mentioned in the introduction it is not necessary to solve
S.V.G. MENON
68
et al.
eq. (1) (as has been done in ref. 7) for determining the density current J,(x, t) and parameters related to them. For computing component F,(x,
1,.(x, t)
we
define
the
second
reduced
n(x, t) and the current
distribution
function
u,? t):
(16) so that
The current component JZ can be expressed and (16) it follows that F, is given by
in a similar
manner.
Using eqs. (8)
F,(x,u,,r)=~d~~~dr~~du~~~,(x,u,,r:~~,u~.~~)e~’Q~(x’.u~.I’) ,, -7 0 I
x
u.,, t; 0, II:, t’) e ’ u:F,(O,
dn; G,(x,
u:, t’),
(17)
x
0 where
Q;(x,
Normally,
u,,
t) =
1 1 du,
it is assumed
that
the
(18)
u, t) .
duz u,Q(x,
source
has
azimuthal symmetry so that u, = (ut + u:)“~. Then Qi
Q = Q(x, u,> u,, t) and W= W(x, II,, u,, t) where
and F, will vanish and Jy = Jz = 0. The energy density E(x, t) in eq. (12) can be rewritten x
x
E(x,
as
t> = I du, uf.,f
-x
where we have introduced defined as
u,, r) +
I -7.
the third
du, F&t, u,, t) ,
reduced
du, (u,‘; + u;)W(x,
FAX, q, t> = -r
-z
distribution
u, t) .
(19) function
K(x,
uY, t)
(20)
GREEN’S
FUNCTION
APPROACH
TO FOKKER-PLANCK
EQUATION
69
From eqs. (8) and (20) it can be seen that F2 satisfies
F2(x, uX, t) = i dt ~jdx’jdu:G~(x,u,,t;x’,u:,t’) -cc
0
0
x [2bAQ,(x’,
u:, t’) + e-*’ Qtl
dv; G,(x, u,, t; 0, u;, t’) 0
-m
x u32bAf(O, u:, t’) + em2’ F,(O, u:, t’)] ,
(21)
where
Qf
(22)
u, t) .
--P
Thus we see that f(x, v,, t), Fl(x, uX, t) and F2(x, uX, t) can be obtained once their boundary values are known. It may be remarked that determining these boundary values from eqs. (15), (17) and (21) (after taking the limit x-0) is simpler than solving for W(0, u, t) from eq. (8) with its three velocity components (though eq. (8) contains information regarding all other velocity moments as well). In the next section we derive the expressions for particle density, current and energy density.
4. Moments of distribution function Integrating eq. (15) over uX (--M < u, < m) we obtain the expression particle density n(x, t):
for
.(x,t)=jdt’jdxrjdu;G~(x,t;x’,u;,t’)Q,(xf,u;,t’) 0
0
-m
du: G,O(x, t; 0, u;, t’)u;f(O,
where
u;, t’) ,
(23)
70
S.V.G. MENON
Gy(x,
t; x’, u:, t’) = (477-y)-“*
exp[-(X
et al.
- u:e,)*/4y], (24)
y=~In a similar expression
$(3-4e-‘+em2’)
.
manner
eq. (15) by u, and integrating
multiplying
we obtain
the
for Jx(x, t):
x VX, 4, T)(A&)Q,(x’,
,
P
dt’
+
I
u:,f>
i
du: G:(x,
t; 0, u:, t’)U(x,
u:, 7)(A/2y)u;f(O,
u:, t’)
-cc
0
,
(25)
where
uw, u:,T) = au: epT + h[X Also, other current (17) as
component
- u:( 1 + e-‘)I J, (and similarly
- b(X - II:)
(26)
J,) can be obtained
from eq.
J,(x,~)=jdr’jdx’jdu~~~(x,l;x’;u~.r’)e-’Q:(x’.u:,f’) 0
0
-x
I +
i
r dt’
i
du: Gz(x, t; 0, u:, t’) emT u:F,(O,
(27)
u:, t’) .
Paz
0
Using eqs. (15), (19) and (21) we can obtain density E(x, t):
the following
expression
for energy
E(x,t)=jdl’jdx.idu:G:(x,r;x’,u:,i’) 0 0 -= x [{2e, + H(X, I +
I 0
u:, 7)}Q,(x’,
u:, t’) + ezT
Q~(JL’, u:, t’)]
m dt’
I -=
+ e-“&(0,
du; G;(x,
II:, t’)] ,
t; 0, u;, f)[{2e,
+ wx,
u:>
T)>f(0,
u:, t’) (28)
GREEN’S FUNCTION
APPROACH
TO FOKKER-PLANCK
EQUATION
71
where H(X, u;, T) = (A/2y)[l
+ (AI2y)U(X,
u;, T)] .
(29)
In the case of 1-D FPE, energy density E,(x, t) is given by the first term on the rhs of eq. (19) and therefore we have
E&,t)=/dt’[dx’[
du: G:(x, t; x’, u:, t’)H(X,
0 I +
u;, T)Q,(x’,
u;, t’)
-m
0 m
dt’
i
du; G,O(x, t; 0, u:, t’)H(x,
i
u;, t’) .
u:, +:f(O,
(30)
-m
0
For obtaining the moments it is necessary to determine the boundary distributionsf(0, u,, t), F,(O, u,, t) and F,(O, u,, t). Integral equations for these boundary distributions can be obtained from eqs. (15), (17) and (21) by taking the limit x+0. For instance f(0, u,, t) satisfies (--co < u, 5 0) f(O,u~,t)=jdtfjdx’jdu~G,(O,u,.t;~r,u~,t~)Q~(~~,u~,t’) 0
-m
0 0
I +
i
dt’
_I
du: [G,(O, u,, t; 0, u;, t’) - KGJO,
ux, t; 0, -II;,
t’)]
0 x
uLf(O,
u:,
t’)
(31)
.
In obtaining eq. (31) we have used the BC
which follows from definition of f eq. (13) and eq. (7). It is also possible to obtain an integral equation of first kind for f(0, u,, t) (0 < u, < 00) and is given by O=idf’ -m
0
0 , +
I 0
x
0
dt’
I
du: [G,(O, u,, t; 0, u;, t’) - KG,(O,
u,, t; 0, -u;,
t’)]
-m
u;f(0,
u;, t’) .
(33)
72
S.V.G.
Either
particle
In the
density
5. Diffusion
next
section
A procedure discussed
to extract
the diffusion
in ref. 13. There
I
+
nels.
I 0
scheme,
distribution
to be discussed
diffusion
approximation
to
conditions.
length
approximation
it has been
shown
n,(x,
t) to n(x, t) has
that for large time and far
eq. (23) for n(x, t) can be approximated
as
7z
dt’
du: [k(x, T) + h(x, T)U; /2]u;f(O,
I -2
u.1, t’) ,
exp(-X2/47) and h(X, 7) = Xk/T are diffusion shown that n, satisfies the diffusion equation
k(X, r) = (47~)“~ It has also been
dn,ldt
in
u,, t) and then eq. (33) to verify the
and Milne extrapolation
away from the boundary,
where
we discuss
n(x, t) and its boundary
approximation
the boundary
in our numerical
6, we have used eq. (31) to computef(0,
its accuracy.
been
et al
eq. (31) or eq. (33) can be used to determine
f(0, u, , t) (-CC < u, 5 0). However section
MENON
= d’n,ldx2
+ QP(x, t) ,
(34)
ker-
(35)
where
Eq. (35) is to be solved using the mixed boundary radiation boundary condition)
anDlax+ A(t)n,
= 0
at x = 0,
condition”)
(also known
as
(36)
where
(37)
which
may be rewritten
as
GREEN’S FUNCTION APPROACH TO FOKKER-PLANCK
EQUATION
73
A(t) = [Cl- K)l(l+ WI In the terminology of neutron transport theory A-‘(t) is known as the Milne extrapolation length. From eq. (37) we note that A-‘(t) is equal to the ratio of average values of us and uXat the boundary. Therefore it can be obtained purely from the knowledge of boundary distribution.
6. Numerical approximations In this section we discuss the numerical scheme to determine the boundary distributions. We shall illustrate this procedure with reference to eq. (31) though F,(O, u,, t) and F,(O, uX, t) can be determined in a similar fashion. All the numerical calculations have been performed for the case of perfectly absorbing boundary (K = 0). For solving eq. (31) Q, should be known and we have chosen two cases, (1) Q,(x, u,, t) = 8(x - x,)S(u, - r@(t) and (2) Q,(x, uX, t) = S(x - x,)S(u, - u,). The former choice corresponds to a monovelocity (u,) pulsed source (t = 0) located at x = x,, while the latter corresponds to a steady source emitting particles for t Z 0 with velocity us at point x,. The solution of eq. (31) has been effected by discretising time and velocity variables. We use a constant time step S and assume that f(0, u,, t) varies linearly in a time step, i.e.
m
ux>0 = M-l(%)(L - 4 +f,(u,)(t- ~,-dl&
t,_,
wheref,(u, >= IV, u,, n8 ). Further, the integration over velocity variable (rhs of eq. (31)) is replaced by a high order Gauss-quadrature formula (order J) with weight function exp(-ut/2) over (--03 < uX5 0). With these approximations, eq. (31) can be rewritten as (15 i 5 J)
fi(“i)
=
‘nC”i?‘s,
us) + 6 i: 5 m=l j=l
x *dt’((l - t’)f,(uj)
ujwj exp(uT/2)
+t’f,_,(uj)}G,(O,
vi, (n - m + t’)S; 0, uj, 0)
,
I 0
(39)
where ui and wi are the quadrature directions and weights respectively and S, is the contribution from the source term Q, . For the pulsed source case S, is given by
S.V.G. MENON
74
et al.
(40) while for the steady
source
case 1
S,(U,> x,, us) =
S i j- dt' G,(O, u,, ???=I 0
(n - m + t’)S; x,,
us,
0)
(41)
The integration over the variable t’ in eqs. (39) and (41) is carried out by Simpson’s rule to a specified accuracy. Rewriting eq. (39) in matrix notation we get n-1 (I-
C,,lf,
=
s,
+
??I=1 II m m c
LB
f
-1
+
cn-mLJ+ 4.x-L 7
wherefn and S, are Z-dimensional vectors the .Z x .I matrices B, and C,,, are given
From
and Z is a unit matrix. by
lBrnlij=
1dt'
t’G,(O, u,, (m + t’)6; 0, uj, 0) ,
[C,Jij =
1dt'
(1 - t’)G,(O, u,, (m + t’)6; 0, u,, 0).
eq. (42) it is clear thatf,,
can be solved
matrix (I - C,) is to be inverted only once Milne extrapolation length can be calculated
n-‘(ns) =
5 l+v; &/i:
I=1
wT = W, exp(ufi2).
u!wT f,C”i>
,=I
recursively.
(42) The elements
It may be noted
since it is independent as
of
that
of II. The
(43)
T
f(0, u,, t) is known
at time
steps t, = nS, one can compute the function f(x, u,, t) or its moments JX(x, t) and E,(x, t) as given by eqs. (23), (25) and (30).
n(x, t),
where
Once
the distribution
7. Results In figs. l-3 we give the boundary distribution f(0, u,, t), the particle n(x, t) and the mean square velocity (ut),, = E,(x, t) ln(x, t) for a steady (switched at t = 0) at time t = 500. The same parameters calculated pulsed source at large t (-100) were almost indistinguishable and hence
density source for the are not
GREEN’S
FUNCTION
APPROACH
TO FOKKER-PLANCK
EQUATION
75
0.7r
vxFig. 1. Normalised boundary v, = -1.0) at t = 500.
f(0, v,, t) vs v, for steady
distribution
001 0.0
1
1
1.0
2.0
3.0
source
(with x, =20.0
and
I
LD
X--,
Fig. 2. Normalised particle q= -1.0) at t=500.
density
n(x, t) ln(0, t) vs x for the steady
source
(with x, = 20.0 and
repeatedly given. The boundary distribution and density plots are normalised such that n(0, t) = 1.0. In the numerical calculations we have used x, = 20.0, u, = -1.0 and a 15th order Gauss-quadrature formula16) (for the velocity variable integration) for both cases. We have also used the time steps S = 0.1 and 6 = 0.5 for the pulsed source and steady source respectively. It was found necessary to use a high accuracy (-10w4) Simpson rule for time integration in computing the source terms in eq. (41) and the matrices B, and C,. Further reduction of time steps was then found to have negligible effect on the boundary distribution and other parameters. This is expected to be due to the linear variation of f(0, u,, t) assumed in a time step. It can also be shown that f(0, uX, t)
76
S.V.G.
MENON
et al.
1.6
Fig. 3. Mean square
velocity
(ut),,
vs x for steady
source
varies as tm3” for the pulsed
source
source for sufficiently large We note that the boundary lian. For positive values of zero since K = 0. However, eq. (33) for positive u,. The
values of t. distribution in u, the boundary we have plotted slight departure
mainly
due
to the
finite
order
and becomes
(with x, = 20.0 and u, = ~ 1.0) at
independent
I = 500.
oft for the steady
fig. 1 resembles a shifted Maxweldistribution should be identically the numerical values of the rhs of from zero for U, > 0 is found to be
of quadrature
(15th)
employed
for velocity
discretisation. We observe that these deviations are much smaller than those obtained by Burschka and Titulae?) using EEM of 140th order (compare fig. 1 of this paper with fig. 2 of ref. 3). We note that Titulae?) and Mayyah) had shown that the boundary distribution approaches zero as u, + 0. However, the EEM “) and a numerical simulation’) lead to non-zero values (-0.186) while the present calculation boundary distribution and the present
method
yields a value (-0.026) much closer to zero. In fact the for small u, as computed from numerical simulation’) are somewhat
different
corresponding formula given
to numerical simulation have in ref. 7 with the normalisation
unity. These accurate.
findings
suggest
that
GFM
as seen in table
I. The values
been obtained using the fitted that the area under the curve is
results
presented
here
are
more
In table II we give the numerical values of n(0, t) lJ,(O, t), (vf),, at x = 0 and Milne extrapolation length K’(t) for the pulsed source (at t = 100) and steady source (at t = 500). The steady source case for sufficiently large time should closely approximate the time independent FPE and therefore these results can be compared with those available in the literature”.“‘). Values of these parameters obtained by EEM of 140th order and extrapolated infinite order”,‘) and the numerical simulation’) are also given in table II for comparison. For
GREEN’S
FUNCTION
APPROACH
Boundary Numerical (ref. 11)
Velocity (u,) 0.0
TABLE I distribution
f(0, u,,
simulation
0.01 0.02 0.05 0.1 0.2 0.3 0.5 0.7 1.0 1.4 2.0 3.0 4.0
length
and boundary
EQUATION
t).
0.0260 0.0330 0.0886 0.1309 0.2087 0.2935 0.4093 0.4900 0.5875 0.6215 0.5843 0.4378 0.1960 0.0210 0.0008
TABLE II values of density
and mean
square
Type
n-‘(t)
n(0, t) /J,(O, I)
(u?),,
pulsed source at t = 100
1.459
0.937
1.556
steady source at t = 500
1.461
0.936
1.561
1.438
1.054
1.365
time independent FPE (extrapolated infinite order EEM “)
1.461
0.914
1.590
time independent FPE (extrapolated infinite order EEM-revised’)
1.460
0.951
1.535
simulation
1.454 * 0.004
0.942 ” 0.008
1.544
time independent
FPE
77
Steady source at r = 500
0.1857 0.1869 0.1976 0.2094 0.2439 0.2984 0.3958 0.4764 0.5852 0.6274 0.5904 0.4363 0.1920 0.0206 0.0007
0.001
Milne extrapolation
TO FOKKER-PLANCK
velocity.
at x = 0
(140th order EEM ‘))
results’)
instance, our value of 1.461 for A-’ (steady source) is to be compared with 1.438 and 1.461 obtained in the EEM of 140th order and extrapolated infinite order, respectively. Table III indicates the rate of convergence of our results with t for the steady source case. It should be remarked that there are significant differences between the integral parameters of the pulsed source and steady source cases (table II) even at large times. This can be understood from fig. 4 where we have plotted the current JX(x, t)lJx(O, t) for these two cases. For the steady source case for
78
S.V.G.
MENON
TABLE
Mime
Time
et al
III
extrapolation length and boundary values of density and mean square velocity as a function of f for the steady source. (t)
5 10 15 20 30 50 100 200 300 400 500
‘K’(t)
n(0, t)iJX(O, t)
(vf).,
3.0745 1.8903 1.6603 1.5762 1.5139 1.4507 1.4660 I .4620 1.4611 1.4608 1.4607
0.3606 0.6870 0.8093 0.8605 0.9006 0.9226 0.9324 0.935 1 0.9357 0.9359 0.9360
8.5261 2.7514 2.0517 1.8317 1.6810 1.6050 1.5722 1.5634 1.5616 1.5609 1.5606
at x = 0
1.021.01-
-2 3
0.96 -
Fig. 4. Normalised current profiles J,(x, r)lJ,(O, source at t = 500, (3) 1,(x, f)lJ,(O, t) = 1.
r) vs X, (1) pulsed
source
at
f = 100, (2) steady
sufficiently large times the current is almost flat spatially. Thus this case approximates the time independent FPE. However, as the current profile is quite different for the pulsed case and therefore its results cannot be compared directly with those of the time independent calculations.
8. Conclusions
In this paper we have formulated the GFM for solving the time dependent FPE with an absorbing boundary. The method consists of obtaining a representation for the distribution function (satisfying the BC) in terms of the boundary
GREEN’S
FUNCTION
APPROACH
TO FOKKER-PLANCK
EQUATION
79
distribution which should be determined by solving an integral equation. The 3-D velocity dependence of the distribution function has been incorporated in the formalism. It has been shown that for calculating the velocity moments of general interest (density, current and energy density) it is not necessary to solve for the full distribution function; it is enough to determine three reduced distribution functions (having one velocity component only). We have presented a numerical scheme to solve the integral equations for the boundary distributions. Numerical results have shown that this procedure yields results which are far superior to that obtained by the conventional EEM and in fact are closer to the revised value?). This is due to the fact that EEM results converge very slowly as the order of approximation is increased. Since the infinite medium GF for the time independent FPE is not available in closed form, we have treated this problem in the time dependent formalism using a steady source of particles away from the boundary. One of the drawbacks of the present method is that the infinite medium GF for the problem under consideration should be known. For Brownian motion under general external force fields this is not known; however, for the case of quadratic potentials (in space variables) the Green’s functions are available in literature14’15 ). The present formalism can also be extended to solve the FPE in spherical geometry for treating the condensation problem of a spherical droplet”).
References 1) 2) 3) 4) 5) 6) 7) 8) 9) 10) 11) 12) 13) 14) 15) 16)
S. Harris, J. Chem. Phys. 75 (1981) 3103. K. Razi Naqvi, K.J. Mork and S. Waldenstrom, Phys. Rev. Lett. 49 (1982) 304. M.A. Burschka and U.M. Titulaer, J. Stat. Phys. 25 (1981) 569. Y.S. Mayya and D.C. Sahni, J. Chem. Phys. 79 (1983) 2302. U.M. Titulaer, J. Stat. Phys. 37 (1984) 589. Y.S. Mayya, J. Chem. Phys. 82 (1985) 2033. J.V. Selinger and U.M. Titulaer, J. Stat. Phys. 36 (1984) 293. S. Waldenstrom, K.J. Mork and K. Razi Naqvi, Phys. Rev. A 28 (1983) 1659. S. Waldenstrom, K. Razi Naqvi and K.J. Mork, Ark. Fys. Semin. Trondheim, No. 8, ISSNO 365-2459 (1983). D.C. Sahni, Phys. Rev. A 30 (1984) 2056. Vinod Kumar and S.V.G. Menon, J. Chem. Phys. 82 (1985) 917. S. Harris, J. Chem. Phys. 72 (1980) 2659. A 32 (1985) 3832. S.V.G. Menon and D.C. Sahni, Phys. Rev. S. Chandrasekhar, Rev. Mod. Phys. 15 (1943) 1. H. Risken, The Fokker-Planck equation; Method of Solution and Applications, Springer series in Synergetics, H. Hanken, ed. (Springer, Berlin, 1983). N.M. Steen, G.D. Byrne and E.M. Gelbard, Math. Comp. 23 (1969) 665.
GREEN’S FUNCTION APPROACH TO THE SOLUTION OF THE TIME DEPENDENT FOKKER-PLANCK EQUATION WITH AN ABSORBING BOUNDARY S.V.G. MENON, Vinod KUMAR
and D.C. SAHNI
Theoretical Physics Division, Bhabha Atomic Research Centre, Bombay,
400 085, India
Received 27 May 1985 Revised 2 September 1985
The solution of the time dependent multivelocity Fokker-Planck equation in plane geometry with an absorbing boundary is formulated in terms of the infinite medium Green’s function and the boundary distribution. The boundary distribution is to be determined by solving an integral equation. It is shown that the velocity moments such as density, current and energy density can be expressed in terms of three reduced distributions which depend on the longitudinal velocity component only. Numerical results for monovelocity pulsed and steady sources are presented.
1. Introduction Recently there has been a lot of interest in the solution of the Fokker-Planck Equation (FPE) describing Brownian motion near an absorbing boundaryle4). In this connection, Harris’) has applied the method of half-range expansions to the one-dimensional (1-D) FPE and obtained an analytical expression for the density of Brownian particles. However, Razi Naqvi et al.‘) have shown that similar results can be deduced without resorting to half-range treatment. Burschka and Titulaer3) and Mayya and Sahni4) have used an eigenfunction expansion method (EEM) for solving the time independent 1-D FPE. The EEM essentially reduces the problem to solving a set of algebraic equations for the expansion coefficients. However, as pointed out by Burschka and Titulaer3) even a very large order expansion (-140) does not yield converged results for the integral parameters such as Milne extrapolation length, boundary density and energy density. Therefore, Mayya and Sahni4) attempted to sum up the infinite series using an asymptotic form of expansion coefficients after making some adhoc assumptions about the emergent distribution. They have obtained 0378-4371/86/$03.50 0 Elsevier Science Publishers B.V. (North-Holland Physics Publishing Division)
S.V.G.
64
an expression obtained
for density
expressions
ing the emergent exploited to solve
a relation
near
connecting
the 1-D FPE.
In their
et al
the boundary.
for expansion
distribution.
MENON
coefficients
In a recent the albedo work
Later
Titulaer’)
and Mayya’)
for large order without
work, kernel
Selinger
and the boundary
the albedo
prescrib-
and Titulaer’)
kernel,
which
have
distribution expresses
the
velocity distribution attained by a particle (of a given velocity) injected into the system (at the boundary) on its first return, is found via a numerical simulation of the Langevin equation equivalent to the FPE. The EEM has been applied by Waldenstrom et a1.8) for solving FPE in plane gemetry (but in three dimensional (3-D) velocity space) to calculate the Milne extrapolation length and particle density (L,P, approximation in their notation). also been treated using EEM’-‘I) f or calculating particles on a spherical droplet. It has been
FPE in spherical geometry has condensation rate of Brownian found that EEM fails to give
converged results”) for droplets of small radii. In this paper we develop the Green’s function method
(GFM)
for solving
the
time dependent FPE in plane geometry (but in 3-D velocity space) with an absorbing boundary. Our idea is to use an integral representation for the distribution function satisfying the boundary condition (BC) (see section 2) in terms of the infinite medium Green’s function (GF). This integral representation contains the distribution at the boundary which should be determined by solving an integral equation. We then introduce three reduced distribution functions (see section 3) which depend only on the longitudinal (along the space co-ordinate) velocity component. The velocity moments of general interest (density, current and completely determined
energy density) are then shown (see section in terms of the reduced distribution functions
forms one of the new results in this paper. density and Milne extrapolation length obtained
from the 1-D FPE. Therefore
4) to be and this
This reduction shows that the particle are exactly same as that would be
the treatment
of the full 3-D FPE (as has
been done in ref. 7) is needed only if some other parameters For the absorbing barrier problem of 1-D FPE, the integral the distribution function was proposed by Harris]‘) on
are to be studied. representation for physical grounds.
However, we have shown’“) that this integral representation can be derived via that adjoint GF of FPE. In ref. 13 we had also derived an expression for the Milne extrapolation length in terms of the boundary distribution (see section 5) which is explicitly determined in GFM. In section 6 we present, for the first time, a numerical scheme to solve the integral equation for the boundary distribution. In section 7 numerical results for the 1-D FPE are presented for two cases, viz. a pulsed source and a steady source. The steady source case for sufficiently large time closely approximates the solution of time independent FPE ‘). It is found that the results of our numerical procedure converge much faster than those obtained by EEM “). Finally section 8 gives our conclusions.
GREEN’S FUNCTION APPROACH TO FOKKER-PLANCK
65
EQUATION
2. Green’s function formalism for FPE We consider the general function W(r, u, t):
time dependent
FPE 14,1s) for the distribution
Vu *(VW + V”W) + Q(r, u, t) ,
aWlat + u*yW=
(1)
where the space (J), velocity (u) and time (t) variables are in appropriate units. Eq. (1) is to be solved with a specified source distribution Q(r, u, t). We assume that the medium extends over the semi-infinite half space, --cc < y, z < x), O
W(O,Y,
2;
v,,
v,, t) = KW(O,Y,
VY’
-v,,
2;
VY’
v,, t),
u, >o,
(2)
where K is a parameter. The limiting values, K = 0 and K = 1, correspond to perfectly absorbing and specularly reflecting boundaries respectively*). In order to cast eq. (1) and BC of eq. (2) in the integral form, we first consider the infinite medium (-m
+ u .V,G = V,, * (UC + V,G) + 6(r - r’)6(u - u’)8(t
- t’) .
(3)
The GF (solution of eq. (3)) has been given by Chandrasekhar14): G = (87rA))3’2 exp[- 1 {aVi - 2hV, * (R + V) + b(R + V)‘}] .
(4)
In eq. (4) and the following we use the abbreviations V=u-u’-(V,,l$,V,),
R=r-r’=(X,Y,Z),
V, = u - u’ exp(-7)
e, = 1 - exp(-r),
= (V,,, Vdy, V,,), e2 = 1 - exp(-2r)
a = 27/A,
A = 2re, - 4e:,
r=t-t’,
Following the procedure for eq. (1):
h = 2e,lA,
(5)
, b = e,lA .
outlined in ref. 13 we obtain an integral representation
W(r,u,t)=jdr’jdx’jdy’jdl./dufG(r,u,t;r’,u’,t’)Q(r’,u’,t’)
d; j
+ [ dtj 0
-cc
pm
d:‘, (m)
d:
G(r_v.
,I%, y’, t’, u’, t’)v;W(O,
y’, z’, u’, t;&
S.V.G. MENON
66
In eq.
(6),
the first term
distribution
Q(r,
particles
is the contribution
u, t) of particles.
(of velocity
u:) crossing
to W(r, u, t) from
The second
the source
term is the contribution
the boundary
crossing being uiW(O, y, z, v’, t’). We now specialise eq. (6) to the plane Q(x,
et al
due to
(X = 0) at time t’, the rate of
symmetric
problems,
i.e., Q(r,
u, t) =
u, r), W(r, u, t) = W(x, u, t) and the BC is
wo, u,, uy, u,, Therefore
t) = KW(O, -u,, uy, u*, t),
eq. (6) reduces
u,
>o
to
W(x, u, t) = / dt ‘/dl’j0 0
du’ G(x, u, t; x’, u’, c’)Q(x’,
u’, t’)
(-1 du’ G(x, u, t; 0, u’, t’)u;W(O,
0
(7)
(8)
u’, t') ,
(ml
where x
X
(9)
dz’ G(r, u, t; r’, u’, t') -s Integrations
-x in eq. (9) can be easily
C? = G,(2vAb)-’ where
G., is defined
exp[-(Vi,.
carried
+ V5,)/2bA],
(10)
as
G., = (271-A)-“* exp[ - : { uVf, - 2hV,,(X We observe
out to obtain
+ V,) + b(X + Vl)2}] .
that the BC of eq. (7) can be imposed
exactly
(11)
in the representation
of eq. (8) for W. Taking the limit x-+0 and using the BC of eq. (7) one can get an integral equation for the boundary distribution W(0, u, f) u, < 00). In the next section we shall consider certain (-~
3. Reduced
distribution
The velocity
moments
functions of distribution
function
which are of general
interest
GREEN’S
FUNCTION
APPROACH
TO FOKKER-PLANCK
EQUATION
67
are the density n(x, t), current J(x, t) and energy density E(x, t) which are defined as n(x, t) =
i (=)
E(x, t) =
I (=I
du W(x, u, t),
J(x, t) =
d u uM’(x, u, t), i Cm)
du u%V(x, u, t) .
(12)
In order to compute the density n(x, t) and current component JX(x, t), it is enough to consider the reduced distribution function f(x, u,, t) defined as Cc
m
Ax, u,, t) = I dUYI du, W(x, u, t) . -02 --m
(13)
For instance, the density is given by
6
I
4 = du, f(x, u,, t> , -m
and the current component
can be expressed as
(14) Integrating eq. (8) with reference to uYand u, we find that the reduced distribution function f(x, u,, t) is given by
f(x,Y,.‘)=jdr’jdx’jdu:G,(x,u,,r;x’,u:,Ir)a,(x~,u~,~‘) 0
-m
0 I
cr
du: G,(x, u,, t; 0, u:,q$f(o, 0
where
Q,
Dz
--m
u:,t’)
7
(15)
-m
m
-a
du, Qk
u, 4.
We note that starting from the 1-D FPE one can obtain eq. (15) in an analogous manner13). Thus, as mentioned in the introduction it is not necessary to solve
S.V.G. MENON
68
et al.
eq. (1) (as has been done in ref. 7) for determining the density current J,(x, t) and parameters related to them. For computing component F,(x,
1,.(x, t)
we
define
the
second
reduced
n(x, t) and the current
distribution
function
u,? t):
(16) so that
The current component JZ can be expressed and (16) it follows that F, is given by
in a similar
manner.
Using eqs. (8)
F,(x,u,,r)=~d~~~dr~~du~~~,(x,u,,r:~~,u~.~~)e~’Q~(x’.u~.I’) ,, -7 0 I
x
u.,, t; 0, II:, t’) e ’ u:F,(O,
dn; G,(x,
u:, t’),
(17)
x
0 where
Q;(x,
Normally,
u,,
t) =
1 1 du,
it is assumed
that
the
(18)
u, t) .
duz u,Q(x,
source
has
azimuthal symmetry so that u, = (ut + u:)“~. Then Qi
Q = Q(x, u,> u,, t) and W= W(x, II,, u,, t) where
and F, will vanish and Jy = Jz = 0. The energy density E(x, t) in eq. (12) can be rewritten x
x
E(x,
as
t> = I du, uf.,f
-x
where we have introduced defined as
u,, r) +
I -7.
the third
du, F&t, u,, t) ,
reduced
du, (u,‘; + u;)W(x,
FAX, q, t> = -r
-z
distribution
u, t) .
(19) function
K(x,
uY, t)
(20)
GREEN’S
FUNCTION
APPROACH
TO FOKKER-PLANCK
EQUATION
69
From eqs. (8) and (20) it can be seen that F2 satisfies
F2(x, uX, t) = i dt ~jdx’jdu:G~(x,u,,t;x’,u:,t’) -cc
0
0
x [2bAQ,(x’,
u:, t’) + e-*’ Qt
dv; G,(x, u,, t; 0, u;, t’) 0
-m
x u32bAf(O, u:, t’) + em2’ F,(O, u:, t’)] ,
(21)
where
Qf
(22)
u, t) .
--P
Thus we see that f(x, v,, t), Fl(x, uX, t) and F2(x, uX, t) can be obtained once their boundary values are known. It may be remarked that determining these boundary values from eqs. (15), (17) and (21) (after taking the limit x-0) is simpler than solving for W(0, u, t) from eq. (8) with its three velocity components (though eq. (8) contains information regarding all other velocity moments as well). In the next section we derive the expressions for particle density, current and energy density.
4. Moments of distribution function Integrating eq. (15) over uX (--M < u, < m) we obtain the expression particle density n(x, t):
for
.(x,t)=jdt’jdxrjdu;G~(x,t;x’,u;,t’)Q,(xf,u;,t’) 0
0
-m
du: G,O(x, t; 0, u;, t’)u;f(O,
where
u;, t’) ,
(23)
70
S.V.G. MENON
Gy(x,
t; x’, u:, t’) = (477-y)-“*
exp[-(X
et al.
- u:e,)*/4y], (24)
y=~In a similar expression
$(3-4e-‘+em2’)
.
manner
eq. (15) by u, and integrating
multiplying
we obtain
the
for Jx(x, t):
x VX, 4, T)(A&)Q,(x’,
,
P
dt’
+
I
u:,f>
i
du: G:(x,
t; 0, u:, t’)U(x,
u:, 7)(A/2y)u;f(O,
u:, t’)
-cc
0
,
(25)
where
uw, u:,T) = au: epT + h[X Also, other current (17) as
component
- u:( 1 + e-‘)I J, (and similarly
- b(X - II:)
(26)
J,) can be obtained
from eq.
J,(x,~)=jdr’jdx’jdu~~~(x,l;x’;u~.r’)e-’Q:(x’.u:,f’) 0
0
-x
I +
i
r dt’
i
du: Gz(x, t; 0, u:, t’) emT u:F,(O,
(27)
u:, t’) .
Paz
0
Using eqs. (15), (19) and (21) we can obtain density E(x, t):
the following
expression
for energy
E(x,t)=jdl’jdx.idu:G:(x,r;x’,u:,i’) 0 0 -= x [{2e, + H(X, I +
I 0
u:, 7)}Q,(x’,
u:, t’) + ezT
Q~(JL’, u:, t’)]
m dt’
I -=
+ e-“&(0,
du; G;(x,
II:, t’)] ,
t; 0, u;, f)[{2e,
+ wx,
u:>
T)>f(0,
u:, t’) (28)
GREEN’S FUNCTION
APPROACH
TO FOKKER-PLANCK
EQUATION
71
where H(X, u;, T) = (A/2y)[l
+ (AI2y)U(X,
u;, T)] .
(29)
In the case of 1-D FPE, energy density E,(x, t) is given by the first term on the rhs of eq. (19) and therefore we have
E&,t)=/dt’[dx’[
du: G:(x, t; x’, u:, t’)H(X,
0 I +
u;, T)Q,(x’,
u;, t’)
-m
0 m
dt’
i
du; G,O(x, t; 0, u:, t’)H(x,
i
u;, t’) .
u:, +:f(O,
(30)
-m
0
For obtaining the moments it is necessary to determine the boundary distributionsf(0, u,, t), F,(O, u,, t) and F,(O, u,, t). Integral equations for these boundary distributions can be obtained from eqs. (15), (17) and (21) by taking the limit x+0. For instance f(0, u,, t) satisfies (--co < u, 5 0) f(O,u~,t)=jdtfjdx’jdu~G,(O,u,.t;~r,u~,t~)Q~(~~,u~,t’) 0
-m
0 0
I +
i
dt’
_I
du: [G,(O, u,, t; 0, u;, t’) - KGJO,
ux, t; 0, -II;,
t’)]
0 x
uLf(O,
u:,
t’)
(31)
.
In obtaining eq. (31) we have used the BC
which follows from definition of f eq. (13) and eq. (7). It is also possible to obtain an integral equation of first kind for f(0, u,, t) (0 < u, < 00) and is given by O=idf’ -m
0
0 , +
I 0
x
0
dt’
I
du: [G,(O, u,, t; 0, u;, t’) - KG,(O,
u,, t; 0, -u;,
t’)]
-m
u;f(0,
u;, t’) .
(33)
72
S.V.G.
Either
particle
In the
density
5. Diffusion
next
section
A procedure discussed
to extract
the diffusion
in ref. 13. There
I
+
nels.
I 0
scheme,
distribution
to be discussed
diffusion
approximation
to
conditions.
length
approximation
it has been
shown
n,(x,
t) to n(x, t) has
that for large time and far
eq. (23) for n(x, t) can be approximated
as
7z
dt’
du: [k(x, T) + h(x, T)U; /2]u;f(O,
I -2
u.1, t’) ,
exp(-X2/47) and h(X, 7) = Xk/T are diffusion shown that n, satisfies the diffusion equation
k(X, r) = (47~)“~ It has also been
dn,ldt
in
u,, t) and then eq. (33) to verify the
and Milne extrapolation
away from the boundary,
where
we discuss
n(x, t) and its boundary
approximation
the boundary
in our numerical
6, we have used eq. (31) to computef(0,
its accuracy.
been
et al
eq. (31) or eq. (33) can be used to determine
f(0, u, , t) (-CC < u, 5 0). However section
MENON
= d’n,ldx2
+ QP(x, t) ,
(34)
ker-
(35)
where
Eq. (35) is to be solved using the mixed boundary radiation boundary condition)
anDlax+ A(t)n,
= 0
at x = 0,
condition”)
(also known
as
(36)
where
(37)
which
may be rewritten
as
GREEN’S FUNCTION APPROACH TO FOKKER-PLANCK
EQUATION
73
A(t) = [Cl- K)l(l+ WI In the terminology of neutron transport theory A-‘(t) is known as the Milne extrapolation length. From eq. (37) we note that A-‘(t) is equal to the ratio of average values of us and uXat the boundary. Therefore it can be obtained purely from the knowledge of boundary distribution.
6. Numerical approximations In this section we discuss the numerical scheme to determine the boundary distributions. We shall illustrate this procedure with reference to eq. (31) though F,(O, u,, t) and F,(O, uX, t) can be determined in a similar fashion. All the numerical calculations have been performed for the case of perfectly absorbing boundary (K = 0). For solving eq. (31) Q, should be known and we have chosen two cases, (1) Q,(x, u,, t) = 8(x - x,)S(u, - r@(t) and (2) Q,(x, uX, t) = S(x - x,)S(u, - u,). The former choice corresponds to a monovelocity (u,) pulsed source (t = 0) located at x = x,, while the latter corresponds to a steady source emitting particles for t Z 0 with velocity us at point x,. The solution of eq. (31) has been effected by discretising time and velocity variables. We use a constant time step S and assume that f(0, u,, t) varies linearly in a time step, i.e.
m
ux>0 = M-l(%)(L - 4 +f,(u,)(t- ~,-dl&
t,_,
wheref,(u, >= IV, u,, n8 ). Further, the integration over velocity variable (rhs of eq. (31)) is replaced by a high order Gauss-quadrature formula (order J) with weight function exp(-ut/2) over (--03 < uX5 0). With these approximations, eq. (31) can be rewritten as (15 i 5 J)
fi(“i)
=
‘nC”i?‘s,
us) + 6 i: 5 m=l j=l
x *dt’((l - t’)f,(uj)
ujwj exp(uT/2)
+t’f,_,(uj)}G,(O,
vi, (n - m + t’)S; 0, uj, 0)
,
I 0
(39)
where ui and wi are the quadrature directions and weights respectively and S, is the contribution from the source term Q, . For the pulsed source case S, is given by
S.V.G. MENON
74
et al.
(40) while for the steady
source
case 1
S,(U,> x,, us) =
S i j- dt' G,(O, u,, ???=I 0
(n - m + t’)S; x,,
us,
0)
(41)
The integration over the variable t’ in eqs. (39) and (41) is carried out by Simpson’s rule to a specified accuracy. Rewriting eq. (39) in matrix notation we get n-1 (I-
C,,lf,
=
s,
+
??I=1 II m m c
LB
f
-1
+
cn-mLJ+ 4.x-L 7
wherefn and S, are Z-dimensional vectors the .Z x .I matrices B, and C,,, are given
From
and Z is a unit matrix. by
lBrnlij=
1dt'
t’G,(O, u,, (m + t’)6; 0, uj, 0) ,
[C,Jij =
1dt'
(1 - t’)G,(O, u,, (m + t’)6; 0, u,, 0).
eq. (42) it is clear thatf,,
can be solved
matrix (I - C,) is to be inverted only once Milne extrapolation length can be calculated
n-‘(ns) =
5 l+v; &/i:
I=1
wT = W, exp(ufi2).
u!wT f,C”i>
,=I
recursively.
(42) The elements
It may be noted
since it is independent as
of
that
of II. The
(43)
T
f(0, u,, t) is known
at time
steps t, = nS, one can compute the function f(x, u,, t) or its moments JX(x, t) and E,(x, t) as given by eqs. (23), (25) and (30).
n(x, t),
where
Once
the distribution
7. Results In figs. l-3 we give the boundary distribution f(0, u,, t), the particle n(x, t) and the mean square velocity (ut),, = E,(x, t) ln(x, t) for a steady (switched at t = 0) at time t = 500. The same parameters calculated pulsed source at large t (-100) were almost indistinguishable and hence
density source for the are not
GREEN’S
FUNCTION
APPROACH
TO FOKKER-PLANCK
EQUATION
75
0.7r
vxFig. 1. Normalised boundary v, = -1.0) at t = 500.
f(0, v,, t) vs v, for steady
distribution
001 0.0
1
1
1.0
2.0
3.0
source
(with x, =20.0
and
I
LD
X--,
Fig. 2. Normalised particle q= -1.0) at t=500.
density
n(x, t) ln(0, t) vs x for the steady
source
(with x, = 20.0 and
repeatedly given. The boundary distribution and density plots are normalised such that n(0, t) = 1.0. In the numerical calculations we have used x, = 20.0, u, = -1.0 and a 15th order Gauss-quadrature formula16) (for the velocity variable integration) for both cases. We have also used the time steps S = 0.1 and 6 = 0.5 for the pulsed source and steady source respectively. It was found necessary to use a high accuracy (-10w4) Simpson rule for time integration in computing the source terms in eq. (41) and the matrices B, and C,. Further reduction of time steps was then found to have negligible effect on the boundary distribution and other parameters. This is expected to be due to the linear variation of f(0, u,, t) assumed in a time step. It can also be shown that f(0, uX, t)
76
S.V.G.
MENON
et al.
1.6
Fig. 3. Mean square
velocity
(ut),,
vs x for steady
source
varies as tm3” for the pulsed
source
source for sufficiently large We note that the boundary lian. For positive values of zero since K = 0. However, eq. (33) for positive u,. The
values of t. distribution in u, the boundary we have plotted slight departure
mainly
due
to the
finite
order
and becomes
(with x, = 20.0 and u, = ~ 1.0) at
independent
I = 500.
oft for the steady
fig. 1 resembles a shifted Maxweldistribution should be identically the numerical values of the rhs of from zero for U, > 0 is found to be
of quadrature
(15th)
employed
for velocity
discretisation. We observe that these deviations are much smaller than those obtained by Burschka and Titulae?) using EEM of 140th order (compare fig. 1 of this paper with fig. 2 of ref. 3). We note that Titulae?) and Mayyah) had shown that the boundary distribution approaches zero as u, + 0. However, the EEM “) and a numerical simulation’) lead to non-zero values (-0.186) while the present calculation boundary distribution and the present
method
yields a value (-0.026) much closer to zero. In fact the for small u, as computed from numerical simulation’) are somewhat
different
corresponding formula given
to numerical simulation have in ref. 7 with the normalisation
unity. These accurate.
findings
suggest
that
GFM
as seen in table
I. The values
been obtained using the fitted that the area under the curve is
results
presented
here
are
more
In table II we give the numerical values of n(0, t) lJ,(O, t), (vf),, at x = 0 and Milne extrapolation length K’(t) for the pulsed source (at t = 100) and steady source (at t = 500). The steady source case for sufficiently large time should closely approximate the time independent FPE and therefore these results can be compared with those available in the literature”.“‘). Values of these parameters obtained by EEM of 140th order and extrapolated infinite order”,‘) and the numerical simulation’) are also given in table II for comparison. For
GREEN’S
FUNCTION
APPROACH
Boundary Numerical (ref. 11)
Velocity (u,) 0.0
TABLE I distribution
f(0, u,,
simulation
0.01 0.02 0.05 0.1 0.2 0.3 0.5 0.7 1.0 1.4 2.0 3.0 4.0
length
and boundary
EQUATION
t).
0.0260 0.0330 0.0886 0.1309 0.2087 0.2935 0.4093 0.4900 0.5875 0.6215 0.5843 0.4378 0.1960 0.0210 0.0008
TABLE II values of density
and mean
square
Type
n-‘(t)
n(0, t) /J,(O, I)
(u?),,
pulsed source at t = 100
1.459
0.937
1.556
steady source at t = 500
1.461
0.936
1.561
1.438
1.054
1.365
time independent FPE (extrapolated infinite order EEM “)
1.461
0.914
1.590
time independent FPE (extrapolated infinite order EEM-revised’)
1.460
0.951
1.535
simulation
1.454 * 0.004
0.942 ” 0.008
1.544
time independent
FPE
77
Steady source at r = 500
0.1857 0.1869 0.1976 0.2094 0.2439 0.2984 0.3958 0.4764 0.5852 0.6274 0.5904 0.4363 0.1920 0.0206 0.0007
0.001
Milne extrapolation
TO FOKKER-PLANCK
velocity.
at x = 0
(140th order EEM ‘))
results’)
instance, our value of 1.461 for A-’ (steady source) is to be compared with 1.438 and 1.461 obtained in the EEM of 140th order and extrapolated infinite order, respectively. Table III indicates the rate of convergence of our results with t for the steady source case. It should be remarked that there are significant differences between the integral parameters of the pulsed source and steady source cases (table II) even at large times. This can be understood from fig. 4 where we have plotted the current JX(x, t)lJx(O, t) for these two cases. For the steady source case for
78
S.V.G.
MENON
TABLE
Mime
Time
et al
III
extrapolation length and boundary values of density and mean square velocity as a function of f for the steady source. (t)
5 10 15 20 30 50 100 200 300 400 500
‘K’(t)
n(0, t)iJX(O, t)
(vf).,
3.0745 1.8903 1.6603 1.5762 1.5139 1.4507 1.4660 I .4620 1.4611 1.4608 1.4607
0.3606 0.6870 0.8093 0.8605 0.9006 0.9226 0.9324 0.935 1 0.9357 0.9359 0.9360
8.5261 2.7514 2.0517 1.8317 1.6810 1.6050 1.5722 1.5634 1.5616 1.5609 1.5606
at x = 0
1.021.01-
-2 3
0.96 -
Fig. 4. Normalised current profiles J,(x, r)lJ,(O, source at t = 500, (3) 1,(x, f)lJ,(O, t) = 1.
r) vs X, (1) pulsed
source
at
f = 100, (2) steady
sufficiently large times the current is almost flat spatially. Thus this case approximates the time independent FPE. However, as the current profile is quite different for the pulsed case and therefore its results cannot be compared directly with those of the time independent calculations.
8. Conclusions
In this paper we have formulated the GFM for solving the time dependent FPE with an absorbing boundary. The method consists of obtaining a representation for the distribution function (satisfying the BC) in terms of the boundary
GREEN’S
FUNCTION
APPROACH
TO FOKKER-PLANCK
EQUATION
79
distribution which should be determined by solving an integral equation. The 3-D velocity dependence of the distribution function has been incorporated in the formalism. It has been shown that for calculating the velocity moments of general interest (density, current and energy density) it is not necessary to solve for the full distribution function; it is enough to determine three reduced distribution functions (having one velocity component only). We have presented a numerical scheme to solve the integral equations for the boundary distributions. Numerical results have shown that this procedure yields results which are far superior to that obtained by the conventional EEM and in fact are closer to the revised value?). This is due to the fact that EEM results converge very slowly as the order of approximation is increased. Since the infinite medium GF for the time independent FPE is not available in closed form, we have treated this problem in the time dependent formalism using a steady source of particles away from the boundary. One of the drawbacks of the present method is that the infinite medium GF for the problem under consideration should be known. For Brownian motion under general external force fields this is not known; however, for the case of quadratic potentials (in space variables) the Green’s functions are available in literature14’15 ). The present formalism can also be extended to solve the FPE in spherical geometry for treating the condensation problem of a spherical droplet”).
References 1) 2) 3) 4) 5) 6) 7) 8) 9) 10) 11) 12) 13) 14) 15) 16)
S. Harris, J. Chem. Phys. 75 (1981) 3103. K. Razi Naqvi, K.J. Mork and S. Waldenstrom, Phys. Rev. Lett. 49 (1982) 304. M.A. Burschka and U.M. Titulaer, J. Stat. Phys. 25 (1981) 569. Y.S. Mayya and D.C. Sahni, J. Chem. Phys. 79 (1983) 2302. U.M. Titulaer, J. Stat. Phys. 37 (1984) 589. Y.S. Mayya, J. Chem. Phys. 82 (1985) 2033. J.V. Selinger and U.M. Titulaer, J. Stat. Phys. 36 (1984) 293. S. Waldenstrom, K.J. Mork and K. Razi Naqvi, Phys. Rev. A 28 (1983) 1659. S. Waldenstrom, K. Razi Naqvi and K.J. Mork, Ark. Fys. Semin. Trondheim, No. 8, ISSNO 365-2459 (1983). D.C. Sahni, Phys. Rev. A 30 (1984) 2056. Vinod Kumar and S.V.G. Menon, J. Chem. Phys. 82 (1985) 917. S. Harris, J. Chem. Phys. 72 (1980) 2659. A 32 (1985) 3832. S.V.G. Menon and D.C. Sahni, Phys. Rev. S. Chandrasekhar, Rev. Mod. Phys. 15 (1943) 1. H. Risken, The Fokker-Planck equation; Method of Solution and Applications, Springer series in Synergetics, H. Hanken, ed. (Springer, Berlin, 1983). N.M. Steen, G.D. Byrne and E.M. Gelbard, Math. Comp. 23 (1969) 665.