information ... The LaTeX file is compatible with: EmTex version 3.0 [3a], LaTeX version 2.09, TCI LaTeX. The actual LaTeX file (placed in the "BODY" case) has the checksum value =60.362 143 (obtained by the checksum command 'sum' from UNIX System V/386 Release 3.2) ---------------------------------------------- BODY \documentstyle[12pt]{article} \setlength{\topmargin}{0in} \textheight=21.60 true cm \setlength{\textwidth}{14.0cm} \newtheorem{theorem}{Theorem} \newtheorem{lemma}[theorem]{Lemma} \newtheorem{proposition}[theorem]{Proposition} \newtheorem{conjecture}[theorem]{Conjecture} \newtheorem{example}[theorem]{Example} \newtheorem{definition}[theorem]{Definition} \newtheorem{corollary}[theorem]{Corollary} \renewcommand{\baselinestretch}{1} \setlength{\oddsidemargin}{+1.0cm} \setlength{\evensidemargin}{0in} \begin{document} \parindent=0.7 true cm \topskip =1.5 true cm \begin{center} {\bf {ON THE NUMERICAL SIMULATION OF A CLASS OF REACTIVE BOLTZMANN TYPE EQUATIONS}} \vskip 24 pt {C. P. Gr\"unfeld*} \vskip 12 pt {IGSS, Institute of Atomic Physics,\\ Bucharest-Magurele, P.O. Box, MG-6, RO-76900, Romania,\\ E-mail grunfeld@roifa.bitnet, grunfeld@ifa.ro} \vskip 24 pt {D. Marinescu} \vskip 12 pt {Institute of Applied Mathematics,\\ Calea 13 Septembrie, No13, Sector 5, Bucharest, Romania} \vskip 48 pt \end{center} \centerline{\bf {ABSTRACT}} \medskip We study the numerical approximation of solutions to a class of nonlinear kinetic equations for reacting gas mixtures. We first prove the existence and uniqueness of solutions and then provide a time-discretized version of the equations. We finally obtain a probabilistic-convergent, numerical scheme by extending simulation techniques for the solutions of the classical Boltzmann equation. \vspace{24pt} \begin{center} {\bf 1. INTRODUCTION} \end{center} \setlength{\baselineskip}{18pt} Recently, Babovsky and Illner$^{1,\,2}$ have proposed an efficient numerical scheme, consistent with the classical Boltzmann equation. Starting from Nambu's ideas$^3$ they have used the properties of the Boltzmann equation, time discretization (and local space homogenization) to obtain a convenient approximate form of the equation. To provide an algorithm, with small numerical effort, they have introduced numerical simulation. Moreover, they have proved the convergence almost sure, in some sense, of the resulting approximation scheme. In the Nambu-Babovsky-Illner scheme the resulting approximating solutions of the Boltzmann equation appear as suitable (homogeneous) sums of point measures$^1$ defined as follows. Let $\left\{ \delta (x-x_{i,n})\,\mid 1\leq i\leq n,\;n\in {\bf N}^{*}\right\} $ be a family of Dirac measures on ${\bf R}^m$ and \begin{equation} \label{pri}\sigma _n:=\frac{a_n}n\sum_{i=1}^n\delta (x-x_{i,n}),\;\;\;n\in {\bf N}^{*}=\{1,2,...\}, \end{equation} where $a_n>0$, for all $\in {\bf N}^{*}$. Then $\sigma _n$ is called a homogeneous sum of point measures (HSPM) approximating $0\leq $$f\in L^1($ $% {\bf R}^m;dx)$ if $\sigma _n\stackrel{w}{\rightarrow }f(x)dx$ as $% n\rightarrow \infty $, the convergence being in the weak sense of measures (actually, in many applications, $a_n=$ constant $=\int_{{\bf R}^m}f(x)dx$ ). The HSPM approximation techniques seem particularly useful for numerical kinetic theory. Indeed, computational and numerical reasons$^{4-6}$ show that the HSPM approximation methods are very convenient for problems where the solution is a probability measure on ${\bf R}^m$ and where the main interest is to calculate moments as well as expectation values of some random variables with respect to the solution. In particular, the control of the approximation can be made by means of the Koksma-Hlavka inequality$% ^{4,\,5}$ in terms of discrepancy. We recall that if $\mu $ and $\nu $ are finite nonnegative measures on ${\bf R}^m$, then, by definition, their discrepancy $D(\mu ,\nu )$ is$^{2,\,4\,,\,5}$% \begin{equation} \label{dis}D(\mu ,\nu ):=\sup \limits_{a\in {\bf R}^m}\left| \mu (\Lambda (a))-\nu (\Lambda (a))\right| . \end{equation} Here\ $\Lambda (a):=\left\{ x\in {\bf R}^m\;\left| x\preceq a\right. \right\} $, where $\preceq $ denotes the usual order on ${\bf R}^m$ (i.e. $% x=(x_1,...,x_m)\preceq $ $y=(y_1,...,y_N)$ iff $x_i\leq y_i$ for all $% i=1,...,m$). We also recall$^2$ that a sequence of measures $\mu _n$ is said to converge to $\mu $ with respect to discrepancy if $D(\mu _n,\mu )\rightarrow 0$ as $n\rightarrow \infty $. If $\mu $ is absolutely continuous with respect to the Lebesgue measure on ${\bf R}^m$, then it is known$^2$ that the convergence of $\mu _n$ to $\mu $ with respect to discrepancy is equivalent to the weak convergence in the sense of measures. In the present paper, we extend some of the methods of Babovsky and Illner$% ^{1,\,2}$ to obtain a convergent simulation scheme for a general class of kinetic equations, studied in Ref.\thinspace 7,\thinspace 8 modelling a reacting gas, without external forces, composed of distinct species of mass points with (possibly) several reaction partners. The equations of Wang, Chang and Uhlenbeck$^{9,\,10}$ Ludwig and Heil$^{11}$, as well as the Krook and Wu$^{12}$ multicomponent Boltzmann equations can be written in forms belonging to the class considered in Ref.\thinspace 7,\thinspace 8. Our study will be limited to the space homogeneous case, which has its own interest when chemical reactions are present (however the space dependent situation seems also tractable by adapting the local homogenization techniques of Ref.\thinspace 2). The paper is organized as follows. In the next section we introduce the setting, the basic kinetic equations and we formulate the approximation problem. In Section 3 we study the initial value problem for the kinetic equations of Section 2, formulated in a suitable function space. Then we obtain an iterative time-discretized approximation of our equations. In Section 4 we generalize certain probabilistic selection results of Ref.\thinspace 1,\thinspace 2. This is possible due to some clarifications with respect to the nature of the convergence introduced by Babovsky and Illner$^{1,\,2}$: the probabilistic part of the convergence proof of Ref.\thinspace 1,\thinspace 2 is based on the central limit theorem for row-wise i.i.d. random variables and the Borel-Cantelli Lemma. In our case, the argument follows from a simple version of strong law of large numbers for arrays of (not necessarily identically distributed) row-wise independent, random variables, which is straightforward from the Chebyshev inequality and the Borel-Cantelli Lemma. In Section 5, the results of Section 4 are applied to the iterative scheme obtained before. Thus we construct the numerical algorithm for the original Cauchy problem and obtain our main result, namely the convergence of the numerical scheme. In the end of this section it should be recalled that the probability of multi-particle collisions is zero in some sense$^{13}$ in the dynamics of the classical hard sphere gas with elastic collisions (which plays an essential role in the validation of the classical Boltzmann equation).\ Consequently, in certain cases modelling the simple {\it continuous} gas with% {\it \ multiparticle} {\it elastic} collisions would appear somehow artificial. However the situation seems completely different in the case of the reacting gas: reaction processes{\it \ producing} {\it more} {\it than two }particles may be important to the gas evolution. The methods of this paper have been tested on the Krook-Wu multicomponent Boltzmann equation$^{14}$. They may be used to provide numerical solutions to other kinetic equations$^{15}$.\medskip\ \begin{center} {\bf 2. THE SETTING AND THE PROBLEM} \end{center} Consider the space-homogeneous version of the kinetic equations introduced in Ref.\thinspace 7,\thinspace 8 to describe a reacting gas consisting of $N$ distinct species of point particles with one-state internal energy and at most $M\geq 2$ identical partners in each reaction (collision) channel (in APPENDIX we briefly show how these equations arise in modelling the reacting mixtures of {\it real} gases, composed of particles with several internal states). \begin{equation} \label{unu}\partial _tf_k(t,{\bf v})=P_k(f)(t,{\bf v})-S_k(f)(t,{\bf v}% )\quad k=1,...,N, \end{equation} The unknown functions $f_k:\left[ 0,\infty \right) \times {\bf R}^3{\bf % \rightarrow }\left[ 0,\infty \right) $ are the one particle distribution functions describing species $k=1,...,N$ of particles (here $f_k=f_k(t,{\bf v% })$ with $t$-time, ${\bf v}$-velocity). The operators $P_k(f)$ and $S_k(f)$ are expressed as sums of terms, multi-linear in the components of $% f=(f_{1,}...,f_N)$:% $$ S_k(f)(t,{\bf v}) $$ \begin{equation} \label{e2'}:=\sum_{\alpha ,\beta \in {\cal M}}\alpha _k\int_{{\bf R}^{3\mid \alpha \mid }}\delta ({\bf w}_{k,\alpha _k}-{\bf v})\;d{\bf w}\int_{{\bf % \Omega }_\beta }{\bf \,}r_{\beta ,\alpha }({\bf w},{\bf n})\cdot f_\alpha (t,% {\bf w}){\bf \,}d{\bf n,} \end{equation} $$ P_k(f)(t,{\bf v}) $$ \begin{equation} \label{e2}:=\sum_{\alpha ,\beta \in {\cal M}}\alpha _k\int_{{\bf R}^{3\mid \alpha \mid }}\delta ({\bf w}_{k,\alpha _k}-{\bf v})\;d{\bf w}\int_{{\bf % \Omega }_\beta }\,p_{\beta ,\alpha }({\bf w},{\bf n})\cdot f_\beta (t,u_{\beta ,\alpha }({\bf w},{\bf n}))\,d{\bf n.} \end{equation} Here ${\cal M:}=\left\{ \gamma =(\gamma _k)_{1\le k\le N}\mid \gamma _k\in \{0,1,\ldots ,M\},\;2\leq \sum_{k=1}^N\gamma _k\right\} $ designates the set of collision channels (with $\gamma _k$ particles of species $k=1,...,N$ in channel $\gamma $). Further ${\bf \Omega }_\gamma $ is unit sphere in ${\bf R% }^{3\left| \gamma \right| -3}$ (the set of impact parameters) with $\left| \gamma \right| :=\sum_{k=1}^N\gamma _k$. Also $r_{\beta ,\alpha }$, $% \,p_{\beta ,\alpha }:{\bf R}^{3\mid \alpha \mid }\times {\bf \Omega }_\beta \rightarrow \left[ 0,\infty \right) $ and $u_{\beta ,\alpha }:{\bf R}^{3\mid \alpha \mid }\times {\bf \Omega }_\beta \rightarrow {\bf R}^{3\mid \beta \mid }$ are continuous given (collision law) functions. Moreover for some $% \gamma \in {\cal M}$, \begin{equation} \label{fg}f_\gamma (t,{\bf w}):=\prod_{k\in {\cal N}_\gamma }\prod_{i=1}^{\gamma _k}f_k(t,{\bf w}_{k,i}), \end{equation} where ${\cal N}_\gamma :=\left\{ k\left| 1\leq k\leq N,\right. \gamma _k\geq 1\right\} $ and ${\bf w}=(({\bf w}_{k,i})_{i\in \left\{ {1,...,\gamma _k}% \right\} })_{k\in {\cal N}_\gamma }\in {\bf R}^{3\mid \gamma \mid }$ -the velocity space of channel $\gamma $, where ${\bf w}_{k,i}\in {\bf R}^3$, for all $i,k$. The following properties of $p_{\beta ,\alpha }$, $r_{\beta ,\alpha }$ are assumed: A$_1$) Let $m_k>0$ be the mass of a particle belonging to species $k=1,...,N$% .\ If for some $\alpha ,\beta \in {\cal M}$,% $$ \sum_{k=1}^N\alpha _km_k\neq \sum\limits_{k=1}^N\beta _km_k\;, $$ then $p_{\beta ,\alpha }\equiv r_{\beta ,\alpha }\equiv 0$ (microscopic mass conservation). A$_2$) For each pair $\alpha ,\beta \in {\cal M}$, $\varphi \in C_b($${\bf R}% ^{3\left| \alpha \right| })\ $and $f\in C_c({\bf R}^{3\mid \beta \mid })$% $$ \int_{{\bf R}^{3\mid \alpha \mid }\times {\bf \Omega }_\beta }\varphi ({\bf w% })\cdot p_{\beta ,\alpha }({\bf w},{\bf n})\cdot f(u_{\beta ,\alpha }({\bf w}% ,{\bf n}))\,d{\bf w\,}d{\bf n} $$ \begin{equation} \label{si}=\int_{{\bf R}^{3\mid \beta \mid }\times {\bf \Omega }_\alpha }\varphi (u_{\alpha ,\beta }({\bf u,n}))\cdot r_{\alpha ,\beta }({\bf u},% {\bf n})\cdot f({\bf u})\,d{\bf u\,}d{\bf n.} \end{equation} Here $C_b($${\bf R}^m)$ denotes the set of bounded elements of $C($${\bf R}% ^m,{\bf R})$, while $C_c($${\bf R}^m)$ is the space of functions of $C($$% {\bf R}^m,{\bf R})$ with compact support in ${\bf R}^m$. Property A$_2$) contains the total microscopic momentum and energy conservations during collisions, and separability of the mass center motion, in the absence of external fields (for more details, the reader is referred to APPENDIX and Ref.\thinspace 8). In this paper, we also introduce the following hypothesis: A$_3$) $\exists $ $K>0$, constant, such that for each pair $\alpha ,\beta \in {\cal M}$% \begin{equation} \label{bo}\int_{\Omega _\beta }r_{\beta ,\alpha }({\bf w},{\bf n})\,d{\bf n}% 0$, $k=1,...,N$. However, each $f_k^jd{\bf v}$, $j=1,...,\left[ \left[ \frac T{\Delta t}\right] \right] $, $k=1,...,N$, results from iteration as a sum of {\it % weighted} Dirac measures, so that it .has not the form (\ref{pri}). In addition, because of the nonlinearity of the initial problem, each step of the iteration produces a power-like growing number of terms in the sums of point measures expressing $f_k^{\;j}d{\bf v}$. In computations, the numerical effort would also be power-like increasing, so that the algorithm could not be effective at this level. To obtain a convenient approximation of $f_k^{\;j}d{\bf v}$ of the form (\ref{pri}), with non-increasing number of terms, we introduce a homogenization procedure followed by probabilistic selection.\ Then the convergence of the numerical scheme is proved in probabilistic terms.\medskip\ \begin{center} {\bf 3. EXISTENCE THEORY AND TIME DISCRETIZATION} \end{center} Let $X$ be the product of $N$ copies of $L^1({\bf R}^3;d{\bf v)}${\bf -}{\it % real} equipped with the norm% $$ \left\| g\right\| _X:=\sum_{k=1}^Nm_k\left\| g_k\right\| _{L^1}, $$ where $g=(g_1,...,g_N)$ and $g_k\in L^1({\bf R}^3;d{\bf v)}$, $k=1,...,N$. We give a meaning to (\ref{unu}) in $X$. Let $g\in C_c({\bf R}^3)^N$. By A$% _3 $), for each $k=1,...,N$, formula (\ref{e2'}) is well defined (in the sense of distributions) and \begin{equation} \label{pro}S_k(g)({\bf v})=R_k(g)({\bf v})\cdot g_k({\bf v}) \end{equation} for all ${\bf v}\in {\bf R}^3$, where% $$ R_k(g)({\bf v}) $$ \begin{equation} \label{r}:=\sum_{\alpha ,\beta \in {\cal M}}\alpha _k\int_{{\bf R}^{3\mid \alpha \mid -3}\times {\bf \Omega }_\beta }d{\bf w}_{(k)}{\bf \,}d{\bf n\,}% \left[ r_{\beta ,\alpha }({\bf w},{\bf n})\cdot g_{\alpha ,k}({\bf w}% )\right] {\bf _{w_{k,\alpha _k}=v},} \end{equation} $d{\bf w}_{(k)}$ being the Euclidean volume element induced by $d{\bf w}$ on the manifold $\left\{ {\bf w}\in {\bf R}^{3\mid \alpha \mid }\,\mid {\bf w}% _{k,\alpha _k}={\bf v}\right\} \simeq {\bf R}^{3\mid \alpha \mid -3}$ and \begin{equation} \label{fk}g_{\alpha ,k}({\bf w}):=\prod_{i=1}^{\alpha _k-1}g_k({\bf w}% _{k,i})\prod_{k\neq s\in {\cal N}_\alpha ,}\prod_{i=1}^{\alpha _s}g_s({\bf w}% _{s,i}). \end{equation} Also using (\ref{si}) (with $\varphi =1$, $f=g_\beta $) and Fubini's theorem, it follows that for each $k=1,...,N$, formula (\ref{e2}) is well defined, in the sense of distributions, and we can write% $$ P_k(g)({\bf v}) $$ \begin{equation} \label{ec}:=\sum_{\alpha ,\beta \in {\cal M}}\alpha _k{\bf \,}\int_{{\bf R}% ^{3\mid \alpha \mid -3}\times {\bf \Omega }_\beta }\,d{\bf w}_{(k)}{\bf \,}d% {\bf n}\left[ p_{\beta ,\alpha }({\bf w},{\bf n})\cdot g_\beta (u_{\beta ,\alpha }({\bf w},{\bf n}))\right] _{{\bf w}_{k,\alpha _k}={\bf v}}\; \end{equation} for all ${\bf v}\in {\bf R}^3$ (in (\ref{r}) and (\ref{ec}), the terms corresponding to $\alpha _k=0$ vanish by definition). Obviously by A$_3$), the maps, $g\rightarrow S_k(g)$ and $g\rightarrow P_k(g),$ are continuous from $g\in C_0({\bf R}^3)^N\subset X$ to $L^1({\bf R}^3;d{\bf v)}$. respectively). Moreover, $g\rightarrow R_k(g)$ is continuous from $C_c({\bf R% }^3)^N\subset X$ to $C_b({\bf R}^3)$ (equipped with the usual sup norm). Since $C_c({\bf R}^3)^N$ is dense in $X$, the maps $g\rightarrow S_k(g)$, $% g\rightarrow P_k(g)$ and $g\rightarrow R_k(g)$ have continuous extensions to $X$. These extensions (also denoted $S_k$, $P_k$ and $R_k$, respectively) are needed to formulate our Cauchy problem for equation (\ref{unu}) in $X$. With this notations, remark that (\ref{pro}) can be extended to all $g\in X$ if the pointwise equality is replaced by a.e. equality. Define $P,S:X\rightarrow X$ by $P(g)=(P_1(g),...,P_N(g))$ and $% S(g)=(S_1(g),...$, $S_N(g))$, for all $g\in X,$ respectively. The Cauchy problem for equation (\ref{unu}) in $X$ can be now written as $% \frac {}{}$% \begin{equation} \label{e4}d_tf(t)=P(f(t))-S(f(t)),\quad f(0)=f_0. \end{equation} We can prove the following existence and uniqueness result. In (\ref{e4}), let the components of $f_0=(f_{1,0},...,f_{N,0})$ be positive, i.e. $f_{k,0}>0$ for all $k=1,...,N$. \begin{theorem} \label{e-u} For each $T>0$, equation (\ref{e4}), has a unique solution in X on $\left[ 0,T\right] $, with $f_k(t)>0$, $t\in \left[ 0,T\right] $, for all $k=1,...,N.$ Moreover, \begin{equation} \label{cons}\sum_{k=1}^Nm_k\int_{{\bf R}^3}f_k(t,{\bf v})d{\bf v}% =\sum_{k=1}^Nm_k\int_{{\bf R}^3}f_{k,0}({\bf v})d{\bf v},\;\;\forall t\in \left[ 0,T\right] . \end{equation} \end{theorem} {\bf Proof: }The result is a{\bf \ }consequence of the Banach fixed point theorem applied to (\ref{e4}) written in convenient form. Consider the cone $C_T^{+}:=\left\{ f\in C(0,T;X)\mid \,f_k(t)\geq 0,\forall t\in \left[ 0,T\right] \right\} $ with the norm% $$ \left\| f\right\| =\sup _{0\leq t\leq T}\left\| f(t)\right\| _X\,. $$ First, observe that if $f\in C_T^{+}$,\thinspace then $R_k(f)\in C(0,T;C_b(% {\bf R}^3))$; $k\in \left\{ 1,...,N\right\} $. Then the Riemann integral $% \int_s^tR_k(f(\tau ))\,\,\,d\tau $ is well defined in $C_b({\bf R}^3)$ for all $0\leq s\leq t\leq T$; $k\in \left\{ 1,...,N\right\} $. Further since $% r_{\beta ,\alpha }$, $p_{\beta ,\alpha }\geq 0$, for all $\alpha ,\beta \in {\cal M}$, it follows that $R_k(f),P_k(f)\geq 0$. Consequently (\ref{e4}) can be written, in equivalent form, in $C_T^{+}$, as \begin{equation} \label{fix}f=I(f), \end{equation} where the map $\left[ 0,T\right] \ni $$t\rightarrow I(f)(t)\in X$ is given by the components of $I(f)(t)$% $$ I_k(f)(t)=exp\left[ -\int_0^tR_k(f(\tau ))\,\,\,d\tau \right] \cdot f_{k,0} $$ \begin{equation} \label{int}+\int_0^t\,exp\left[ -\int_s^tR_k(f(\tau ))\,\,d\tau \right] \cdot P_k(f(s))\,ds \end{equation} for all $t\in \left[ 0,T\right] $; $k=1,...,N$ (the integration with respect to $ds$ being in the sense of Riemann in $L^1({\bf R}^3;d{\bf v)}$) It follows that $I_k(f)(t)\geq 0$ for all $t\in \left[ 0,T\right] $, $% k=1,...,N$. For some $R>\left\| f_0\right\| $, let ${\cal B}(R):=\left\{ f\in C_T^{+}\mid \left\| f\right\| \leq R,\quad f\left( 0\right) =f_0\right\} $. After a few computations invoking (\ref{r}), (\ref{ec}) and A% $_3$), it follows that there are some numbers $k_1(R)$, $k_2(R)\geq 0$ such that% $$ \left\| I(f)\right\| \leq \left\| f_0\right\| +T\cdot k_1(R), $$ and% $$ \left\| I(f)-I(h)\right\| \leq T\cdot k_2(R)\cdot \left\| f-h\right\| \, $$ for all $f,h\in {\cal B}(R)$. Choosing $T$ small enough, the map $I$ becomes a strict contraction on ${\cal B}(R)$ and has a unique fixed point in ${\cal % B}(R)$. Consequently (\ref{e4}) has a unique positive solution on $\left[ 0,T\right] $. On the other side, the positivity of $f_k,$ implies $$ \left\| f(t)\right\| _X=\sum_{k=1}^Nm_k\int_{{\bf R}^3}f_k(t,{\bf v})\,d{\bf % v,\ }0{\bf \leq }t{\bf \leq }T. $$ Then by (\ref{e4}) we get \begin{equation} \label{zero}d_t\left\| f(t)\right\| _X=\sum_{k=1}^Nm_k\int_{{\bf R}^3}\left[ P_k(f)-S_k(f)\right] d{\bf v}=0, \end{equation} the last equality resulting by means of (\ref{r}), (\ref{ec}) and (\ref{si}) (with $\varphi \equiv 1$), and applying condition A$_1$). This proves (\ref {cons}). Moreover,% $$ \left\| f\right\| =\sup _{0\leq t\leq T}\left\| f(t)\right\| _X=\left\| f_0\right\| _X. $$ Then, by continuation, and uniqueness ,the local solution $f(t)$ can be made time-global, so the proof is complete. $\Box $ {\bf Remark: }For collision laws more general than (\ref{bo}), the existence of solutions has been obtained$^{7,\,8}$ in spaces of functions decaying in position and velocity at infinity. In the rest of this section we introduce a discretized version of (\ref{e4}% ). Let $\Delta t0,\;\;\;\;j=1,...,\left[ \left[ \frac T{\Delta t}\right] \right] ;\;\;k=1,...,N \end{equation} where $f^{\;j}=(f_1^{\;j},...,f_N^{\;j})$ and $f_k^{\;j}=f_k^{\;j}({\bf v})$% . It follows that the mass conservation is always fulfilled. Indeed, by induction and the same argument as for (\ref{zero}) we get \begin{equation} \label{mass}\sum_{k=1}^Nm_k\int_{{\bf R}^3}f_k^{\;j}({\bf v})\,d{\bf v=}% \sum_{k=1}^Nm_k\int_{{\bf R}^3}f_{k,0}({\bf v})\,d{\bf v} \end{equation} for all $j=1,...,\left[ \left[ \frac T{\Delta t}\right] \right] $. The resulting discretized scheme (\ref{it}) may destroy the positivity of the functions $f_k^{\;j}$ for $j\geq 1$. However, it can be shown that for $% \Delta t$ small enough, $f_k^{\;j}$ is positive and close, in some sense, to $f_k^{\;}$ given by Theorem \ref{e-u}$.$ \begin{proposition} \label{iitt} If $\Delta t$ is sufficiently small in the iteration scheme (% \ref{it}), then, for each $k\in \left\{ 1,...,N\right\} $, the function $% f_k^{\;j}=f_k^{\;j}({\bf v})$ is positive. Moreover, for $\Delta t$ small enough, there is some number $k(\left\| f_0\right\| _X)>0$ such that \begin{equation} \label{ineg}\left\| f(t)-f^{\,j}\right\| _X\leq k(\left\| f_0\right\| _X)\cdot \Delta t, \end{equation} for all $t\in \left( (j-1)\Delta t,\,j\Delta t\right] $, $j=1,...,\left[ \left[ \frac T{\Delta t}\right] \right] $. \end{proposition} {\bf Proof: }First we prove (\ref{ineg}). To this end, from (\ref{e4}) and (% \ref{it}), we get $$ \left\| f(j\Delta t)-f^{\,j}\right\| _X\leq \left\| f((j-1)\Delta t)-f^{\,j-1}\right\| _X+\int_{(j-1)\cdot \Delta t}^{j\cdot \Delta t}\left\| P(f(s))-P(f^{\,j-1})\right\| _Xds $$ \begin{equation} \label{ult}+\int_{(j-1)\cdot \Delta t}^{j\cdot \Delta t}\left\| S(f(s))-S(f^{\,j-1})\right\| _Xds \end{equation} for all $j=1,2,...,\left[ \left[ \frac T{\Delta t}\right] \right] $. Denote $% a_j:=\left\| f(j\Delta t)-f^{\,j}\right\| _X\;$. Using the explicit forms of $S$ and $P$ and the conservation relations (\ref{cons}), (\ref{mass}), we find that there is some number $c>0$, depending on $\left\| f_0\right\| _X$ such that $a_j0$ depends only on $\left\| f_0\right\| _X$. Suppose that $t\in \left( (j-1)\Delta t,\,j\Delta t\right] $. Again, from (\ref{e4}), the explicit forms of $S$ and $P$ and (\ref{cons}), we get% $$ \left\| f(t)-f((j-1)\Delta t)\right\| _X $$ \begin{equation} \label{a1}\leq \int_{(j-1)\Delta t}^{j\Delta t}(\left\| P(f(s))\right\| _X\,+\left\| S(f(s))\right\| _X)ds\leq c_0\Delta t, \end{equation} where $c_0$ depends only on $\left\| f_0\right\| _X$. Now (\ref{ineg}) is immediate from (\ref{aj}) and (\ref{a1}). It remains to prove the positivity of $f_k^j$. To this end, we write (\ref {it}) more conveniently. Define ${\cal U}:=\left\{ \gamma =(\gamma _k)_{1\leq n\leq N}\,\;\left| \gamma _k\in \left\{ 0,1,...,NM\right\} ,\;\left| \gamma \right| \geq 2\right. \right\} .$ Using the multinomial formula we have for $\xi =(\xi _1,...,\xi _N)\in {\bf R}^N$, $k=1,...,N$, \begin{equation} \label{mul}\sum_{p=2}^{NM}(\xi _1+...+\xi _N)^{p-1}=\sum_{p=2}^{NM}p^{-1}\partial _{\xi _k}(\xi _1+...+\xi _N)^p=\sum_{\alpha \in {\cal U}}c_\alpha \alpha _kT_{\alpha ,k}(\xi )\,, \end{equation} where% $$ c_\alpha :=\frac{(\left| \alpha \right| -1)!}{\alpha _1!\cdot ...\cdot \alpha _N!}\,\,, $$ and% $$ T_{\alpha ,k}(\xi ):=\left\{ \begin{array}{l} \xi _1^{\alpha _1}\cdot ...\cdot \xi _{k-1}^{\alpha _{k-1}}\cdot \xi _k^{\alpha _k-1}\cdot \xi _{k+1}^{\alpha _{k+1}}\cdot ...\cdot \xi _N^{\alpha _N}\;\;\;if\;\alpha _k\geq 1 \\ 0\;\;if\;\alpha _k=0 \end{array} \right. . $$ If \begin{equation} \label{csi}\xi _1+...+\xi _N=1, \end{equation} then, by (\ref{mul}) we get% $$ MN-1 $$ \begin{equation} \label{ss}=\left[ \frac 1{(M+1)^N-N-1}\sum_{\alpha ,\beta \in {\cal M}% }\alpha _kc_\alpha T_{\alpha ,k}(\xi )\,+\sum_{\alpha \in {\cal U\setminus M}% }\alpha _kc_\alpha T_{\alpha ,k}(\xi )\,\right] \end{equation} For each $k=1,...,N$, let $\xi _k=\mu _kI_k$, where $\mu _k=m_k\left( \sum_{k=1}^Nm_k\right) ^{-1}$ and% $$ I_k=\int_{{\bf R}^3}f_k^{\;j}({\bf v})\;\,d{\bf v.} $$ It follows that (\ref{csi}) is satisfied, due to (\ref{mass}). Consequently by (\ref{ss}) $$ 1=\sum_{\alpha ,\beta \in {\cal M}}\alpha _k\cdot \Gamma _{\alpha ,k}I_1^{\alpha _1}\cdot ...\cdot I_{k-1}^{\alpha _{k-1}}\cdot I_k^{\alpha _k-1}\cdot I_{k+1}^{\alpha _{k+1}}\cdot ...I_N^{\alpha _N} $$ \begin{equation} \label{un}+\sum_{\alpha \in {\cal U\setminus M}}\Lambda _{\alpha ,k}.I_1^{\alpha _1}\cdot ...\cdot I_{k-1}^{\alpha _{k-1}}\cdot I_k^{\alpha _k-1}\cdot I_{k+1}^{\alpha _{k+1}}\cdot ...I_N^{\alpha _N}, \end{equation} where% $$ \Lambda _{\alpha ,k}:=\frac{\alpha _kc_\alpha \mu _1^{\alpha _1}\cdot ...\mu _{k-1}^{\alpha _{k-1}}\cdot \mu _k^{\alpha _k-1}\cdot \mu _{k+1}^{\alpha _{k+1}}\cdot ...\cdot \mu _N^{\alpha _N}}{MN-1} $$ $$ \Gamma _{\alpha ,k}:=\frac{c_\alpha \mu _1^{\alpha _1}\cdot ...\mu _{k-1}^{\alpha _{k-1}}\cdot \mu _k^{\alpha _k-1}\cdot \mu _{k+1}^{\alpha _{k+1}}\cdot ...\cdot \mu _N^{\alpha _N}}{(MN-1)\left[ (M+1)^N-N-1\right] } $$ Multiplying the first term of the right side of (\ref{it}) by (\ref{un}) and using (\ref{r}), equation (\ref{it}) becomes \begin{equation} \label{it2}f_k^{\;j}=Q_k(f^{j-1})+L_k(f^{j-1})+\Delta t\cdot P_k(f^{\;j-1})\,, \end{equation} where $P_k$ is given by (\ref{ec}) and $$ Q_k(f^j)({\bf v}):=\sum\limits_{\alpha ,\beta \in {\cal M}}\alpha _k\int\nolimits_{{\bf R}^{3\mid \alpha \mid -3}}d{\bf w}_{(k)} $$ \begin{equation} \label{d}\times \int_{{\bf \Omega }_\beta }d{\bf n\,}\left[ \left( \frac{% \Gamma _{\alpha ,k}}{V({\bf \Omega }_\beta )}-\Delta t\cdot r_{\beta ,\alpha }({\bf w},{\bf n})\,\right) \ f_\alpha ^j({\bf w})\right] _{{\bf w}_{\alpha ,k}={\bf v}}, \end{equation} \begin{equation} \label{L}L_k(f^j)({\bf v}):=\sum_{\alpha \in {\cal U\setminus M}{\em \ }% }\Lambda _{\alpha ,k}\int_{{\bf R}^{3\mid \alpha \mid -3}}d{\bf w}_{(k)}{\bf % \ }\left[ f_\alpha ^{\,j}({\bf w})\right] _{{\bf w}_{\alpha ,k}={\bf v}} \end{equation} Here $V({\bf \Omega }_\beta )$ is the volume of the unit sphere ${\bf \Omega }_\beta $. Let $d:=\inf _{\alpha ,\beta }\frac{\Gamma _{\alpha ,k}}{V({\bf % \Omega }_\beta )}$. By condition $A_3$), we can choose $\Delta t$ such that \begin{equation} \label{delta}\Delta t\cdot \int_{{\bf \Omega }_\beta }r_{\beta ,\alpha }(% {\bf w},{\bf n})\,d{\bf n}From Proposition 2 we have the following immediate result for the discrepancy $D(f_k(j\Delta t)d{\bf v,\,}f_k^{\,j}d{\bf v)}$ \begin{corollary} \label{cor}% $$ \max _{1\leq k\leq N}\max _{1\leq j\leq \left[ \left[ \frac T{\Delta t}\right] \right] }D(f_k^{\,j}d{\bf v},\;f_k(j\Delta t)d{\bf v)\rightarrow }% 0\;\;\;\;as\;\Delta t\rightarrow 0. $$ \medskip\ \end{corollary} \begin{center} {\bf 4. THE PROBABILISTIC APPROXIMATION SCHEME} \end{center} The central result of this section extends, in some sense, the probabilistic selection methods of Babovsky and Illner$^{1,\,2}$ (see e.g. Lemma 2 of Ref.\thinspace 1). For each $n\in {\bf N^{*}}$, let ${\cal I}_n:=\left\{ 1,2,...,m_n\right\} ^{p_n}$ be an index set (with $m_n$, $p_n\in {\bf N^{*}}$ and $m_n\rightarrow \infty $ as $n\rightarrow \infty $). Consider some given set ${\bf X}\subset {\bf R}^m$ and a given sequence $(F_n)_{n\in {\bf N^{*}}% } $ of functions $F_n:{\bf X}\times {\cal I}_n\rightarrow {\bf R}$. Define $% S_n:{\bf X}\rightarrow {\bf R}$ by \begin{equation} \label{sum}S_n(x)=\frac 1{m_n^{p_n}}\sum_{{\bf j}\in {\cal I}_n}F_n(x,{\bf j}% ) \end{equation} Suppose that there is some function $F:{\bf X\rightarrow R}$ such that, for each $x\in {\bf X}$, \begin{equation} \label{lim}F(x)=\lim _{n\rightarrow \infty }S_n(x). \end{equation} In general, for a given $n$, the sum $S_n$ contains $m_n^{p_n}$ terms. Roughly speaking, our problem is to conveniently diminish the numbers of terms in $S_n$, by random selection of the terms in (\ref{sum}) and ''renormalize'' the resulting sum such that the convergence to $F(x)$ be kept, in some sense. First we construct two suitable families of random multi-indices. Let $\Omega =\left[ 0,1\right) ^\infty $ (in the countable sense) be endowed with the usual product Borel $\sigma -$algebra $\beta _\Omega $ and $P$ the usual product probability induced on $\Omega $ by the uniform distribution of $\left[ 0,1\right) $. For some given $m_n$ and for some $l\in {\bf N^{*}}$, define $c_{n,l}:\Omega \rightarrow \left\{ 1,2,...,m_n\right\} $ by $c_{n,l}(\omega ):=\left[ \left[ \omega _l\cdot m_n\right] \right] +1$, where $\omega _l$ is the $% l^{th}$ component of $\omega =(\omega _1,\omega _2,....)\in \Omega $. Obviously,% $$ P(\left\{ c_{n,l}(\omega )=j\right\} =\frac 1{m_n}\;\;\;\;\forall \;j=1,...,m_n. $$ Consequently, for some $q\in {\bf N^{*}}$ and distinct $i_1,i_2,...,i_q\in {\bf N^{*}}$, the random variables $c_{l_1,n},...,c_{l_{q,n}}$ are i.i.d.. Then for each $n$ $\in {\bf N^{*}}$, we introduce the following families of row-wise independent, random multi-indices. a) Let $p_n\geq 2$ for all $n\in {\bf N^{*}}$. For each $n$, consider the family $({\bf J}_{n,i})_{1\leq i\leq m_n}$ of random vectors ${\bf J}% _{n,i}:\Omega \rightarrow \left\{ 1,2,...,m_n\right\} ^{p_n}$ defined on components by \begin{equation} \label{prima}{\bf J}_{n,i}(\omega ):=(i,\;\;c_{n,\,(i-1)p_n+1}(\omega ),\;\;c_{n,\,(i-1)p_n+2}(\omega ),...,\;\;c_{n,\,ip_n-1}(\omega )) \end{equation} for all $\omega =(\omega _1,\omega _2,....)\in \Omega $; $i\in \left\{ 1,...,m_n\right\} $. b) Let $p_n\geq 1$ for all $n\in {\bf N^{*}}$ and $\left\{ k_n\right\} _{n\in {\bf N^{*}}}$ a given sequence of positive integers. For each $n$, consider the family $({\bf C}_{n,i})_{1\leq i\leq k_n}$ of random vectors $% {\bf C}_{n,i}:\Omega \rightarrow \left\{ 1,2,...,m_n\right\} ^{p_n}$ defined (on components) by \begin{equation} \label{fam}{\bf C}_{n,i}(\omega ):=(c_{n,\,(i-1)p_n+1}(\omega ),\;\;c_{n,\,(i-1)p_n+2}(\omega ),...,\;\;c_{n,\,ip_n}(\omega )) \end{equation} for all $\omega =(\omega _1,\omega _2,....)\in \Omega $; $i\in \left\{ 1,...,k_n\right\} $. The row-wise independence of $({\bf J}_{n,i})_{1\leq i\leq m_n}$ and $({\bf C% }_{n,i})_{1\leq i\leq k_n}$, respectively, is immediate observing that $% ip_n+j=i^{\prime }p_n+j^{\prime }$ iff $i=i^{\prime }$ and $j=j^{\prime }$, for all $i,i^{\prime }\in N^{*}$ and $j,j^{\prime }=1,2,...,p_n$. Suppose that one of the following conditions is fulfilled: 1) ${\bf X}$ is at most countable. 2) ${\bf X}$ is the whole ${\bf R}^m$, the function $F$ is continuous and each $F_n(\cdot ,{\bf j})$ is increasing with respect to the order of ${\bf R% }^m$. Define \begin{equation} \label{cnx}a_n(x):=\max _{{\bf j}\in {\cal I}_n}\left| F_n(x,{\bf j})\right| . \end{equation} \begin{proposition} \label{sel-gen} a) Let $p_n\geq 2$ for all $n\in {\bf N^{*}}$ and $% \sum\limits_{n=1}^\infty \frac{a_n(x)^4}{m_n^2}<\infty $ for all $\,\,x\in {\bf X}$. Then for each $x\in {\bf X}$, with probability one, \begin{equation} \label{con1}\lim _{n\rightarrow \infty }\frac 1{m_n}\sum_{i=1}^{m_n}F_n(x,\cdot )\circ {\bf J}_{n,i}=F(x) \end{equation} b) Let $k_n\leq m_n$ for all $n\in {\bf N}^{*}$, with $k_n\rightarrow \infty $ as $n\rightarrow \infty $ and $\sum\limits_{n=1}^\infty \frac{a_n(x)^4}{% k_n^2}<\infty $ for all $\,x\in {\bf X}$. Then for each $x\in {\bf X}$, with probability one,% $$ \lim _{n\rightarrow \infty }\frac 1{k_n}\sum_{i=1}^{k_n}F_n(x,\cdot )\circ {\bf C}_{n,i}=F(x). $$ \end{proposition} This is a consequence of the following strong-law of large numbers for arrays of row-wise independent random variables. Let $(\Omega ,\beta ,P)$ be a probability space. For some real random variable $X$, by $\langle X\rangle $ we denote its mean with respect to $P$. \begin{proposition} \label{slln} Let $((X_{n,i})_{1\leq i\leq q_n})_{n\in {\bf N^{*}}}$ be a given array of row-wise independent random variables on $\Omega $, with zero mean. Denoting $A_n:=\sup _{1\leq i\leq q_n}\langle X_{n,i}^4\rangle $, suppose that $\sum\limits_{n=1}^\infty \frac{A_n}{q_n^2}<\infty $. Then, with probability one, $\frac 1{q_n}\sum\limits_{i=1}^{q_n}X_{n,i}\rightarrow 0$, as $n\rightarrow \infty $. \end{proposition} {\bf Proof: }Our result is a simple generalization of the strong law of large numbers for i.i.d. random variables with bounded fourth momentum (see e.g. Theorem IV.\S 3-1 in Ref 16, p.363). Therefore we shall only sketch the proof. According to a version of the Borel-Cantelli Lemma, it is sufficient to show that for each $\epsilon >0$,% $$ \sum_{n=1}^\infty P(\left| \frac 1{q_n}\sum_{i=1}^{q_n}X_{n,i}\right| >\epsilon )<\infty \,. $$ To this end, by Chebyshev's inequality, we obtain% $$ P(\left| \sum_{i=1}^{q_n}X_{n,i}\right| >\epsilon \cdot q_n)\leq \frac 1{\epsilon ^4q_n^4}\langle \left| \sum_{i=1}^{q_n}X_{n,i}\right| ^4\rangle $$ Expanding the fourth power, we invoke the independence of $X_{n,i}$ and use the fact that $\langle X_{n,i}\rangle =0$. Then a simple computation shows that $\forall \epsilon >0$,% $$ 0\leq \sum\limits_{n=1}^\infty P(\frac 1{q_n}\left| \sum_{i=1}^{q_n}X_{n,i}\right| >\epsilon )\leq \frac 3{\epsilon ^4}\sum_{n=1}^\infty \frac{A_n}{q_n^2}<\infty . $$ \ This concludes the proof. $\Box $ {\bf Proof }of{\bf \ }Proposition \ref{sel-gen} a){\bf \ }First remark that it is sufficient to consider the case in which all functions $F_n$ are positive. $a_1$) Case ${\bf X}$ countable. Let $x\in {\bf X\ }$be fixed. For each $i\in \left\{ 1,...,m_n\right\} $, define $Y_{n,i}:=F_n(x,\cdot )\circ {\bf J}_{n,i}$. Clearly, the independence of $({\bf J}_{n,i})_{1\leq i\leq m_n}$ implies that $% (Y_{n,i})_{1\leq i\leq n}$ is a row-wise independent family of random variables.\ Let ${\bf j}=(j_1,...,j_{p_n})\in {\cal I}_n$. From the definition of ${\bf J}_{n,i}$ we get easily $$ P\left\{ {\bf J}_{n,i}(\omega )={\bf j}\right\} =\left\{ \begin{array}{l} m_n^{1-p}\quad if\;i=j_1 \\ 0\quad if\;i\not =j_1 \end{array} \right. \;\;\;\;\forall \;{\bf j}\in {\cal I}_n. $$ Consequently,% $$ \langle Y_{n,i}\rangle =\frac 1{m_n^{p-1}}\sum_{j_2,...,j_{p_n}=1}^{m_n}F_n(x,(i,\;j_2,...,j_{p_n})), $$ so that \begin{equation} \label{media}\frac 1{m_n}\sum_{i=1}^{m_n}\langle Y_{n,i}\rangle =\frac 1{m_n^p}\sum_{{\bf j}\in {\em I}_n{\em .}}F_n(x,{\bf j})=S_n(x). \end{equation} Define $X_{n,i}:=Y_{n,i}-\langle Y_{n,i}\rangle $. Then the family $(($ $% X_{n,i})_{1\leq i\leq m_n})_{n\in {\bf N^{*}}}$ satisfies the conditions of Proposition \ref{slln}, with $A_n\leq (2a_n(x))^4$. Therefore, $x$ being fixed, by (\ref{media}) and (\ref{sum}) we get the limit (\ref{con1}). For each $x\in {\bf X}$, let $\Omega _x$ be the subset of $\Omega $ where the limit (\ref{con1}) holds. Define $\Omega {\bf _X}:=\bigcap_{x\in {\bf X}% }\Omega _x.$ Since ${\bf X\ }$is countable, we have $P(\Omega _{{\bf X}})=1$% , so that the argument is complete. $a_2$) Case ${\bf X=R}^m$. Observe that the argument of $a_1$) is valid on the countable set ${\bf Q}^m$ of the vectors of ${\bf R}^m$ with rational components. Further, remark that for any $x\in {\bf R}^m{\bf \setminus Q}^m$ and $\epsilon >0$, by the continuity of $F$ and the monotonicity of $F_n$, we can find two elements $x^{-},x^{+}\in {\bf Q}^m$, with $x^{-}\leq x\leq x^{+}$ such that% $$ F(x^{+})-\frac 1{m_n}\sum_{i=1}^{m_n}F_n(x^{+},\cdot )\circ {\bf J}% _{n,i}(\omega )-\epsilon $$ $$ \leq F(x)-\frac 1{m_n}\sum_{i=1}^{m_n}F_n(x,\cdot )\circ {\bf J}% _{n,i}(\omega ) $$ $$ \leq F(x^{-})-\frac 1{m_n}\sum_{i=1}^{m_n}F_n(x^{-},\cdot )\circ {\bf J}% _{n,i}(\omega )+\epsilon $$ for all $\omega $. No we approximate $x$ by two sequences $\left\{ x_p^{+}\right\} _{p\in {\bf N}}$, $\left\{ x_p^{-}\right\} _{p\in {\bf N}% }\subset {\bf Q}^m$, with $x_p^{-}\leq x\leq x_p^{+}$, for $p$ sufficiently large. Then to conclude the proof of $a_2$) we refer to the result $a_1)$. b) With $Y_{n,i}:=F_n(x;\cdot )\circ {\bf C}_{n,i}$, the previous argument works similarly. $\Box $ In the following, if $\mu $ is some finite measure on ${\bf R}^m$, then set $% \left| \mu \right| :=\mu ({\bf R}^m)$. The index set ${\cal I}${\em $_n$} being defined as before, let $((\mu _{n,% {\bf j}})_{{\bf j}\in {\cal I}_n})_{n\in {\bf N^{*}}}$ be a bounded family of positive measures on ${\bf R}^m$ (i.e. $\exists $ some constant $a>0$ such that $\left| \mu _{n,{\bf j}}\right| \leq a$ for all ${\bf j}\in {\cal I% }_n$, $n\in {\bf N^{*}}$). Suppose that there is a positive measure $\mu $ on ${\bf R}^m$, absolutely continuous with respect to the Lebesgue measure on ${\bf R}^m$ such that \begin{equation} \label{limas}\frac 1{m_n^{p_n}}\sum_{{\bf j}\in {\cal I}_n}\mu _{n,{\bf j}}% \stackrel{w}{\rightarrow }\mu ,\,\;as\;n\rightarrow \infty . \end{equation} Let $(\Omega ,\beta _\Omega ,P)$, $({\bf J}_{n,i})_{1\leq i\leq m_n}$ and $(% {\bf C}_{n,i})_{1\leq i\leq k_n}$be as in Proposition \ref{sel-gen}. \begin{theorem} \label{sel}a) Let $p_n\geq 2$ for all $n\in {\bf N^{*}}$ and define $\mu _{n,i}(\omega ):=\mu _{n,{\bf j}}\mid _{_{{\bf j}={\bf J}_{n,i}(\omega )}}$ for all $\omega \in \Omega $. If $\sum\limits_{n=1}^\infty m_n^{-2}<\infty $% , then for $P-$almost all $\omega $, \begin{equation} \label{s1}\sigma _{1,n}(\omega ):=\frac 1{m_n}\sum\limits_{i=1}^{m_n}\mu _{n,i}(\omega )\stackrel{w}{\rightarrow }\mu ,\,\;as\;n\rightarrow \infty . \end{equation} b) Let $p_n\geq 1$ for all $n\in {\bf N^{*}}$ and define $\mu _{n,i}(\omega ):=\mu _{n,{\bf j}}\mid _{_{{\bf j}={\bf C}_{n,i}(\omega )}}$ for all $% \omega \in \Omega $. If $k_n\leq m_n$ for all $n$ and $\sum\limits_{n=1}^% \infty k_n^{-2}<\infty $, then for $P-$almost all $\omega $, \begin{equation} \label{s2}\sigma _{2,n}(\omega ):=\frac 1{k_n}\sum\limits_{i=1}^{k_n}\mu _{n,i}(\omega )\stackrel{w}{\rightarrow }\mu ,\,\;as\;n\rightarrow \infty . \end{equation} \end{theorem} {\bf Proof: }Define% $$ F_n(x,{\bf j):=}\int\nolimits_{y\leq x}d\mu _{n,{\bf j}}(y),\quad \forall x\in {\bf R}^m $$ and% $$ F(x):=\int\nolimits_{y\leq x}d\mu (y),\quad \forall x\in {\bf R}^m. $$ Then it is sufficient to observe that $F$ and $F_n(x,{\bf j)\ }$satisfy the conditions of Proposition \ref{sel-gen}, (with $a_n(x)=a$) and the family $% \left\{ y\in {\bf R}^m\,\left| \,y\leq x\right. \right\} _{x\in {\bf R}^m}$ is determining (see Ref.\thinspace 16) for the weak convergence of the measures $\mu _{n,{\bf j}}$. $\Box $ It can be easily seen that Babovsky Lemma (see Lemma 2 of Ref.\thinspace 1) is a particular case of Theorem \ref{sel}a) with ${\cal I}_n${\em $:=\left\{ 1,...,n\right\} ^2$}, and $\mu _{n,{\bf j}}$ given by a product of two point measures. Theorem \ref{sel} will be the basic point of the probabilistic part of our numerical scheme for the solutions of (\ref{unu}) in the next section. As we have mentioned in Section 1, our purpose is to approximate the solutions of (\ref{unu}) by sums of Dirac measures of the form (\ref{pri}). However, except the case of (\ref{unu}) modelling the one component gas with purely elastic collisions, a certain step of the numerical scheme destroys the homogeneity of the sums of Dirac measures. This difficulty will be surmounted by using the following simple homogenization result. Let $((a_{n,i})_{1\leq i\leq n})_{n\in {\bf N^{*}}}$ be a family of positive numbers and $((\mu _{n,i})_{1\leq i\leq n})_{n\in {\bf N^{*}}}$ a family of positive finite measures on ${\bf R}^m$. For some $c>1$ fixed, set $% s_n:=\sum\limits_{i=1}^n\left[ \left[ n^c\cdot a_{n,i}+1\right] \right] $. Also define the function $\Lambda $$:\left\{ 1,...,s_n\right\} \rightarrow \left\{ 1,...,n\right\} $ as \begin{equation} \label{li}\Lambda (i):=\min _{1\leq j\leq n}\left\{ j\left| i\leq \left( \sum\limits_{k=1}^j\left[ \left[ n^c\cdot a_{n,k}+1\right] \right] \right) \right. \right\} . \end{equation} Suppose that there is a positive sequence $\left\{ b_n\right\} _{n\in {\bf % N^{*}}}$ such that $\left| \mu _{n,i}\right| =b_n$, for all $1\leq i\leq n$, $n\in {\bf N^{*}}$. \begin{lemma} \label{omo} If for some finite measure $\mu $, \begin{equation} \label{lunu}\sum\limits_{i=1}^na_{n,i}\cdot \mu _{n,i}\stackrel{w}{% \rightarrow }\mu \quad as\;\,n\rightarrow \infty , \end{equation} then \begin{equation} \label{ldoi}\frac{\left| \mu \right| }{s_nb_n}\sum\limits_{i=1}^{s_n}\cdot \mu _{n,\Lambda (i)}\stackrel{w}{\rightarrow }\mu \quad as\;n\rightarrow \infty . \end{equation} \end{lemma} {\bf Proof:} Let $\varphi \in C_b({\bf R}^m)$. From the definition of $% \Lambda $, we obtain $$ \left| \frac 1{n^c}\sum\limits_{i=1}^{s_n}\cdot \mu _{n,\Lambda (i)}(\varphi )-\sum\limits_{j=1}^na_{n,j}\mu _{n,j}(\varphi )\right| $$ $$ \leq \left| \sum\limits_{i=1}^n\frac{\left[ \left[ n^c\cdot a_{n,i}\right] \right] +1}{n^c}-a_{n,i}\right| b_n\left\| \varphi \right\| \leq \frac 1{n^{c-1}}b_n\left\| \varphi \right\| . $$ Now it is sufficient to observe that $\sum_{i=1}^na_{n,i}b_n\rightarrow \left| \mu \right| $ as $n\rightarrow \infty $, implying $\frac{s_nb_n}{n^c}% \rightarrow \left| \mu \right| $ as $n\rightarrow \infty $.$\Box $\medskip\ \begin{center} {\bf 5.\ THE MAIN RESULT} \end{center} >From the considerations of the previous sections, it turns out that it would be convenient to approximate $f_k^{\,j}({\bf v})d{\bf v}$, by suitable weakly convergent sums of discrete measures. So we need a weak form of equation (\ref{it2}). Denote% $$ (\varphi ,h)=\int_{{\bf R}^3}\varphi ({\bf v})h({\bf v})d{\bf v\,,\ }\varphi \in C_b({\bf R}^3),\;h\in L^1({\bf R}^3;d{\bf v).} $$ Then (\ref{it2}) leads to \begin{equation} \label{weak}\left( \varphi ,f_k^{\;j}\right) =(\varphi ,Q_k(f\,^{j-1}))+(\varphi ,L_k(f\,^{j-1}))+\Delta t\cdot (\varphi ,P_k(f^{\;j-1})) \end{equation} for all $\varphi \in C_b({\bf R}^3)$; $j=1,...,\left[ \left[ \frac T{\Delta t}\right] \right] $, $k=1,...,N$. In equation (\ref{weak})% $$ (\varphi ,Q_k(f\,^j)):=\sum_{\alpha ,\beta \in {\cal M}}\alpha _k\int_{{\bf R% }^{3\mid \alpha \mid }\times {\bf \Omega }_\beta }{\bf (}\varphi \circ i_{\alpha ,k})({\bf w}) $$ \begin{equation} \label{dw}\times \left( \frac{\Gamma _{\alpha ,k}}{V({\bf \Omega }_\beta )}% -\Delta t\cdot r_{\beta ,\alpha }({\bf w},{\bf n})\,\right) \ f\,_\alpha ^j(% {\bf w})\,d{\bf w\,}d{\bf n,} \end{equation} \begin{equation} \label{lw}L_k(f\,^j)({\bf v}):=\sum_{\alpha \in {\cal U\setminus M}{\em \ }% }\Lambda _{\alpha ,k}\int_{{\bf R}^{3\mid \alpha \mid }}{\bf (}\varphi \circ i_{\alpha ,k})({\bf w})f_\alpha ^{\,j}({\bf w})\,d{\bf w} \end{equation} for all $\varphi \in C_b({\bf R}^3)$, where $i_{\gamma ,k}:$ ${\bf R}^{3\mid \gamma \mid }\rightarrow $ ${\bf R}^3$ is the projection map $i_{\gamma ,k}(% {\bf w})={\bf w}_{k,\gamma _k}$, ${\bf \gamma }\in {\cal M}$; $k=1,...,N$. Moreover, using (\ref{e2}) and hypothesis $A_2$) we get, for any $\varphi \in C_b({\bf R}^3)$,% $$ (\varphi ,P_k(f\,^j)) $$ \begin{equation} \label{pw}=\sum_{\alpha ,\beta \in {\cal M}}\beta _k\int_{{\bf R}^{3\mid \alpha \mid }\times {\bf \Omega }_\beta }\varphi \circ i_{\beta ,k}(u_{\beta ,\alpha }({\bf w},{\bf n}))\,r_{\beta ,\alpha }({\bf w},{\bf n})\ f_\alpha ^j({\bf w})\,d{\bf w\,}d{\bf n.} \end{equation} Consequently we can write (\ref{weak}) as an equation for measures. For $\mu _1$,..., $\mu _N$ positive measures on ${\bf R}^3$ and ${\bf \gamma }\in {\cal M}$, define the measure $\mu _\gamma $ on ${\bf R}^{3\mid \gamma \mid } $ by \begin{equation} \label{mpro}d\mu _\gamma ({\bf w})=\otimes _{k\in {\cal N}_\gamma }\otimes _{i=1}^{\gamma _k}d\mu _k({\bf w}_{k,i}). \end{equation} Set $d\mu _k^j\left( {\bf v}\right) :=f_k^j\left( {\bf v}\right) d{\bf v}$, $% j=0,1,...,\left[ \left[ \frac T{\Delta t}\right] \right] $, $k=1,...,N$. Using (\ref{dw}), (\ref{lw}) and (\ref{pw}), it is not very difficult to see that we can find some $L\in {\bf N}^{*}$, $\left\{ \alpha (l)\right\} _{1\leq l\leq L}\subset $${\cal M}$ and $\left\{ q(l)\right\} _{1\leq l\leq L}\subset {\bf N}^{*}$, such that (\ref{weak}) takes the form \begin{equation} \label{ecm}\int_{{\bf R}^3}\varphi ({\bf v})\,d\mu _k^j({\bf v}% )=\sum_{l=1}^L\int_{{\bf R}^{3\left| \alpha (l)\right| +q(l)}}R_{k,l}({\bf z)% }(\varphi \circ h_{k,l})({\bf z})\,d(\mu _{\alpha (l)}^{j-1}\otimes \pi _l)(% {\bf z)}\,, \end{equation} $\forall \varphi \in C_b({\bf R}^3)$, $j=1,...,\left[ \left[ \frac T{\Delta t}\right] \right] $. Here $R_{k,l}\in C({\bf R}^{3\left| \alpha (l)\right| +q(l)},\left[ 0,\infty \right) )$ and $h_{k,l}\in C({\bf R}^{3\left| \alpha (l)\right| +q(l)},{\bf R}^3)$ $k=1,...,N$, $l=1,...,L$, are given functions. Also each $\pi $$_l$, $l=1,...,L$, is a finite nonnegative measure, absolute continuous with respect to the Lebesgue measure on ${\bf R}^{q(l)}$. Equation (\ref{ecm}) is the basic object for the numerical approximation scheme solving equation (\ref{unu}). Therefore, for a given timestep $\Delta t$, we write the weak form (\ref{ecm}) of the time discretized version (\ref {it}) of equation (\ref{e4}). Then we proceed as follows. Let $(\Omega ,\beta _\Omega ,P)$ be as in Theorem \ref{sel}. a) For each $l=1,...,L$, we approximate $\pi $$_l$ by a convenient HSPM of the form (\ref{pri}), containing $n$-terms, $\pi _{l,n}\stackrel{w}{% \rightarrow }\pi _l$ as $n\rightarrow \infty $ (this can be done, e.g., by means of low discrepancy, well distributed sequences$^{4,\,5}$). b) The initialization of the scheme is done by giving $n$-terms HSPM approximations $\nu _{k,n}^0$ of the initial data, $\mu _k^0$; $k=1,...,N$. c) The $n$-terms HSPM approximations $\nu _{k,n}^1$of $\mu _k^1$, $k=1,...,N$% , resulting from the scheme, can be obtained as follows: Step 1 (First selection ). We introduce $\nu _{k,n}^0$, $k=1,...,N$, in the product formula (\ref{mpro}% ) for $\gamma =\alpha (l)$; $l=1,...,L$. Consequently, for each $l=1,...,L$, we obtain a sequence of finite measures $\nu _{\alpha (l),n}^0\stackrel{w}{% \rightarrow }\mu _{\alpha (l)}^0$ as $n\rightarrow \infty $, implying $\nu _{\alpha (l),n}^0\otimes \pi _{l,n}\stackrel{w}{\rightarrow }$ $\mu _{\alpha (l)}^0\otimes \pi _l$ as $n\rightarrow \infty $. Obviously, each $\nu _{\alpha (l),n}^0\otimes \pi _{l,n}$ is a sum of point measures of the form (% \ref{ecm}), containing $n^{\left| \alpha (l)\right| +1}$ terms. We apply the selection algorithm cf. Theorem \ref{sel} a) (with $m_n=n$ and $p_n=\left| \alpha (l)\right| +1$) to construct $n$-terms HSPM approximations for all $% \nu _{\alpha (l),n}^0\otimes \pi _{l,n}$. Thus, by Theorem \ref{sel} a), for each $l=1,...,L$, there exists some set $\Omega _l\subset \Omega $, with $% P(\Omega _l)=1$, such that from $\nu _{\alpha (l),n}^0\otimes \pi _{l,n}$, one can extract a $n$-terms HSPM approximation (of the form (\ref{s1})) $% \sigma $$_{1,l,n}(\omega _l)\stackrel{w}{\rightarrow }$$\mu _{\alpha (l)}^0\otimes \pi _l$ as $n\rightarrow \infty $, for all $\omega ^l\in \Omega _l$. Step 2 (Homogenization). In the right side of (\ref{ecm}), written for $j=1$, we replace each $\mu _{\alpha (l)}^0\otimes \pi _l$ by the corresponding $\sigma $$_{1,l,n}$. Consequently the right side of (\ref{ecm}) defines a measure $M_k$ on ${\bf R% }^3$, \begin{equation} \label{mk}dM_{k,n}({\bf v})=\frac 1n\sum_{l=1}^L\sum_{i=1}^na_l(R_{k,l}\circ {\bf z}_{l,n,i})(\omega ^l)\delta ({\bf v}-h_{k,l}({\bf z}_{l,n,i}(\omega ^l)), \end{equation} for some ${\bf z}_{l,n,i}(\omega ^l)\in {\bf R}^{3\left| \alpha (l)\right| +q(l)}$ and $a_l>0$; $l=1,...,L$, $i=1,...,n$. By Step 1, it follows that $% M_{k,n}\stackrel{w}{\rightarrow }\mu _k^1$ as $n\rightarrow \infty $, for all $\omega ^1\in \Omega _1,\,.\omega ^2\in \Omega _2,..,\,\omega ^L\in \Omega _L$; $k=1,...,N$. Now it can be easily seen that (\ref{mk}) can be written as a sum (containing, at most $L\cdot n$ terms) of point measures of the form (\ref{lunu}). Then by Lemma \ref{omo}, a measure $M_{k,n}^{\prime }$% , of the form (\ref{ldoi}), containing $s_{k,n}=s_{k,n}(\omega ^1,...,\omega ^L)=\sum_{l=1}^L\sum_{i=1}^n\left[ \left[ n^{c-1}a_l(R_{k,l}\circ {\bf z}% _{l,n,i})(\omega ^L)+1\right] \right] $ terms, can be associated to each $% M_{k,n}$. Step 3 (Second selection). We fix, for the moment, some $\omega ^1\in \Omega _1,...,\omega ^L\in \Omega _L$, so that $M_{k,n}\stackrel{w}{\rightarrow }\mu _k^1$ (and implicitly $% M_{k,n}^{\prime }\stackrel{w}{\rightarrow }\mu _k^1$) as $n\rightarrow \infty $; $k=1,...,N$. We apply the selection algorithm cf. Theorem \ref{sel} b) (with $p_n=1$, $m_n=s_{k,n}$, and $k_n=n$). Then for each $k$, there is some $\Omega _{L+k}\subset \Omega $ with $P(\Omega _{L+k})=1$, $k=1,...,N$, such that from $M_{k,n}^{\prime }$, we obtain a $n$-terms HSPM approximation (of the form (\ref{s2})) $\sigma _{2,k,n}(\omega ^{L+k};\ \omega ^1,...,\omega ^L)\stackrel{w}{\rightarrow }\mu _k^1$ as $n\rightarrow \infty $, for all $\omega ^{L+k}\in \Omega _{L+k}$. Set $\tilde \nu _{k,n}^1(\omega ^1,...,\omega ^{L+k}):=\sigma _{2,k,n}(\omega ^{L+k};\omega ^1,...,\omega ^L) $. Therefore for each $\mu _k^1$ in (\ref{ecm}),we obtain a corresponding $n$-terms HSPM approximation $\tilde \nu _{k,n}^1\stackrel{w}{% \rightarrow }\mu _k^1$ as $n\rightarrow \infty $, for all $\omega $$^1\in \Omega _1,...,\omega ^{L+k}\in \Omega _{L+k}$; $k=1,...,N$. e) The procedure can be repeated, with the entering data $\tilde \nu _{k,n}^1 $, to obtain similar approximations $\tilde \nu _{k,n}^2(\omega ^1,...,\omega ^{2L+N+k})$ of $\mu _k^2$, $k=1,...,N$. f) Repeating this procedure over and over, after $j$ timesteps, we provide the $n$-terms HSPM approximations $\tilde \nu _{k,n}^j(\omega ^1,...,\omega ^{jL+(j-1)N)+k})\stackrel{w}{\rightarrow }\mu _k^j$ for all $\omega $$^1\in \Omega _1,\,\omega ^2\in \Omega _2,...,\,\omega ^{jL+(j-1)N)+k}\in \Omega _{jL+(j-1)N)+k}$, $j=1,...,\left[ \left[ \frac T{\Delta t}\right] \right] $, $k=1,...,N$, where $\Omega _l\subset \Omega $ with $P(\Omega _l)=1$, $% l=1,...,\left[ \left[ \frac T{\Delta t}\right] \right] (L+N)$. Now observe that we can find a family $\left\{ Q_l\right\} _{l\in {\bf N^{*}}% }$ of measurable maps $Q_l:\Omega \rightarrow \Omega $, with $% P(Q_l^{-1}(A))=1$, for all $A\subset \Omega $ with $P(A)=1$. For instance, we can consider $U,V:\Omega \rightarrow \Omega $. given by $$ U(\omega )=U(\omega _1,\omega _2,...,\omega _{2n-1},\omega _{2n},...):=(\omega _1,\omega _3,...\omega _{2n-1},\omega _{2n+1},...) $$ $$ V(\omega )=V(\omega _1,\omega _2,...,\omega _{2n-1},\omega _{2n},...):=(\omega _2,\omega _4,...,\omega _{2n},\omega _{2n+2},...) $$ respectively, for all $\omega =$$(\omega _1,\omega _2,...,\omega _{2n-1},\omega _{2n},...)\in \Omega $. Then it is sufficient to put $Q_1=U$ and $Q_l:=U\circ V^{l-1}$, $l=2,3,...$. Let $$ \Omega _{\Delta t}:=\bigcap_{l=1}^{\left[ \left[ \frac T{\Delta t}\right] \right] (L+N)}Q_l^{-1}(\Omega _l). $$ Since $P(Q_l^{-1}(\Omega _l))=1$ for all $l=1,...,\left[ \left[ \frac T{\Delta t}\right] \right] (L+N)$, clearly $P(\Omega _{\Delta t})=1$. Define $\nu $$_{k,n}^j(\omega ):=\tilde \nu _{k,n}^j(Q_1(\omega ),...,Q_{jL+(j-1)N+k}(\omega ))$ for all $\omega \in $$\Omega $; $% j=1,...,\left[ \left[ \frac T{\Delta t}\right] \right] $, $k=1,...,N$. Then it follows that $\nu _{k,n}^j(\omega )\stackrel{w}{\rightarrow }\mu _k^j$ as $n\rightarrow \infty $, for all $\omega \in $$\Omega _{\Delta t}$; $% j=1,...,\left[ \left[ \frac T{\Delta t}\right] \right] $, $k=1,...,N$. In particular, if $D(\cdot ,\cdot )$ is the discrepancy introduced in Section 1, then \begin{equation} \label{dae}\lim _{n\rightarrow \infty }\max _{1\leq k\leq N}\max _{1\leq j\leq \left[ \left[ \frac T{\Delta t}\right] \right] }D\left( \nu _{k,,n}^j(\omega ),\mu _k^j\right) =0 \end{equation} for almost all $\omega \in \Omega $. All these and Corollary \ref{cor} lead to our main result: Let $\,f(t)=(f_1(t),...,f_N(t))$ be the solution of equation (\ref{e4}), provided by Theorem \ref{e-u}. Consider some family $\left\{ \Delta t_p\right\} _{p\in {\bf N}}$ of discretization timesteps as in Section 3. For each $\Delta t_p$, consider the solutions $\mu _{k,p}^j$ of (\ref{ecm}), corresponding to the initial data $d\mu _k^0({\bf v})=f_k(0,{\bf v})d{\bf v}$% , $j=1,...,\left[ \left[ \frac T{\Delta t_p}\right] \right] $, $k=1,...,N$. Suppose that for each $\mu _{k,p}^j$, there is given a $n$-terms HSPM approximation $\nu _{k,p,n}^j$ by the aforementioned scheme. \begin{theorem} \label{main}For each sequence of timesteps $\Delta t_p\rightarrow 0$ as $% p\rightarrow \infty $, there is a sequence of positive integers $% n(p)\rightarrow \infty $ as $p\rightarrow \infty $, such that $$ \lim _{p\rightarrow \infty }\max _{1\leq j\leq \left[ \left[ \frac T{\Delta t_p}\right] \right] }D\left( \nu _{k,p,n(p)}^j(\omega ),\,f\left( j\cdot \Delta t_p,{\bf v}\right) d{\bf v}\right) =0 $$ for almost all $\omega \in $$\Omega $. \end{theorem} {\bf Proof: }Denote% $$ d_{p,n}(\omega ):=\max _{1\leq k\leq N}\max _{1\leq j\leq \left[ \left[ \frac T{\Delta t_p}\right] \right] }D\left( \nu _{k,p,n}^j(\omega ),\mu _{k,p}^j\right) . $$ Consider some positive sequence $\epsilon _p\downarrow 0$ as $p\rightarrow \infty $. By (\ref{dae}), for each $p$, $\lim _{n\rightarrow \infty }P(d_{p,n}>\epsilon _p)=0$. Then for each $p$, we can choose $n=n(p)$ such that $P(d_{p,n(p)}>\epsilon _p)\leq \frac 1{p^2}$. Consequently $$ \sum_{p=1}^\infty P(d_{p,n(p)}>\epsilon _p)<\infty , $$ so that $\lim _{n\rightarrow \infty }d_{p,n(p)}(\omega )=0$ for almost all $% \omega \in \Omega $. Now the proof can be easily concluded by using Corollary \ref{cor}. $\Box $ This theorem can be considered as a (space homogeneous) reactive correspondent to the main result in the Babovsky-Illner simulation scheme for the classical Boltzmann equation (Theorem 7.1 of Ref.\thinspace 2).% \medskip\ \begin{center} {\bf 6.\ CONCLUSIONS} \end{center} For a given number of terms in the sum (\ref{lunu}) the approximation errors introduced by Lemma \ref{omo} can be controlled arbitrarily well by suitable choice of $c$, for a numerical effort of order at most $n\log n$ (the last being introduced by the calculation of $\Lambda (i)$). It can be seen that, in the worst case, the numerical effort of the method is at the most of order $n\log n.$ However this situation can be improved, depending on the form of the collision laws. Thus for equation (\ref{unu}) modelling a gas where all the inelastic collisions are velocity independent, the numerical effort is of order $n$. In particular, this is true for the gas with purely elastic (possibly multiple) collisions. Moreover, in the case of the classical Boltzmann equation for the simple one-component gas, the scheme of the previous section reduces, in some sense, to that one introduced in Ref.\thinspace 1. It should be remarked that in Ref.\thinspace 1 the simulation is applied to a convenient functional form of the discretized Boltzmann equation, which seems not obtainable in the reactive case, because of the presence of reactions thresholds. The qualities and the drawbacks of the method as well as the errors, introduced by the probabilistic part of the argument of Section 4, are similar to those in Ref.\thinspace 1,\thinspace 2. Thus, good numerical results are expected if in Theorem \ref{main}, $n(p)$ increases much faster than $\left[ \left[ \frac T{\Delta t_p}\right] \right] $ (see also Ref.\thinspace 2). A space dependent version of Section 3 can be obtained in the frame of the results of Ref.\thinspace 7,\thinspace 8. Then the result of Theorem \ref {main} appears as generalizable to the space dependent case by using the cell homogenization procedure introduced by Babovsky and Illner$^2$ for the case of the classical Boltzmann equation. The method presented here has been tested$^{14}$ on the exact solutions provided by Krook and Wu$^{12}$ on the multicomponent Boltzmann equation. The discrepancies show values close to those presented in the case in which the methods of Ref 2 are used to simulate the solutions of the classical Boltzmann equation.\medskip\ \begin{center} {\bf ACKNOWLEDGEMENTS} \end{center} Part of this work was written at the Institute of Mathematics Jussieu-CNRS, Paris VII University. One of us (C.P.G) is particularly indebted to Prof. Anne Boutet de Monvel for the warm hospitality during his visit at the Laboratory of Mathematical Physics and Geometry. \begin{center} {\bf APPENDIX} \end{center} According to Ref.\thinspace 9-11, kinetic equations with common features can be obtained for different models of real gas with inelastic collisions, using the classical Boltzmann argument, starting from the molecular chaos hypothesis. An abstract form of these equations can be obtained$^{7,\,8}$ by the same Boltzmann-type argument, if one takes into account the following observation. In certain cases, it is possible to treat different internal states of a real gas particle with internal structure as distinct structureless point particles with given mass and unique internal energy state$^{7,\,8}$. Consequently we can think of the gas of particles with internal structure as a gas mixture of different mass-points, with unique internal energy, and re-write the original kinetic equations for the real gas, by re-labeling the original distribution functions (each original distribution function, describing a succession of internal states of a particle with internal structure, is replaced by a sequence of distribution functions associated to each internal state). The space homogeneous version of the equations of Ref.\thinspace 7,\thinspace 8 can be formally written as% $$ \partial _tf_k(t,{\bf v})=P_k(f)(t,{\bf v})-S_k(f)(t,{\bf v})\quad k=1,...,N, $$ where% $$ S_k(f)(t,{\bf v}):=\sum_{\alpha ,\beta \in {\cal M}}\alpha _k\int_{{\bf R}% ^{3\mid \beta \mid }\times {\bf R}^{3\mid \alpha \mid }}\ \ f_\alpha (t,{\bf % w})K_{\alpha ,\beta }({\bf w},{\bf u}) $$ $$ \times \ \delta ({\bf w}_{k,\alpha _k}-{\bf v})\ \delta (V_\beta ({\bf u}% )-V_\alpha ({\bf w}))\ \delta (W_\beta ({\bf u})-W_\alpha ({\bf w}))\ d{\bf % u\,}d{\bf w},\;\eqno (A.1) $$ $$ P_k(f)(t,{\bf v}):=\sum_{\alpha ,\beta \in {\cal M}}\ \alpha _k\ \int_{{\bf R% }^{3\mid \beta \mid }\times {\bf R}^{3\mid \alpha \mid }}\,f_\beta (t,{\bf u}% )\,K_{\beta ,\alpha }({\bf u},{\bf w})\ $$ $$ \times \ \delta ({\bf w}_{k,\alpha _k}-{\bf v})\ \delta (V_\beta ({\bf u}% )-V_\alpha ({\bf w}))\ \delta (W_\beta ({\bf u})-W_\alpha ({\bf w}))\ d{\bf % u\,}d{\bf w},\;\eqno (A.1^{\prime }) $$ Here, besides the notations of Section 1, $V_\gamma ({\bf w})$ and $W_\gamma ({\bf w})$ are the classical mass center velocity and the total energy, respectively, for the particles in channel $\gamma $, i.e.,% $$ V_\gamma ({\bf w}):=(\sum_{n=1}^N\gamma _nm_n)^{-1}\sum_{n\in {\cal N}% _\gamma }\sum_{i=1}^{\gamma _n}m_n{\bf w}_{n,i},\;\eqno (A.2) $$ $$ W_\gamma ({\bf w}):=\sum_{n\in {\cal N}_\gamma }\sum_{i=1}^{\gamma _n}(2^{-1}m_n{\bf w}_{n,i}^2+E_n).\ \eqno (A.2^{\prime }) $$ where $E_n\in {\bf R}$ the internal energy of the particle belonging to species $n=1,...,N$. Further, $(K_{\alpha ,\beta })_{(\alpha ,\beta )\in {\cal M}\times {\cal M}}$ is the family of given (transition) functions $% K_{\alpha ,\beta }:{\bf R}^{3\mid \alpha \mid }\times {\bf R}^{3\mid \beta \mid }\rightarrow \left[ 0,\infty \right) $, with the following properties. k$_1$) If for some $\alpha $, $\beta \in {\cal M}$,$\;\sum_{n=1}^N\alpha _nm_n\not =\sum_{n=1}^N\beta _nm_n$, then $K_{\alpha ,\beta }\equiv 0$ (microscopic mass conservation). k$_2$) For each ${\bf w}$, ${\bf u}$ and $n\in {\cal N}_\alpha $ fixed, $% K_{\alpha ,\beta }({\bf w},{\bf u})$ is invariant at the interchange of components ${\bf w}_{n,1},...,{\bf w}_{n,\alpha _n}$ of ${\bf w}$ and similarly with respect to the interchange of the components of ${\bf u}$ (indistinguishability of identical particles). k$_3$) For each $a\in {\bf R}^3$, define the map ${\bf w\rightarrow }T(a)% {\bf w}$ by setting $(T(a){\bf w})_{n,i}:={\bf w}_{n,i}+a$, for all $n,i$; then $K_{\alpha ,\beta }({\bf w},{\bf u})\equiv K_{\alpha ,\beta }(T(a){\bf w% },{\bf u})\equiv K_{\alpha ,\beta }({\bf w},T(a){\bf u})$, $\forall $ $a\in {\bf R}^3$, $({\bf w,u)}\in {\bf R}^{3\left| \alpha \right| }\times {\bf R}% ^{3\mid \beta \mid }$ and $\alpha $, $\beta \in {\cal M}$ (absence of external fields)$.$ In (A.1) and (A.1'), the Dirac $\delta $-functions exhibit the conservation of the total momentum and energy (provided that condition k$_1$) is fulfilled). We introduce the following correspondent of the boundedness assumption of Ref.\thinspace 1,\thinspace 2 on the collision law for the classical Boltzmann equation. For some $\gamma \in {\cal M}$, let% $$ W_{r,\gamma }({\bf w}):=W_\gamma ({\bf w})-2^{-1}(\sum_{n=1}^N\gamma _nm_n)V_\gamma ({\bf w})^2-\sum_{n=1}^N\gamma _nE_n,\;\eqno (A.3) $$ be the corresponding mass center kinetic energy. Obviously, $W_{r,\gamma }(% {\bf w})\ge 0$.\ {\bf ASSUMPTION A.}{\it -a) For each pair }$\alpha $, $\beta \in {\cal M}$,% {\it \ the transition law} $K_{\alpha ,\beta }$ {\it is continuous on the set }$\{({\bf w},{\bf u})\in {\bf R}^{3\mid \alpha \mid }\times {\bf R}% ^{3\mid \beta \mid }\mid W_{r,\alpha }({\bf w})>0,W_{r,\beta }({\bf u})>0\}$. b) $\exists $ $C>0$,$\;${\it constant,} {\it such that} $$ K_{\alpha ,\beta }({\bf w},{\bf u})\le C\,\left( W_{r,\alpha }({\bf w}% )^{(3\mid \alpha \mid -5)/2}+W_{r,\beta }({\bf u})^{(3\mid \beta \mid -5)/2}\right) ^{-1}.\;\eqno (A.4) $$ We give a precise meaning to (A.1) and (A.1'). Let $\delta _\epsilon :{\bf R}% \rightarrow \left[ 0,\infty \right) $, $\epsilon >0$, be an even mollifier with $supp\,\delta _\epsilon =[-\epsilon ,\epsilon ]$ (i.e. $\delta _\epsilon (t)=:\epsilon ^{-1}J(t/\epsilon )$, for some even, positive function $J\in C_c({\bf R})$, with $supp\,J=[-1,1]$ and $\left\| J\right\| _{L^1}=1$). Set $\delta _\epsilon ^3(y):=\delta _\epsilon (y_1)\cdot \,\delta _\epsilon (y_2)\cdot \,\delta _\epsilon (y_3)$, where $% y=(y_1,y_2,y_3)\in {\bf R}^3$. Define% $$ S_{k\epsilon \eta }(g)({\bf v}):=\sum_{\alpha ,\beta \in {\cal M}}\alpha _k\int_{{\bf R}^{3\mid \beta \mid }\times {\bf R}^{3\mid \alpha \mid -3}}\ d% {\bf u\,}d{\bf w}\ \ \left[ g_\alpha ({\bf w})\right. $$ $$ \left. \times \,K_{\alpha ,\beta }({\bf w},{\bf u)\,}\delta _\epsilon ^3(V_\beta ({\bf u})-V_\alpha ({\bf w}))\,\delta _\eta (W_\beta ({\bf u}% ))-W_\alpha ({\bf w}))\right] _{{\bf w}_{k,\alpha _k}={\bf v}},\eqno (A.5) $$ $$ P_{k\epsilon \eta }(g)({\bf v}):=\sum_{\alpha ,\beta \in {\cal M}}\alpha _k\int_{{\bf R}^{3\mid \beta \mid }\times {\bf R}^{3\mid \alpha \mid -3}}\ d% {\bf u\,}d{\bf w}\ \left[ g_\beta ({\bf u)}\right. $$ $$ \left. \times \,K_{\beta ,\alpha }({\bf u},{\bf w})\,\delta _\epsilon ^3(V_\beta ({\bf u})-V_\alpha ({\bf w}))\,\delta _\eta (W_\beta ({\bf u}% ))-W_\alpha ({\bf w}))\right] _{{\bf w}_{k,\alpha _k}={\bf v}},\eqno % (A.5^{\prime }) $$ for all $g\in C_c({\bf R}^3)^N$; ${\bf v}\in {\bf R}^3$, $1\le k\le N$. We show that if $g\in C_c({\bf R}^3)^N$, then $S_k$ and $P_k$ can be defined as limits of (A5) and (A5'), respectively. To this end we need the following transformation. Let $n$ be a non-negative integer and $a_1,\ldots ,a_n>0$, constants. Consider a positive quadratic form $T:=T({\bf v}_1,\ldots ,{\bf v}% _n)=\sum_{i=1}^na_i{\bf v}_i^2$ on ${\bf R}^{3n}$, $r_i\in {\bf R}^3$, $1\le i\le n$. Consider the Jacobi-type transformation% $$ \qquad {\bf R}^{3n}\ni ({\bf v}_1,\ldots ,{\bf v}_n)\rightarrow (\underline{V% },\xi )\in {\bf R}^3\times {\bf R}^{3n-3},\eqno (A.6) $$ defined by% $$ \underline{V}:=(\sum_{i=1}^na_i)^{-1}\sum_{i=1}^na_i{\bf v}_i, $$ $$ \xi :=(\xi _1,\ldots ,\xi _{n-1}), $$ where% $$ \xi _i:=\mu _i^{1/2}\left[ {\bf v}_{i+1}-(\sum_{j=1}^ia_j)^{-1}% \sum_{j=1}^ia_j{\bf v}_j\right] ,\quad \;i=1,\ldots ,n-1, $$ $$ \mu _i^{-1}:=a_{i+1}^{-1}+(\sum_{j=1}^ia_j)^{-1},\qquad i=1,\ldots ,n-1. $$ By transformation (A.6), the form $T$ becomes% $$ T=T(\underline{{\sl V}},\xi )=(\sum_{i=1}^na_i)\underline{{\sl V}}^2+\xi ^2, $$ We have the following result (for a more general case see Ref.\thinspace 7,\thinspace 8). \smallskip\ {\bf Proposition A 1. }{\it If} $g\in C_c({\bf R}^3)${\it , then for each} $k=1,...,N$, {\it there exist the limits}% $$ P_k(g)({\bf v)}=\lim \limits_{\eta \rightarrow 0}\lim \limits_{\epsilon \rightarrow 0}P_{k\epsilon \eta }(g)({\bf v})\;and\;R_k(g)({\bf v})=\lim \limits_{\eta \rightarrow 0}\lim \limits_{\epsilon \rightarrow 0}R_{k\epsilon \eta }(g)({\bf v}),\eqno (A.7) $$ {\it for all} ${\bf v\in R}^3$. {\bf Proof.} Consider the form $T_\beta ({\bf u})=W_\beta ({\bf u}% )-\sum_{n=1}^N\beta _nE_n$ on ${\bf R}^{3\mid \beta \mid }$, and a Jacobi-type transformation of the form (A.6), ${\bf R}^{3\mid \beta \mid }\ni {\bf u}\rightarrow $ $(\underline{{\sl V}},\xi )\in {\bf R}^3\times {\bf R}^{3\mid \beta \mid -3}$ with $\xi :=(\xi _1,\ldots ,\xi _{\mid \beta \mid -1})$, $\,\xi _i\in {\bf R}{^3}$,${\ \;}i=1,\ldots ,\mid \beta \mid -1$% . Let $\xi $ be represented in spherical coordinates on ${\bf R}^{3\mid \beta \mid -3}$, $\xi =r{\bf n}$ , with $(r,{\bf n})\in \left[ 0,\infty \right) \times \Omega _\beta $. In (A.5) and (A.5'), choose $(\underline{V}% ,r,{\bf n})$ as new integration variables such that ${\bf u}={\bf u}( \underline{V},r,{\bf n})$. Then the limits follow by repeated application of Lebesgue's dominated convergence theorem, using the properties of $K_{\alpha ,\beta }$, $\delta _\epsilon ^3$ and $\delta _\eta $. $\Box $\ The proof of Prop. A.1 provides the limits (A.5), (A.5') in explicit form. Indeed, define% $$ H_{r,\alpha }({\bf w}):=W_{r,\alpha }({\bf w})+\sum_{n=1}^N(\alpha _n-\beta _n)E_n, $$ $$ t_{\beta ,\alpha }({\bf w})=\left\{ \begin{array}{l} H_{r,\alpha }( {\bf w})^{1/2}\;\;\;\,if\;H_{r,\alpha }({\bf w})\geq 0,\medskip \\ 0,\;\;otherwise. \end{array} \right. . $$ If $H_{r,\alpha }({\bf w})\geq 0$, then set ${\bf u}_{\beta \alpha }({\bf w},% {\bf n}):={\bf u}(\underline{V},r,{\bf n})_{\mid \underline{V}=V_\alpha (% {\bf w}),\;r=t_{\beta ,\alpha }({\bf w})}$. For the sake of simplicity, $% {\bf u}_{\beta \alpha }$ will replace the notation ${\bf u}_{\beta \alpha }(% {\bf w},{\bf n})$. Define% $$ p_{\beta \alpha }({\bf w},{\bf n}):=\left\{ \begin{array}{l} 2^{-1}\Delta _\beta \cdot t_{\beta ,\alpha }( {\bf w})^{3\mid \beta \mid -5}K_{\beta ,\alpha }({\bf u}_{\beta \alpha },% {\bf w})\quad \;if\;H_{r,\alpha }({\bf w})\geq 0\medskip\ \\ 0,\;otherwise. \end{array} \right. \; $$ $$ r_{\beta \alpha }({\bf w},{\bf n}):=\left\{ \begin{array}{l} 2^{-1}\Delta _\beta \cdot t_{\beta ,\alpha }( {\bf w})^{3\mid \beta \mid -5}K_{\alpha ,\beta }({\bf w},{\bf u}_{\beta \alpha })\quad \;if\;H_{r,\alpha }({\bf w})\geq 0\medskip \\ 0,\;otherwise. \end{array} \right. , $$ where the constant $\Delta _\beta $ is introduced by the Jacobian of the transformation ${\bf w\rightarrow }(\underline{V},r,{\bf n})$. With these notations we have the following result. \smallskip\ {\bf Proposition A.2}{\it \ If} $\varphi \in C_b($${\bf R}^{3\left| \alpha \right| })\ ${\it and} $f\in C_c({\bf R}^{3\mid \beta \mid })${\it , then}% $$ \lim \limits_{\eta \rightarrow 0}\lim \limits_{\epsilon \rightarrow 0}\int_{% {\bf R}^{3\mid \alpha \mid }\times {\bf R}^{3\mid \beta \mid }}\varphi ({\bf % w})\,K_{\beta ,\alpha }({\bf u},{\bf w}) $$ $$ \times \delta _\epsilon ^3(V_\beta ({\bf u})-V_\alpha ({\bf w}))\,\delta _\eta (W_\beta ({\bf u}))-W_\alpha ({\bf w}))\,\,f({\bf u})\,d{\bf w\,}d{\bf % u} $$ $$ =\int_{{\bf R}^{3\mid \alpha \mid }\times {\bf \Omega }_\beta }\varphi ({\bf % w})\,\,p_{\beta ,\alpha }({\bf w},{\bf n})\,f(u_{\beta ,\alpha }({\bf w},% {\bf n}))\,d{\bf w\,}d{\bf n} $$ $$ =\int_{{\bf R}^{3\mid \beta \mid }\times {\bf \Omega }_\alpha }\varphi (u_{\alpha ,\beta }({\bf u,n}))\,\,r_{\alpha ,\beta }({\bf u},{\bf n})\,\,f(% {\bf u})\,d{\bf u\,}d{\bf n,}\eqno (A.8). $$ {\bf Proof.} The first equality in (A.8) is a by-product of Prop. A.1. The second equality results by the same argument, but starting from the Jacobi-type transformation ${\bf R}^{3\mid \alpha \mid }\ni {\bf w}% \rightarrow (\underline{{\sl V}},\xi )\in {\bf R}^3\times {\bf R}^{3\mid \alpha \mid -3}$, corresponding to the form $T_\alpha ({\bf w})=W_\alpha (% {\bf w})-\sum_{n=1}^N\alpha _nE_n$ defined on ${\bf R}^{3\mid \alpha \mid }$.% $\Box $ {\bf Remark:} Assumption A) implies that (\ref{bo}) is verified by $% \,r_{\alpha ,\beta }$ and $p_{\beta ,\alpha }$. Finally, as a consequence of Prop. A.1 and Prop. A.2, the limits (A.7) can be written as (\ref{r}), and (\ref{ec}). \medskip\ \setlength{\baselineskip}{12pt} \begin{thebibliography}{99} \bibitem{b} H. Babovsky, European Journal of Mechanics B-Fluids, \underline{% 8}, 41 (1989). \bibitem{bi} H. Babovsky and R. Illner, SIAM J. Num. Anal., \underline{26}, 45 (1989). \bibitem{na} K.\ Nambu, J.\ Phys. Soc. Japan, \underline{49} 2042 (1980). \bibitem{gr} F.\ Gropengiesser, H. Neunzert, J. Struckmeier, ''Computational Methods for the Boltzmann Equation, Bericht Nr. 43, Fachbereich Mathematik,Kaiserslautern (1990). \bibitem{ni1} L. Kuipers, H. Niederreiter, ''Uniform Distribution of Sequences'', JohnWiley \& Sons (1974). \bibitem{ni2} L.\ Niederreiter, ''Random Number Generation and Quasi-Monte Carlo Methods'', SIAM (1992). \bibitem{p1} C. P. Grunfeld, C. R. Acad. Sci. Paris Tome I, \underline{316}% , 953 (1993). \bibitem{p2} C.P. Grunfeld and E. Georgescu, ''On a Class of Kinetic Equations for Reacting, Gas Mixtures'', Preprint No 604/11/93, Bielefeld Universit\"at (1993). \bibitem{wa} C.S.G. Wang Chang, E. Uhlenbeck and J. De Boer, in ''Studies in Statistical Mechanics'', J. De Boer and G.E. Uhlenbeck Eds, vol\underline{ 2}, Part C, p.241, North Holland Amsterdam (1964). \bibitem{cer} C. Cercignani, ''The Boltzmann Equation and Its Applications'', in Applied Mathematical Science, No 67, Springer (1988). \bibitem{lu} G Ludwig, M. Heil, in Advances in Applied Mechanics vol. \underline{6}, p.39, Academic Press (1960). \bibitem{kr} M.\ Krook, T. T. Wu, Phys. Rev Lett.,{\bf \ }\underline{38}, 991{\bf \ }(1977). \bibitem{ai} I.\ Aizenman, Duke Math.\ J., \underline{45}, 809 (1978). \bibitem{test} E. Georgescu, C. P. Grunfeld, D. Marinescu, Tests of a Simulation Method for the Multicomponent Boltzmann Equation, to appear in Balkan Phys.Lett.. \bibitem{mar} C. P. Grunfeld, D. Marinescu, On a Simulation Method in the Kinetic Theory of Reacting Gas Mixtures, Second International Conference in Applied and Industrial Mathematics, Oradea (1994). \bibitem{sh} A. N. Shiryayev, ''Probability'', in Graduate Texts in Mathematics, vol. \underline{95}, Springer (1984). \end{thebibliography} \end{document}