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}%
<K,\quad \forall \,\,{\bf w\in R}^{3\left| \alpha \right| }.
\end{equation}

This is the reactive correspondent of the conditions imposed on the
collision law in the numerical study of the classical Boltzmann equation$%
^{1,\,2}$. By the same considerations as in Ref.\thinspace 1,\thinspace 2,
hypothesis A$_3$) does not appear to represent a severe limitation {\it to
numerical purposes}.

Our problem is to approximate the solutions $f_k(t,{\bf v})$, $k=1,...,N$ of
(\ref{unu}) by HSPM of the form (\ref{pri}), provided that the initial data $%
f_k(0,{\bf v})$, $k=1,...,N$ are also approximated by (given) sums of point
measures of the form (\ref{pri}). Following ideas similar to those of
Ref.\thinspace 1,\thinspace 2, we are led to the following approximation
procedure.

First, using the mass conservation and choosing some fixed timestep $\Delta
t<T$, one obtains a time discretized version of (\ref{unu}). Taking $%
f_k^{\;0}=f_k(0,{\bf v})$, $k=1,...,N$, as initial data in the discretized
equation, an iteration procedure provides a set of probability densities $%
f_k^{\;j}=f_k^{\;j}({\bf v)}$, $j=1,...,\left[ \left[ \frac T{\Delta
t}\right] \right] $ (here $\left[ \left[ x\right] \right] $ denotes the
integer part of $x\in {\bf R}$). It is shown that for each $k=1,...,N$, the
function $f_k^{\;j}$ approximates the solution $f_k$ of equation (\ref{unu})
on the interval $\left( (j-1)\Delta t,\,j\Delta t\right] $, $j=1,...,\left[
\left[ \frac T{\Delta t}\right] \right] $, in some sense, with an error of
order $\Delta t$.

The iteration scheme is therefore initialized by taking suitable HSPM
approximations: $f_{k,n}^{\;0}d{\bf v:=}\frac{a_k}n\sum\limits_{i=1}^n\delta
({\bf v}-{\bf v}_{k,n})$ with $f_{k,n}^{\;0}d{\bf v}\stackrel{w}{\rightarrow
}f_k(0,{\bf v})d{\bf v}$, as\ $n\rightarrow \infty $; $a_k>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 t<T$ be a fixed timestep. Consider the following iteration
scheme%
$$
f_k^{\;j}=f_k^{\;j-1}+\Delta t\cdot \left[
P_k(f^{\;j-1})-S_k(f^{\;j-1})\right] ,
$$
\begin{equation}
\label{it}f_k^0=f_{k,0}>0,\;\;\;\;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_j<a_{j-1}\cdot (1+c\Delta t)$ for all $j=2,...,\left[ \left[
\frac T{\Delta t}\right] \right] $ and $a_1\leq c\cdot \Delta t$. Then
\begin{equation}
\label{aj}a_j\leq a_1\cdot (1+c\Delta t)^{\left[ \left[ \frac T{\Delta
t}\right] \right] }\leq c_1\Delta t,
\end{equation}
where $c_1>0$ 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}<d,\;\;\;\;\forall \,{\bf w\in R}^{3\left| \alpha
\right| },\alpha ,\beta \in {\em M.}
\end{equation}

Then, under condition (\ref{delta}), the positivity of the components of $\
f^j$, for all $j=1,...,\left[ \left[ \frac T{\Delta t}\right] \right] $,
follows by straightforward induction.$\Box $

As a consequence of Proposition 1, $f_k^{\,j}({\bf v})d{\bf v}$ and $%
f_k(j\Delta t,{\bf v})d{\bf v}$ are close, for small $\Delta t$, in the weak
sense of the measures.

>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}
ENDBODY
