Monte Carlo Methods and Appl., Vol.3, No.4, pp. 313-326 (1997)
A MONTE CARLO APPROACH TO THE SMOLUCHOWSKI EQUATIONS FLAVIUS GUIAS Abstract. Convergence of a family of Markov jump processes evolving according to coagulation-fragmentation dynamics to the corresponding deterministic equations is proved on every nite interval of existence. After obtaining local existence by the Banach xed-point theorem, we can prove global existence using the positivity of the limit function on the local existence interval.
0. Introduction The so-called Smoluchowski equations (coagulation-fragmentation equations): 1 1 X d u1 (t) = ?u1 X 1 u + 1u dt =2 =1 1 1 ?1 X X d u (t) = 1 X ? u + u u ? u dt 2 =1 ? =1 = +1 ;i
i
i
i;
i
i
k
k
i;k
i
i
i
k
i
k
k;i
i
i
i
k
u (0) = u0 0; k
k
?1 X
u ? 21 u; =1 k = 2; 3; : : : for all k 0 (1) k
i;k
i
k;i
k
i
are used in several elds such as chemistry, meteorology or astronomy, in order to describe the evolution of a system of particles with a size parameter (e.g. mass) k = 1; 2; : : : which can aggregate or break-up. Here u (t) stands for the mean concentration of k-fold particles, and the four terms stand for : () Production of particles of size k by coalescence of particles of sizes i and k ? i, (i < k). ( ) Disappearance of particles of size k by coalescing with particles of an arbitrary size i = 1; 2; : : : . ( ) Production of particles of size k by splitting o from particles of size i, (i > k). () Disappearance of particles of size k by splitting into particles of sizes i and k ? i, (i < k). Used originally by Smoluchowski (1917) for the case of pure coagulation with constant kernels, equation (1) was analyzed e.g. in [1] and [8] and was extended to the case of continuous size parameter e.g. in [7] and [9]. k
1991 Mathematics Subject Classi cation. 60J75, 65C05. Key words and phrases. In nite Systems of ODE's, Markov Jump Processes, Monte Carlo Methods. The author wishes to thank W.Wagner and K.Sabelfeld for their valuable suggestions and comments on the manuscript. 1
2
FLAVIUS GUIAS
Altough all these descriptions of coagulation-fragmentation dynamics are made by deterministic equations involving the mean concentrations, it is very natural to consider eective particles which aggregate or break-up according to stochastic dynamics. The question is now to show the consistency of the two models: stochastic and deterministic, i.e. to prove a convergence result of a density function derived from the particle system to the solution of (1), w.r.t the parameter of the family of stochastic models. For the case of nite-dimensional systems of ODE's, results in this sense can be found in [4] and [5]. In [6] is considered a system of coagulating brownian spheres, whose density function is shown to converge in the Boltzmann-Grad limit (number of spheres radius of a sphere = const) to the solution of a coagulation equation with diffusion, while in [3] a result concerning stochastic particle systems which undergo coagulations with constant kernels is obtained. All this convergences are of stochastic type: in probability (in [6] ), and in mean (in [3]). In our approach we consider particle systems with stochastic coagulation and fragmentation, whose density functions can be described in terms of a family of Markov jump process. A global existence result for equation (1) for some particular classes of unbounded kernels can be found in [8], but essentially the same estimates used for the dierence of the solution and the stochastic process can be used to give a local existence proof based on the Banach xed-point theorem (for the case of bounded kernels). We will prove afterwards mean-square convergence of the family of processes to the solution of (1) on every nite existence interval. This yields positivity and uniqueness of the solution, which will lead nally to a global existence result. For the case of monodisperse initial conditions, we will give also estimates of the order of convergence in the mean square norm.
1. The local existence result Consider the scale of Banach spaces fX g, q 2 [0; 1], endowed with the norms kak = P1=1 k ja j, for every a = (a )1=1 2 X . Consider the in nite matrices with positive elements = ( )1 =1 , and = ( )1 =1 , which are supposed to satisfy the following hypotheses: (i) is symmetric, i.e. = , for all i; j = 1; 2 : : : (ii) is uniformly bounded, i.e. there exists a constant M , such that M , for all i; j = 1; 2 : : : (iii) = 0 for i j and = ? for i > j . Hypothesis (iii) and the fact that the expected size of a particle after splitting cannot exceed its initial size, imply: q
q
q
k
k
q
k
k
i;j
i;j
i;j i;j
i;j
j;i
i;j
i;j
i;j
i;i
?1 X
k
=1
i
k;i
j
= k2
?1 X
k
=1
i
k;i
B
i
for any k > 1, with B a positive constant (see also [8]). On the spaces X de ne the bilinear operator [; ] and the linear operator L by: q
A MONTE CARLO APPROACH TO THE SMOLUCHOWSKI EQUATIONS
[f; g] (Lf )
= 12
k
k
=1
? f g ? i;k
i
i
k
i
1 X
=
k
(X ?1
=k+1
i;k
i
i
?f
k
=1
k
i
g ?g i
k;i
k
k;i
=1
f
1 X =1
f k;i
)
(2)
i
i
i
?1 X
f ? 12
1 X
3
(3)
k
i
for any f; g 2 X . Using the product formula for series, the boundedness of and the inequality (i + j ) i + j , one gets for all q 2 [0; 1]: q
q
q
q
k[f; g]k M (kf k0 kgk + kf k kgk0) 2M kf k kgk q
q
q
q
and by the hypotheses on and from:
1 1 X X k =1 = +1 q
i
k
k
i;k
1 X 1 1 X 1 X X i k jf j = f =1 =1 =1 =1 1 1 ?1 X X X q
i
i
k
follows
k jf j
=1
k
k
=1
i
k;i
B
k;i
jf j = k
k
=1
1 X ?1 X k
k
=1 i=1
k jf j = B kf k q
kLf k 32B kf k q
q
i
k
q
k
i;k
i
k
q
(4)
q
i
q
k;i
jf j k
(5)
(6)
q
for any f; g 2 X and q 2 [0; 1]. With these notations, equation (1) can be written as: q
@ u = [u; u] + Lu u(0) = u0 t
or in the equivalent integral form:
Z
t
u(t) = u(0) + f[u(s); u(s)] + Lu(s)g ds: 0 For an arbitrary T > 0, let X = C (0; T ; X1) be the Banach space of continuous functions de ned on the interval [0; T ] with values in X1 , endowed with the norm T
kf k T = sup kf (t)k1;
X
t
T
for all f 2 X : T
If we de ne on X the operator A by: T
Af (t) = u0 +
Z
t
0
f[f (s); f (s)] + Lf (s)g ds;
the coagulation-fragmentation equation on the interval [0; T ] can be then written as:
u = Au
in X : The local existence result is stated in the next theorem. T
(7)
4
FLAVIUS GUIAS
Theorem 1.1. De ne Q = max(4M; 3B=2) and
(w) = Q(w + ku k +w1)(w + ku k ) 0 1 0 1
for any w 0. If T = (r) for some r > 0 and u0 6= 0, then equation (7) has a unique solution u 2 B (u0 ) = fu 2 X : ku ? u0k T rg. r
T
X
Proof
For u 2 B (u0 ) we have: r
kAu ? u0 k T = sup kAu(t) ? u0 k1 = sup k
X
Z
sup
t
t
0
T
t
0
T
0
f[u(s); u(s)] + Lu(s)g dsk1
2M ku(s)k21 + 32B ku(s)k1 ds
t
sup
T
t
fk[u(s); u(s)]k1 + kLu(s)k1g ds
Z
T
t
t
Z
T 2M sup ku(t)k21 + 32B sup ku(t)k1 = T 2M kuk2 T + 32B kuk
t
T
t
T
2M (ku0k1 + r)2 + 3B (ku0 k1 + r)
T 2 TQ (r + ku0k1 + 1) (r + ku0k1 ) = r;
X
X
T
the last equality being our assumption in the hypothesis. Thus the operator A leaves the ball B (u0 ) invariant. In order to obtain existence and uniqueness by the Banach xed-point theorem, we will show next that A is a contraction on B (u0 ). For u; v 2 B (u0 ) we have: r
r
r
kAu ? Avk
X
T
= sup k
t
T
= sup k
t
T
sup
t
T
sup
t
T
Z
t
f[u(s); u(s)] ? [v(s); v(s)] + L(u(s) ? v(s))g dsk1 =
t
f[u(s) ? v(s); u(s) + v(s)] + L(u(s) ? v(s))g dsk1
0
Z
0
Z
t
0
fk[u(s ? v(s)); u(s) + v(s)]k1 + kL(u(s) ? v(s))k1 g ds
Z t
0
2M ku(s) ? v(s)k1 ku(s) + v(s)k1 + 32B ku(s) ? v(s)k1 ds
T 2M ku + vk T + 32B ku ? vk T 3 B T 2M (kuk T + kvk T ) + 2 ku ? vk T TQ (ku0 k1 + r + 1) ku ? vk T = r + kru k ku ? vk X
X
X
X
X
X
0 1
X
T
the last equality being again a consequence of our hypothesis. Since u0 6= 0, the operator A de nes a contraction on B (u0 ) and by the Banach xed-point theorem r
A MONTE CARLO APPROACH TO THE SMOLUCHOWSKI EQUATIONS
5
it has a unique xed point u 2 B (u0 ), which is the local solution of the coagulationfragmentation equation. r
2
Remark It is not trivial to obtain a global solution by applying succesively the previous result. This would be possible if one could nd a sequence fr g, such that k
the series for the correspondent time intervals diverges, which is equivalent with the divergence of the series 1 X k
=0
P r 2 : k
=0 rj
k j
But for r = g(k + 1), for any monotone decreasing continuous function g, a comparison with the Riemann integral of k
g(x)=(
Z
x
1
g(s)ds)2
shows that the series converges. However, if the solution of equation (1) would be positive on the existence interval, one can show that it has the mass conserving property, which nally leads to the existence of a global solution. The positivity will turn out as an immediate consequence of the approximation result of the solution by a family of Markov jump processes proved in the next section. 2. Stochastic modelling of coagulation-fragmentation dynamics Consider in a closed volume a nite number of particles of sizes 1; 2 : : :; k; : : : which undergo coagulations and fragmentations at rates governed by the kernels and . Let n (t) denote the number of k-fold particles at time t, u (t) = h n (t) be the corresponding mean concentration, where h 1=V is a parameter inverse proportional volume. Pto the considered If N = k n (t) is P the total number P of monomers (which is conserved), the total mass is given by: k u (t) = k h n (t) = h N; thus we have also h 1=N . In order to simplify things, we take from now on the total mass equal to 1 in both cases: the Smoluchowski equations (i.e. ku0 k1 = 1), and the stochastic particle model. The maximal possible size of a cluster being N , let us denote by e the element in R with the k-th component equal to 1 and all others equal to 0. The transitions of the density process are then the following: k
k h
k
k
k h
k
k
N
(for j i), at rate
K
ij
u ! u ? he ? he + he + (coagulation) h
h
i
j
i
j
= h n (t)n (t) ? h2 n (t)n (t) + n (t) = h1 u u ? 21h u u + 21 u i;j
i;j
i
i h
j
ij
i;j
j
h
ij
i h
i;j
i
j
i;j
j
i;j
h
i h
u ! u ? he + he ? + he (fragmentation) h
(for j [i=2]), at rate
h
i
F = ij
i
i;j
j
j
n (t) = h1 i
i;j
u
i h
i
6
FLAVIUS GUIAS
We made here the natural assumptions that the coagulation rates of the density process are directly proportional to the numbers of the corresponding particles (proportional to n (n ? 1)=2 if the particles are of the same type) and inverse proportional to the volume, and that the fragmentation rates are directly proportional to the number of respective particles. u (t) is thus an R -valued Markov jump process with in nitesimal operator acting on the bounded, measurable functions on R by: i
i
N
h
h
N
( g)(u) = (u) h
where h
(g(v) ? g(u)) (u; dv) h
1 0 X X X F = 12 @ K + F A + K K + 6= 1 0 [ 2] X X X 1 @ uu + u A? 1 u X
ij
ij
j
RN
h
X
(u) =
Z
i
j
= 2h
i
i;j
j
2
j
ii
i
j
i
i
i;j
i;j
ij
ij
j
i=
i
i;i
i
denotes the parameterPof the exponentially distributed waiting time in the state u 2 E = fv 2 R+ : k v = 1g, the values of outside E being de ned as 0; while the discrete measure N
N
k
h
k
N
K = (u)
on fu ? he ? he + he + g F = (u) on fu ? he + he ? + he g :
(u; dv) puts mass h
ij
h
i
j
ij
h
i
i
i
j
j
j
Details concerning this approach in de ning Markov jump processes can be found e.g. in [2], p.333. From the martingale characterization of Markov processes (see [4], p.162) known also as Dynkin's formula, it follows that for any bounded, Borel measurable function g de ned on R , the process N
M (t) = g(u (t)) ? g(u (0)) ? g
h
h
Z
h
0
t
( g)(u (s))ds h
h
(8)
is a martingale with respect to the -algebra generated by the Markov jump process u (t). For v 2 R and r > 0, de ne h
N
p (v) = (v ); k r
r
k
(9)
where the functions are de ned on R for all r > 0 by: r
(x) = r
x for jxj r 2r sgn(x) ? r2 x?1 for jxj > r:
p are thus bounded C 1 -functions de ned on R . k r
N
(10)
A MONTE CARLO APPROACH TO THE SMOLUCHOWSKI EQUATIONS
Since u 2 E we have p1 (u ) = u , thus: X p1 (u ) = (?he ? he + he + ) K + h
h
k
k
N
h
k i
k h
h
k j
k i j
X
ij
[ 2]
7
(?he + he + he ? ) F k i
k j
k i j
X X = 21 (?e ? e + e + ) u u + 21 (?e + e + e ? ) j
i
j
k i
k j
k i j
i;j
X ? h2 (?2e + e2 ) u
j
i h
k i
h
i;j
k i
k i
i;i
ij
i=
k j
k i j
i;j
u ? i h
j
i h
i
? ?1 X X X u ? 12 u + u u ?u = 21 =1 =1 = +1 += +h ; where for k = 2k0 + r with k0 > 0 and r 2 f0; 1g, (t) = 12 2 u (t) ? (1 ? r ) u (t) : Hypotheses (ii) and the fact that u 1=k for all k, yield: j (t)j 2M :
X
N
i
j
j
i h
i;j
k
k;j
k;j
h
j
k
N
k
j
k h
h
k h
j
i
i;k
u + i h
k
(11)
k h
k
k
k h
(12)
0
k h
k;k
k h
k0 ;k0
k
k h
k h
(13)
k
Consider u 2 E as an element of X1 by de ning the coordinates of index larger than N as 0. This fact, together with the property that u u = 0 if i + j > N , enables us to consider the indices j and i in the second, respectively the fourth sums in the last equation of (11) as running up to in nity. From (8) we obtain then for all k = 1; 2; : : :; N the representation: h
N
j
i h
Z u (t) = u (0) + [u (s); u (s)] t
k h
k h
h
0
h
k
h
+ (Lu (s)) + h (s) ds + M (t); k
h
k h
k h
(14)
where M (t) is a martingale w.r.t. the -algebra generated by the process. By the same argument as above, one obtains for all k > N : k h
Z [u (s); u (s)] 0 = u (t) = u (0) + t
k h
k h
h
0
h
k
+ (Lu (s)) ds k
h
(15)
If we consider M (t) = (M (t)1=1 as an element of X1 by de ning the components of index larger than N as 0, we have the following result: Lemma 2.1. M (t) 2 X1 satis es the estimate: k h
h
h
k
E sup jM (t)j2 1
C h
=1 1 u
t
T
k h
(16)
T
k
where C is a constant depending only on T; M and B , and where the conditional expectation is taken with respect to the initial value of u . T
h
Proof
Since for all k we have:
M (t) = u (t) ? u (0) ? k h
k h
k h
Z [u (s); u (s)] t
0
h
h
k
+ (Lu (s)) + h (s) ds h
k
k h
8
FLAVIUS GUIAS
(where for k > N the 's are taken as 0), using (13) it follows that: k h
kM (t)k1 ku (t)k1 + ku (0)k1 + T sup 2M ku (t)k21 + 32B ku (t)k1 + hk k1 2 + T R a.s. (17) h
h
h
h
t
h
h
T
where R = 4M + 3B=2. Using a result from [5], one has that the coupled process
Y (t) = (u (t); M (t)) =1 k h
h
N k
k h
is Markov with values in R R and for any bounded function g(u; m) g(m ) of class C 1 which depends only on the k-th coordinate of the second component, the pair N
g(m ); (u) k
h
N
k
Z ? g (v ? u + m ) ? g(m ) ? (v ? u ) @ g(m ) (u; dv) N k
k
k
k
k
R
k
k
k
h
lies in the full generator of Y (as de ned in [4], p.23), that is: h
Z ? g(v ? u (s) + M (s)) ? g(M (s))? RN 0 ?(v ? u (s)) @ g(M (s) (u (s); dv) ds:
E [g(M (t)] = k h
u
k
k h
Z
t
E (u (t)) h
u
k
k h
k
h
h
h
k h
k h
k h
A MONTE CARLO APPROACH TO THE SMOLUCHOWSKI EQUATIONS
9
By taking g(u; m) = (p2+ (m))2 , where the functions p were de ned in (9), we get for all k N : Z h1 X 2 E [jM (t)j ] = E 2 j ? he ? he + he + + M (s)j2 ? jM (s)j2 0 ?2(?he ? he + he + ) M (s) h?1 u (s)u (s) + X j ? he + he + he ? + M (s)j2 ? jM (s)j2 ? + 21 ?2(?he + he + he ? ) M (s) h?1 u (s) ? X ? 21 j ? 2he + he2 + M (s)j2 ? jM (s)j2 k
k r
TR
t
k h
u
k i
u
k j
k i j
k h
k h
i;j
k i
k j
k i j
k h
k i
k j
k j
k i j
k h
k i
k h
i;j
k i j
j
i h
k h
h
k h
j
k i
k i
i
i;j
i
i h
k h
?2(?2he + he2 ) M (s) u (s) ds = k i
Z
k i
k h
i;i
i h
h X E h2 j ? e ? e + e + j2 u (s)u (s) + 0 X + h j ? e + e + e j2 u (s) ?
=
t
k i
u
k j
k i j
?
i;j
i;j
j
i h
h
i;j
2
k i
h2 X
= h
+
k i j
j ? 2e + e2 j2 u (s) ds = k i
h E 12 u
0
X
i;k
k i
X
i
t
i h
i
j
?2
Z
k j
+j =k
i
u (s)u (s) + i;j
j
i h
1X
u (s) + 2 i h
i h
i;i
X
k h
k;k
h
k;j
u (s) ? k h
j