key: cord-0966164-5v5egz64 authors: Lotova, G. Z.; Mikhailov, G. A. title: Numerical-Statistical and Analytical Study of Asymptotics for the Average Multiplication Particle Flow in a Random Medium date: 2021-09-19 journal: Comput DOI: 10.1134/s0965542521060075 sha: 5655332df4d90f3fbda88bbee3374d71351260c8 doc_id: 966164 cord_uid: 5v5egz64 It is well known that, under rather general conditions, the particle flux density in a multiplying medium is asymptotically exponential in time [Formula: see text] with a parameter [Formula: see text], i.e., with an exponent [Formula: see text]. If the medium is random, then [Formula: see text] is a random variable, and the time asymptotics of the average number of particles (over medium realizations) can be estimated in some approximation by averaging the exponent with respect to the distribution of [Formula: see text]. Assuming that this distribution is Gaussian, an asymptotic “superexponential” estimate for the average flux expressed by an exponential with the exponent [Formula: see text] can be obtained in this way. To verify this estimate in a numerical experiment, a procedure is developed for computing the probabilistic moments of [Formula: see text] based on randomizations of Fourier approximations of special nonlinear functionals. The derived new formula is used to study the COVID-19 pandemic. We study the time asymptotics of the average flow of scattering and multiplying particles in a random medium. For this purpose, an averaging procedure is developed for corresponding analytical and numerical statistical estimates (i.e., Monte Carlo estimates) obtained for deterministic realizations of the medium. It is well known that, under rather general conditions, the particle flux density in a multiplying medium occupying a domain is asymptotically exponential in time (see, e.g., [1] ): The time constant is the leading eigenvalue of the corresponding homogeneous stationary kinetic equation (1.1) Here, is the stationary flux density (characteristic function of Eq. (1.1)), ; is the total cross section (attenuation coefficient); , where is the scattering cross section, is the fission (multiplication) cross section ( and are the corresponding phase functions), and is the absorption cross section; is the number of particles leaving the multiplication point; is the velocity vector, where is the unit direction vector and ; and is a spatial point. To construct and study Monte Carlo algorithms, as a mathematical model of the transport process corresponding to Eq. (1.1) (see, e.g., [2] ), we use a terminating (with probability 1) homogeneous Markov chain of states that are the phase points of successive "collisions of a particle with substance elements," i.e., points at which the particle velocity changes instantaneously. This Markov chain , , …, is considered in the state space of coordinates, velocities, and time, i.e., , where is the point of the th collision, is the velocity immediately before the collision, and is the lifetime of the colliding particle. The considered chain is determined by the distribution density of the initial collision and by the transition density from the state to ; moreover, it is assumed that (1.2) i.e., the chain terminates with probability 1 and the average number of transitions is finite. Condition (1.2) is satisfied, for example, for a bounded system (see [2, 3] ). A substochastic kernel is obtained from the transport process characterization (see [2, 3] ), which also determines Eq. (1.1). The ratios , , and are the probabilities of scattering, fission, and absorption immediately after a collision at the phase point , and the distribution density of the free path from to is , where is the optical path (see [2, 3] ). It can be seen that the collision density , where is the distribution density of collisions of index , represents a Neumann series for the second-kind integral equation where is an integral operator with kernel . As was indicated, for example, in [2, 3] , despite the generalized multipliers and present in the kernel , the operator can be treated as acting from to , even more so that all functions in the problem are nonnegative. Under condition (1.2), we have and, hence, the spectral radius of the operator satisfies . Used usually in transport theory, the radiation intensity (particle flux density) is connected with the collision density by the relation . As a rule, the Monte Carlo method is used to estimate linear functionals of the form , . Weighted Monte Carlo algorithms are constructed using a Markov chain with an initial density and a transition density containing the indicated generalized multipliers, for example, a chain of collisions with other values of the parameters , , , and . The auxiliary weights are introduced using the formulas Under the "unbiasedness conditions" we have (see [2, 3] ) If, additionally, , where is an operator with kernel , and , then . The random variable is called the collision estimator for the functional . Note that the new velocity is usually modeled with probability according to the indicatrix and with probability according to . To construct a weighted modification, the index of the modeling type (i.e., of the collision type) is added to the set of coordinates of the state space (see [2, 3] ). In this section, we consider two approaches for estimating the time constant of particle multiplication, which determines the asymptotics of the function in a deterministic medium (see, e.g., [1] ). It is well known (see [1] ) that, if ( is the particle velocity) is added to the absorption cross section, then the system becomes critical, i.e., is the eigenvalue of Eq. (1.1). In the case of standard spectral properties of the resolvent operator (see [1] ), we have as , which yields the following limit form of Eq. (2.1): Thus, To estimate the derivatives of the linear functional with respect to , we model a chain of collisions with parameters and ; moreover, , where is the total time required for the particle to move to and is a weight taking into account fissions for the simplest process modification without branching (see [2] ). It follows that Since is added to the absorption cross section, the quantity modified by the substitutions and can be made smaller than unity by choosing a sufficiently large . The unbiasedness of the corresponding estimators of and the finiteness of their variances for can easily be shown by applying the vector approach; the conditions that the variances are finite are less restrictive if branching is modeled (see [3] ). The considered approach was formulated in [4] and was elaborated in [5] . Its shortcoming is that it involves expensive computations of multiple derivatives , whose variances grow rather significantly with increasing in real-world problems. Nevertheless, by additionally computing the derivatives of estimates of and using the values of a piecewise constant random function , rather accurate test estimates of and for a model spherically symmetric system were obtained in [5] ; they will be used in Section 4 below. Below, we formulated a more universal approach for estimating based on Green's function with respect to the time parameter . Let be the collision density (in ) produced by a single collision at the point , i.e., for . The functional can be represented in the form (see [5] ) where Assume that and, hence, for . In what follows, the symbol denotes the th derivative of the function with respect to . [2] Note that, in [2] , an estimate of type (i.e., for ) is used to solve nonstationary optical sensing problems with actual radiation sources. In this paper, following [5] , we use and with an auxiliary density to estimate the parameter for the exponential time asymptotics of the particle flux. It is well known, that under rather general conditions in the case of a source localized at the point , the following asymptotic relation holds as (see [1] ): where is the leading eigenvalue of Eq. (1.1). Specifically, these conditions hold for one-speed particle transport in a bounded medium with a source density rapidly decreasing with time. Relying on what was said above, we can construct an estimate of . In the case , the same result is obtained by introducing additional absorption with a coefficient . This completes the proof of Theorem 1. To simplify the presentation, we consider one-speed particle transport. It is assumed that is a random field and the ratios and and the scattering and fission phase functions are fixed. If , where is the indicator of the domain , then the functional is the total particle flux in at the given time . Assuming that the random variable is Gaussian and the passage to the limit is uniform (in ), we can estimate the asymptotics of the function as : where and are the expectation and variance of the corresponding variable. Additionally, we assume that the multipliers and in the asymptotics are weakly correlated; hence, . The assumptions made are natural, for example, for a small random perturbation of a multi- , < < +∞, ρ < , Note that formulas (3.1) and (3.2) can be used as a basis for numerical studies of particular problem instances. Determined by formula (3.1), the growth law for the average number of particles can be called "superexponential." A more general and practically convenient definition of the superexponential property can be associated with an increase in the growth coefficient with increasing . In real-world problems, ; therefore, the resulting asymptotic representation seems to approximate only in some interval . Accordingly, an additional numerical study is required. Practically important are the quantities E and D . The conceptually simplest (direct) Monte Carlo technique for their estimation is to obtain fairly accurate estimates of for a sample of sufficiently large size. However, this technique may be extremely expensive for realistic medium models and transport processes. For this reason, a randomized algorithm is constructed in [7] using the following relation ( is replaced by ): where is a preliminary statistical estimator of , which can be made deterministic in the given case. The simplest (elementary) unbiased randomized estimator of based on Lemma 1 is constructed by realizing conditionally independent particle trajectories for a chosen value of , where A similar decomposition is used to construct an estimator of (see [7] ). It can be seen that in (3.3) is a statistical estimator of . In this case, and are simultaneously estimated using the double randomization method (see [3] ). For definiteness, the result is stated as the following lemma. Lemma 2. Under the conditions of Lemma 1, (3.4) where and are the corresponding statistical estimators produced by the double randomization method. For each realization of the medium, we can construct only one particle trajectory by using (2.3) for . If the complexity of modeling the medium is substantially higher than the complexity of modeling the trajectory, then it is reasonable to construct conditionally independent trajectories for each medium realization. An optimal value of is estimated by analogy with the parameter of the splitting method (see [9] ). Note that the variance of is estimated from above with the help of linearization of the fraction using the formula (3.5) For a medium with a random piecewise constant density , an algorithm for estimating and based on the Taylor series expansions of with respect to at the point is constructed in [5] . More specifically, truncated series up to the first and second orders inclusive are used for the variance and the mean, respectively. Assuming that the random variables are independent, we have (3.6) The variance estimate can be improved using second-order terms. In the example considered below, this improvement was found insignificant. Thus, estimate (3.6) is reduced to estimating the first and second derivatives of with respect to . After replacing by an iterative resolvent approximation (see Subsection 2.1), such estimates are determined by the mixed derivatives (3.7) Quantities (3.7) can be computed as indicated in [5] by applying the Monte Carlo method with auxiliary weights corresponding to variations in the parameters and , . The results thus obtained can be independently verified using more expensive computations based on (3.3). To perform test computations, we considered one-speed particle transport in a spherically symmetric medium with a piecewise constant random density within a ball of radius with macroscopic sections , , and , where To construct a realization of the medium, the ball is divided into spherical layers of identical volume. In each layer, the random variable is drawn independently and uniformly on the interval . To construct an efficient Monte Carlo algorithm based on Lemma 2, into the formulated model, we introduced absorption with a constant nonrandom coefficient , which led to the substitution . Note that this technique is universal and can significantly improve the efficiency of the weighted method, eliminating the necessity of branching modeled trajectories (see below). By using the transport equation (see [2, 3] ), we can make the substitution In the numerical experiment, the transport process was modeled with constants , , and . The auxiliary weights were determined by the formula , where is the number of fissions preceding the given collision (see [5] ). We used the value , for which and, thus, (see [2, 3] ). The distribution density of the first collisions was specified as where is the improved diffusion approximation of the spatial characteristic function for and (see [8] ). To estimate and , we also set , where and (see [8] ) Moreover, , i.e., a functional of the particle flux was computed. The computations showed that these functional parameters of the algorithm significantly improve the convergence in (2.6) as compared with and even with the variant in which is defined by formula (4.1) and . It is easy to see that the conditions of Lemma 1 hold for the above-formulated characteristics of the computational model. The corresponding estimates of and are presented in the first two columns of Table 1 . They were obtained using distributed computations at the Siberian Supercomputer Center of the Siberian Branch of the Russian Academy of Sciences. These estimates were determined on stabilizing two significant digits with increasing and taking into account the statistical error. The final results were obtained by averaging the estimates produced at the times for and and for and . This corresponds to improved properties of the estimates as . The mean-square deviations are given as errors of these estimates. For the problem under consideration, formulas (3.6) were also implemented in [5] . In view of , they have the form The results thus obtained are given in the last two columns of Table 1 . They can be regarded as test results for and at . Figure 1 shows the estimated logarithmic derivative can be computed by applying the well-known formulas (see [10] ) In the considered problem, the random variables are positively correlated, so An analysis of the resulting estimated variances with the incomplete dependence of taken into account showed that the mean-square errors of and can be estimated using the formulas and . The analytical and regression estimates for the coefficients of the linear approximation to the function are given in Table 2 . The required values of and were determined using formulas (4.2) (see Table 1 and Table 7 in [3] ). Thus, we conclude that the analytical estimate for the exponent of the exponential time asymptotics of the particle flux is satisfactory in the interval . For , the statistical estimates of the coefficients agree better with their analytical values, which can be explained by the better Gaussian property of the distribution of for this variant of the problem, according to the theory of small perturbations. To conclude, we note that the results obtained above have a wider range of applications. They show that the average number of particles of an arbitrary nature (e.g., microorganisms) in the case of their distributed multiplication in a random medium can have a superexponential growth rate (see Section 3). This corresponds to an increasing growth coefficient of the numbers of particles at fixed equidistant times , i.e., to an increase in the ratio . However, if such an increase is observed experimentally, this may suggest a distributed nature of the multiplication "source." Specifically, according to WHO statistics (see [11] ), such behavior was exhibited by the COVID-19 pandemic developing around the world from March 9, 2020, to March 21, 2020. The corresponding number of cases (by days) is approximated within a 2% error by a formula of type (3.1), namely, Individual numbers of cases are compared in Table 3 . Neutron Transport Theory The Monte Carlo Methods in Atmospheric Optics Numerical Statistical Modeling: Monte Carlo Methods (Akademiya New Monte Carlo methods for computing critical parameter values in a particle transport equation Moments of the critical parameters of the transport of particles in a random medium Tables of Indefinite Integrals Monte Carlo algorithms for estimating time asymptotics of multiplication particle flow in a random medium Exact solutions of one-velocity equation and their application in the computation of diffusion problems (improved diffusion method)," in Investigation of Critical Parameters of Reactor Systems Improvement of multidimensional randomized Monte Carlo algorithms with 'splitting Mathematical Statistics r r r ae ρ = π σ + νσ ρ σ + νσ + .