Physics Letters A 179 (1993) 385-388 North-Holland
PHYSICS LETTERS A
Reaction probability derived from a stochastic interpolation formula Tetsuya Misawa
Department of Mathematics, Collegeof GeneralEducation, Nagoya City University,Mizuho-ku, Nagoya467, Japan Received 19 April 1993; revised manuscript received 16 June 1993; accepted for publication 23 June 1993 Communicated by J.P. Vigier
In a simulation involving a diffusion-controlled reaction with discrete time-steps, lack of information regarding the trajectory may result in a failure to count the number of reactions of the particles within the step. In order to rectify this, an interpolated diffusion process is used. The process is derived from a stochastic interpolation formula recently developed by the author. In this method, the probability that a reaction has occurred during the time-step given the initial and final positions of the particles is calculated. The result confirms an improvement over Clifford and Green's work on the same problem.
It is well-known that m a n y reactions in solution dealt with in radiation chemistry are diffusion-controlled [ 1,2 ]. Such reactions are d e t e r m i n e d by the rate at which the reaction particles encounter each other in the course o f diffusion. Hence, in chemical reaction models, diffusion processes are i m p l e m e n t e d to determine the trajectories o f the particles, which are often assumed to be governed by the following ( o n e - d i m e n s i o n a l ) stochastic differential equation [ 3 ], d X ( t ) = b ( X ( t ) , t) dt+tr(X(t), t) d W ( t ) ,
X(0) =Xo,
(1)
where W ( t ) is a Brownian m o t i o n ( W i e n e r process). In this framework, reactions are c o m b i n e d with reactive b o u n d a r y conditions - either an absorbing condition or an elastic one. However, it is difficult to find a solution to eq. ( 1 ) under such conditions analytically, and hence, one must perform a numerical simulation for the stochastic equation by using a finite difference a p p r o x i m a t i o n to eq. ( 1 ) as follows [4],
Xn+l = X , +b(X~, nSt) S t + cr(Xn, nSt)Tnx/'~
( n = 0 , 1, ...),
Xo = X o ,
(2)
where 8t is the discrete time-step, X , = X ( n S t ) and 7~ a r a n d o m variable drawn from a normal distribution N(0, 1 ). Through eq. (2) an a p p r o x i m a t e trajectory for the true solution to eq. ( 1 ) with a b o u n d a r y condition is built up over successive time steps. On the basis o f eq. ( 2 ) with an absorbing boundary, Clifford and Green simulated reactions in radiationinduced spurs [ 1,2]. Their simulation is outlined as follows: The evolution o f interparticle distances in the spurs is described by eq. ( 2 ) . If the distance o f a pair o f particles in the spur is within the reaction distance during the course o f their r a n d o m flight, a reaction occurs and the particles are removed. Thus, the reaction distance plays a role as an absorbing b o u n d a r y for the diffusion process o f the interparticle distance. However, in this framework, one m a y fail to count the n u m b e r o f reactions which would have occurred within the timestep, since we only know the positions o f the particles at the beginning and the end of each step. This failure can be reduced by decreasing the size o f the time-step St; but this is exceedingly c o m p u t e r - t i m e consuming, which is i m p o r t a n t in the simulation m e n t i o n e d above. Therefore, Clifford and Green p r o p o s e d a m e t h o d o f interpolating the initial and final positions o f the particles by a Brownian-bridge process [ 1 ]: By using the interpolation process, they calculated the probability that a reaction has occurred at the absorbing b o u n d a r y 0375-9601/93/$ 06.00 © 1993 Elsevier Science Publishers B.V. All rights reserved.
385
Volume 179, number 6
PHYSICS LETTERS A
23 August 1993
during the time-step, and thereby they counted the number of reactions that would have been missed in the simulation. This idea is very appealing but the method seems to be a bit rough. Indeed, the connection between the original diffusion process ( 1 ) and the Brownian-bridge process is ambiguous. Recently the author has derived a new interpolation formula for the diffusion process satisfying eq. ( 1 ) [ 5 ]. The formula is based on a pinned process which is defined by using the original diffusion process ( 1 ) in a way analogous to that in the derivation of the Brownian-bridge from the standard Brownian process. Such a pinned diffusion process must yield a better interpolation to the original process that the Brownian-bridge one. Therefore, one may expect to improve Clifford and Green's reaction probability by using this formula, and this is the main purpose of this paper. The following is a brief outline of the author's interpolation formula. Let X(t) be a one-dimensional diffusion process satisfying eq. ( 1 ). As the definition of a pinned Brownian bridge [3], we define a pinned process J((u) for X(t) on a given time interval [s, t] as follows, )((u) = t - u x +
t--S
u-s t--s -
t-u [--S
u-s [--s
-oo
--v+--[X(u)-X(s)]---[X(t)-X(u)],
(3)
where x and y are given values on s and t satisfying X ( t ) = y and X(s) = x , respectively. Moreover, in view of the interpolation, it is assumed that the increments X ( u ) - X ( s ) and X ( t ) - X ( u ) in eq. (3) are forward and backward, respectively [5]. Then, if t-s<< 1 holds, X ( u ) - X ( s ) is approximated by'
X(u) - X ( s ) ~b(X(s), s ) ( u - s ) +a(X(s), s)[ W(u) - W ( s ) ] ,
(4)
since d X ( t ) in eq. ( 1 ) arises as a limit of the finite forward increment AX(t) = X ( t + A t ) - X ( t ) ( A t > 0 ) [3 ]. On the other hand, under weaker conditions on b, cr and the probability density function p of X [ 6 ], the diffusion process X(t) satisfies the following backward stochastic differential equation [6-9 ], d . X ( t ) =b.(X(t), t) d.z+a.(X(t), t) d. I 4 , ( t ) ,
(5)
where 1 0
b.(x, t)=b(x, t ) - p~x (a2P) '
~(x, t ) = c r . ( x , t) .
(6)
In eq. (5), W. is a backward standard Wiener process identified through [6] d. ~ , = d W +
10 p0~x (pa) dr.
(7)
Since d . X ( t ) in eq. (5) arises as a limit of the finite backward increment A X ( t ) = X ( t ) - X ( [ - A t ) [ 6 - 9 ] , the backward increment X ( t ) - X ( u ) in eq. (3) is approximately given by'
X(t) - X ( u ) ~b.(X(t), t ) ( t - u ) + a . ( X ( t ) , t) [ D , ( t ) - H,,(u) ] .
(At>0) (8)
Substituting eqs. (4) and (8) into eq. (3) together with eqs. (6) and (7), we obtain the interpolation formula (its approximate version for t-s<< 1; for the exact version, see ref. [7] ) for the original diffusion process X on Is, t] as follows,
f((u)= t - U x + --v+t_s, ( t - u ) ( u - s ) ( b ( X ( s ) ' s ) - b ( X ( t ) ' t ) + ~ 2 r + t-u a(X(s),s)[W(u)-W(s)]t-s
u-S~(x(t),t)[W(t)-W(u)] t-s
, t)) ,
-~
(9)
(The preceding method dealing with the interpolation process )((u) may be a particular case of the so-called "Bernstein diffusions" where arbitrary boundary probability densities are given [5,10]. ) 386
Volume 179, number 6
PHYSICS LETTERSA
23 August 1993
Now, by using eq. (9), we will find the probability that the first passage to a reactive boundary (absorbing boundary), at the separation R, occurs during the time interval t - s , given the positions X(s) and X(t) for the original process X, which would be obtained if the boundary were not present. For this purpose, we rewrite formula (9) as follows: First, for each u on [s, t], the terms containing the differences W ( u ) - W(s) and W ( t ) - W(u) in eq. (9) are replaced by a random variable determined by
t - u a 2 ( X ( s ) , s ) + u-~sa2(X(t), t) ( t - u ) B I--S t--S
(10)
where B is another Wiener process, due to the Gaussian nature of the Wiener process W. Next, we assign x=X(s) and y=X(t), since the values of X at s and t are given, respectively, and rewrite u as s+ vSt, where 8t=_t-s ( << 1 ) and v - ( u - s ) / ( t - s ) (0~< v~< 1 ). Then, eq. (9) is put into the following form,
X(s+ vSt)= ( 1 - v)x+ vy+ v ( 1 - v)St(bl-b2 +c) +x/(1-v)~+v~St(1-v)B
(1-v)St
'
0~
(11)
where bl ~ b(x, s), b2 = -- b(y, t), ffl =-a(x, S), ~2=-if(y, t) and c= (~(O~/Ox) ) (y, t). Equation ( 1 1 ) corresponds to eq. (10) in Clifford and Green's paper [ 1 ], and hence in a way similar to that in their work, we change the variables v and 2 in eq. ( 11 ) as P
r = ( 1 - v)St'
0~
(12)
Z(z) = (1 + r S t ) ) ( ( s + vfit) - x - (r~St)R.
(13)
The probability we seek is that 2 reaches the reactive boundary at the separation R during the time interval St. This is equivalent to the probability that the process Z ( r ) with Z(0) = 0 reaches R - x during the span in time because of transformations (12) and (13). Such a probability is calculated by using a solution of the following Fokker-Planck equation for Z ( r ) , which is derived from eqs. ( I 1 )- (13), with an absorbing boundary condition p ( R - x ) =0,
Op _ Or
O 62 O2p Oz (b'a)+ ~ Oz~,
(14)
where
[(a2/al)2-1]St
~) :11(]18'+ rSt) ) (y / ~ ( z , r ) = 2 [ l + r ( a 2 / a ~ ) 2 8 t ] ( l + r S t ) z + ( 1 - 211 q-r~((aa~/:)~, +(1 -
r[ (a2/o'l)2_ 1 ]St "] (bl - b 2 + c ) (St) 2 (l+rSt)2 ,
2[l+r(a21a~)28t] ] #(r) = N/1
q-r(o'2/rY
1+ r S t
I
R)St (15)
)28t cr~at.
(16)
It is difficult, however, to solve eq. (14) with the boundary condition exactly, since/~ and # depend on the variable r. Hence, as the first approximation for/~ and 6, we take a long-term average for them in the following form, T
/~= lim ~1 I t ) ( z , r ) d r = ( y - R ) S t , T~oo
0
T
•=
l i m ~if O(r) d r = ~
T~oo
St.
(17)
1 ,/ 0
387
Volume 179, number 6
PHYSICS LETTERS A
23 August 1993
By replacing/~ a n d ~ in eq. ( 14 ) b y / ~ a n d 0., respectively, we find an a p p r o x i m a t e f u n d a m e n t a l solution ( t r a n sition p r o b a b i l i t y density ) with a b s o r b i n g b o u n d a r y p (R - . r ) = 0 a n d initial c o n d i t i o n Z ( 0 ) = 0 as p(z,:,0)=q(r,c,
0)-exp
(2(y-R)(x-R)) o1a28/
q(r,z, 2(R-.v)),
(18)
where 1 ( q(r, z, ,~) = 2 x / 2 ~ l a2 ( 8 t ) -; r exp -
[:-5-
r(r-
R)St]2"~
20.102(~l)2r
/] .
( 19 )
By using ( 1 8 ) with ( 1 9 ) [ 1 1 ], we finally get the reaction p r o b a b i l i t y that Z ( r ) reacts at R - x d u r i n g the whole time, that is, A" reacts at R in the t i m e interval 8 / = t - s u n d e r the c o n d i t i o n s ,~(t) =3, a n d )((s) = x as follows, (2(),-R)(.r-R)) 14"= exp
-
0.1 0"2 5 [
(20) "
This is the desired result. As m e n t i o n e d in the i n t r o d u c t o r y part o f this paper, Clifford a n d G r e e n o b t a i n e d the following reaction probability [ l 1. Wc(~=exp(-
2(),-R)(x-R)~~28t J"
(21)
where 0" is a single value o f the diffusion coefficient a ( x , t) in eq. ( 1 ) at initial t i m e s. In contrast, eq. ( 2 0 ) c o n t a i n s two values of 0"(x, t) at initial t i m e s a n d final t i m e t. T h i s indicates that the present result is m o r e precise t h a n Clifford a n d G r e e n ' s . In the future, eq. ( 2 0 ) s h o u l d be c o m p a r e d with eq. (21 ) through a concrete numerical simulation.
References
[ l ] P. Clifford and N.J.B. Green, Mol. Phys. 57 (1986) 123. [ 2 ] P. Clifford, N.J.B. Green and M.J. Pilling, J.R. Star. Soc. B 49 ( 1987 ) 266. [ 3 ] N. lkeda and S. Watanabe, Stochastic differential equations and diffusion processes, 2nd Ed. (North-Holland/Kodansha, Amsterdam/ Tokyo, 1989 ). [4] W. Riimelin, SlAM J. Numer. Anal. 19 (1982) 604. [5]T. Misawa, J. Math. Phys. 34 (1993) 775. [6] U.G. Haussmann and E. Pardoux, Ann. Probab. 14 (1986) 1188. [7 ] E. Nelson, Dynamical theories of Brownian motion (Princeton Univ. Press, Princeton, NJ, 1967 ). [8] J.-C. Zambrini, Int. J. Theor. Phys. 24 (1985) 277. [9] P. Garbaczewski, Phys. Lett. A 143 (1990) 85. [ 10] J.-C. Zambrini, J. Math. Phys. 27 (1986) 2307. [ 1I ] H. Riskin, The Fokker-Planck equation: methods of solution and applications, 2nd Ed. (Springer, Berlin, 1989 ).
388