Cooperation with: S. Orszag (Yale University, USA), T. Vesala (Helsinki University, Finland), P.K. Yeung (Georgian Institute of Technology, USA), O. Kurbanmuradov (Phys. Tech. Institute, Turkmenian Academy of Sciences, Ashkhabad), I.A. Shalimova (Institute of Computational Mathematics and Mathematical Geophysics, Russian Academy of Sciences, Novosibirsk), O. Smidts (Université Libre de Bruxelles, Belgium), G. Wei (Hong Kong Baptist University), M. Tretyakov (University of Wales, Swansea, UK)
Supported by: INTAS Project INTAS-99-1501, Nato Linkage Grant N 971664
Complex physical phenomena are conventionally governed by deterministic differential and integral equations, but the advanced measurements and modern technology require nowadays a development of stochastic research methods. The basis of these methods are stochastic models, so the main effort is often related to the creation of an adequate stochastic model which mimics, in a probabilistic sense, the desired physical phenomena. Next, these models should be efficiently modelled by computer. Mathematical analysis of the model and the implementation algorithm is an important part of research. A new type of stochastic model is developed for solving transport problems in saturated porous media which is the first stochastic Lagrangian model of Langevin type for porous media. We develop also stochastic algorithms for solving some examples of Cauchy and boundary value problems for nonlinear parabolic PDEs. The study is accomplished with an analysis of algorithms for high-dimensional integration.
1. Stochastic simulation models for transport in porous media (K.K. Sabelfeld).
It is well-known that stochastic models are well-developed for solving transport problems in turbulent flows like the transport in the atmospheric boundary layer (see, e.g., , ). This technique was extended to a wide class of flows, in particular, to rivers (), to flows in porous media (see ). In the porous media transport, only one type of stochastic models was used, namely, the random displacement method (RDM) for hydrodynamic dispersion equations. It should be stressed that RDM can be applied only if the displacement covariance tensor is known (e.g., from measurements, or numerical simulation), and cannot be applied if the functionals of interest are evaluated at times comparable with the characteristic correlation scale of the flow. In contrast, the Lagrangian stochastic models based on the tracking particles in a random velocity field extracted from the numerical solution of the flow equation (for brevity, we will call this model DSM, the direct simulation method) are free of these limitations, but the computational resources required are vast. Therefore, we suggest to construct a Langevin-type stochastic model which is an approximation to DSM, and is written in the form of a stochastic differential equation for the position and velocity. It is worth to mention that the same scheme has been carried out in the atmospheric transport problem (see our recent work , , ,  and ). The basis for the Langevin-type approach comes from the Kolmogorov similarity theory of fully developed turbulence () saying that the velocity structure tensor is a linear function in time which is universal in the inertial subrange. The linearity is the necessary condition to derive a Langevin-type equation to mimic the behavior of the real Lagrangian trajectories. Therefore, the crucial point is here to study if in the porous media , this kind of linear law can be observed. This problem is studied by the DSM. Detailed derivation of the Langevin-type model is given. Numerical simulations and comparisons with the random displacement model confirm that the new Lagrangian approach is highly efficient.
In many general flow conditions, the phenomenological Darcy law forms the basis of the theory of flow through porous media. It is a consequence of the linearity of the equations of slow viscous flow which are obtained from the Navier-Stokes equations by neglecting the inertial terms. For time-independent flow conditions and saturated porous media, it is written as
where , , , K, and are all macroscopic variables depending on space vector , is the effective (or kinematic) porosity . This porosity takes into account the volumes of voids effectively concerned with ground-water flow (that is, for instance, without the dead-end pores and the adherence volume of the fluid to the grains). It is upper bounded by total porosity which is the volume occupied by the pores divided by the volume of the bulk medium, is called Darcy's velocity and represents the ground-water flow rate, that is, the volume of water crossing a unit area of porous medium per unit time; is a measurable quantity whereas is the pore velocity, that is, the flow rate per unit area of fluid (which is equivalent to consider that only fluid is present). is the hydraulic potential (or pressure head). It is defined by where p is the fluid pressure, the volumetric mass of the fluid, g the gravitational constant, and z the height. In , the kinematic term is always neglected due to small ground-water velocities in common applications. Finally, K, the proportionality coefficient between Darcy's velocity and the gradient of the hydraulic potential, is called hydraulic conductivity . This parameter (and also permeability) is recognized as a key parameter for ground-water flow. Several experimental techniques (of which mainly: pumping tests, sedimentary analysis, but also seismic, geoelectrical or tracer methods) are used to intensively measure this parameter in the laboratory or in the field scale. These measurements have put in evidence the (highly) heterogeneous behavior of K in space and have suggested the use of stochastic models.
In the stochastic models, the hydrogeological parameters (like K and ) are represented by Random Space Functions (RSF). A RSF h is regarded as a random variable with an infinite number of components. Thus, with the random function depending on the spatial coordinate , hi is defined as the value of h at a point r = ri and the joint probability distribution function (Pdf) ,with , contains the probability information about h. We can also regard as an ordinary function of the space coordinates in each realization and the ensemble underlying F, the Pdf, is that of the values of h at the set of points ().
Law () was the first to apply RSF in porous media and to propose, on the basis of core analysis data from a carbonate oil field reservoir, a log-normal probability density function (Pdf) for K. Since this proposition, there is now a large body of direct evidence to support the statement that the Pdf for hydraulic conductivity is log-normal (), (). Hydraulic log-conductivity is therefore commonly used and assumed to be distributed according to a Gaussian distribution , where mY = <Y>, and is the standard deviation.
Solving the Darcy equation with the random (log-normal) hydraulic conductivity, we have found the statistical characteristics required for the construction of the stochastic Lagrangian model in the form of a Langevin equation. The main difficulty was the unique derivation of this equation.
The results of the numerical simulation on the basis of the developed stochastic Lagrangian model are highly efficient, but what is especially important, is that the new model is free of many limiting assumptions used in the conventional models. In addition, the method calculates not only the mean concentration but also the flux of the concentrations.
2. A probabilistic approach to the solution of boundary value problems. Numerics for problems of stochastic dynamics (G.N. Milstein).
A probabilistic approach to constructing new numerical methods for solving Cauchy and boundary value problems for nonlinear parabolic PDEs is based on making use of the well-known probabilistic representations of solutions to linear PDEs and ideas of the weak sense numerical integration of SDEs. It takes into account a coefficients dependence on the space variables and a relationship between diffusion and advection in an intrinsic manner and allows to obtain a number of new effective numerical algorithms. The approach is applied to the numerical solution of the Dirichlet problem for nonlinear parabolic equations in .
In , some random walks for the general Dirichlet problem for linear elliptic and parabolic PDEs are proposed. They are based on the weak Euler approximation and on linear interpolation. Due to this, it is possible to suggest a number of efficient Monte Carlo algorithms.
Stochastic systems, phase flows of which have integral invariants, are considered in , , and . In particular, the Hamiltonian systems perturbed by additive noise belong to systems with integral invariants. For such systems, numerical methods preserving a number of important features of the original phase flows are constructed. They demonstrate superiority in comparison with nonsymplectic methods in numerical experiments (see Fig. 1).
Some theoretical and numerical investigations of stochastic systems connected with such phenomena as stochastic resonance and noise-induced unidirectional transport are reported in .
Stability properties of stochastic systems are studied in  and .
3. Monte Carlo and Quasi-Monte Carlo methods for integration (P. Mathé).
Many physically relevant quantities are means or require to compute certain means intermediately. Therefore, we investigate numerical methods for the efficient computation of integrals. Precisely, for a fixed probability on some space X, we aim at approximating by means of a sample mean of random variables arising from an ergodic Markov chain with transition kernel K, having as its invariant distribution. Thus, if is an initial distribution, then we use, for a given f on X, the sample mean
where are the consecutive steps of our Markov chain. Typically, the error is measured in mean square sense, i.e.
denotes the individual error of the sample mean at function f. Previous analysis, see , was restricted to uniformly ergodic Markov chains. It is, however, known that on general state space, most Markov chains will not be uniformly ergodic, see . Therefore, as in , we extended our study to V-uniformly ergodic Markov chains. Another goal was to establish the interaction between ergodicity properties of the kernel K, respective classes F of functions to be integrated, and initial distributions .
Let be a real valued function on X. By we denote the Banach space of all functions f on X, for which . The Harris-Markov chain K is called V-uniformly ergodic, if, for the transition operator , we have
We established a series of consequences of this definition, using tools from interpolation theory. The intrinsic rôle of the weight function V became transparent.
Convergence results for the sample mean towards the integral are presented in terms of the following (quadratic) functional
A typical result can be stated as follows. For all initial distributions with and bounded sets F in the weighted Hilbert space L2(V) we have
The convergence results are then extended to classes of functions which are no longer square integrable.
A further aspect was the study of QMC methods for functions over unbounded domains, in particular integrals over with probability weight function . These can be evaluated with the aid of QMC algorithms using a proper decomposition of the domain and arrangement of the low discrepancy points over a series of hierarchical hypercubes. Precisely, we proposed
where m+1 is any number of cubes and , where each nj is the number of points used within cube , and the are the indicator functions of . Furthermore,
The total amount of function evaluations is . The Pj s are some low-discrepancy point sets, consisting of nj points, respectively.
By arranging the number of hypercubes as well as the number of low-discrepancy points properly, the proposed method can be seen to be optimal for power/exponential decaying weights. This was joint work with G. Wei, Hong Kong Baptist University, see . Emphasis was on results which allow application to mathematical finance. Therefore we could assume smooth weights, but less smoothness for the integrands, in order to cover the evaluation of financial products.