Search

In recent years, tremendous progress has been made on numerical algorithms for solving partial differential equations (PDEs) in a very high dimension, using ideas from either nonlinear (multilevel) Monte Carlo or deep learning. They are potentially free of the curse of dimensionality for many different applications and have been proven to be so in the case of some nonlinear Monte Carlo methods for nonlinear parabolic PDEs. In this paper, we review these numerical and theoretical advances. In addition to algorithms based on stochastic reformulations of the original problem, such as the multilevel Picard iteration and the deep backward stochastic differential equations method, we also discuss algorithms based on the more traditional Ritz, Galerkin, and least square formulations. We hope to demonstrate to the reader that studying PDEs as well as control and variational problems in very high dimensions might very well be among the most promising new directions in mathematics and scientific computing in the near future.

One of the key issues in the analysis of machine learning models is to identify the appropriate function space and norm for the model. This is the set of functions endowed with a quantity which can control the approximation and estimation errors by a particular machine learning model. In this paper, we address this issue for two representative neural network models: the two-layer networks and the residual neural networks. We define the Barron space and show that it is the right space for two-layer neural network models in the sense that optimal direct and inverse approximation theorems hold for functions in the Barron space. For residual neural network models, we construct the so-called flow-induced function space, and prove direct and inverse approximation theorems for this space. In addition, we show that the Rademacher complexity for bounded sets under these norms has the optimal upper bounds.

We develop Banach spaces for ReLU neural networks of finite depth $L$ and infinite width. The spaces contain all finite fully connected $L$-layer networks and their $L^2$-limiting objects under bounds on the natural path-norm. Under this norm, the unit ball in the space for $L$-layer networks has low Rademacher complexity and thus favorable generalization properties. Functions in these spaces can be approximated by multi-layer neural networks with dimension-independent convergence rates. The key to this work is a new way of representing functions in some form of expectations, motivated by multi-layer neural networks. This representation allows us to define a new class of continuous models for machine learning. We show that the gradient flow defined this way is the natural continuous analog of the gradient descent dynamics for the associated multi-layer neural networks. We show that the path-norm increases at most polynomially under this continuous gradient flow dynamics.

A new framework is introduced for constructing interpretable and truly reliable reduced models for multiscale problems in situations without scale separation. Hydrodynamic approximation to the kinetic equation is used as an example to illustrate the main steps and issues involved. To this end, a set of generalized moments are constructed first to optimally represent the underlying velocity distribution. The well-known closure problem is then solved with the aim of best capturing the associated dynamics of the kinetic equation. The issue of physical constraints such as Galilean invariance is addressed and an active learning procedure is introduced to help ensure that the dataset used is representative enough. The reduced system takes the form of a conventional moment system and works regardless of the numerical discretization used. Numerical results are presented for the BGK (Bhatnagar--Gross--Krook) model and binary collision of Maxwell molecules. We demonstrate that the reduced model achieves a uniform accuracy in a wide range of Knudsen numbers spanning from the hydrodynamic limit to free molecular flow., Comment: Proceedings of the National Academy of Sciences (2019)

High-dimensional partial differential equations (PDE) appear in a number of models from the financial industry, such as in derivative pricing models, credit valuation adjustment (CVA) models, or portfolio optimization models. The PDEs in such applications are high-dimensional as the dimension corresponds to the number of financial assets in a portfolio. Moreover, such PDEs are often fully nonlinear due to the need to incorporate certain nonlinear phenomena in the model such as default risks, transaction costs, volatility uncertainty (Knightian uncertainty), or trading constraints in the model. Such high-dimensional fully nonlinear PDEs are exceedingly difficult to solve as the computational effort for standard approximation methods grows exponentially with the dimension. In this work we propose a new method for solving high-dimensional fully nonlinear second-order PDEs. Our method can in particular be used to sample from high-dimensional nonlinear expectations. The method is based on (i) a connection between fully nonlinear second-order PDEs and second-order backward stochastic differential equations (2BSDEs), (ii) a merged formulation of the PDE and the 2BSDE problem, (iii) a temporal forward discretization of the 2BSDE and a spatial approximation via deep neural nets, and (iv) a stochastic gradient descent-type optimization procedure. Numerical results obtained using ${\rm T{\small ENSOR}F{\small LOW}}$ in ${\rm P{\small YTHON}}$ illustrate the efficiency and the accuracy of the method in the cases of a $100$-dimensional Black-Scholes-Barenblatt equation, a $100$-dimensional Hamilton-Jacobi-Bellman equation, and a nonlinear expectation of a $ 100 $-dimensional $ G $-Brownian motion., 56 pages, 12 figures

New estimates for the population risk are established for two-layer neural networks. These estimates are nearly optimal in the sense that the error rates scale in the same way as the Monte Carlo error rates. They are equally effective in the over-parametrized regime when the network size is much larger than the size of the dataset. These new estimates are a priori in nature in the sense that the bounds depend only on some norms of the underlying functions to be fitted, not the parameters in the model, in contrast with most existing results which are a posteriori in nature. Using these a priori estimates, we provide a perspective for understanding why two-layer neural networks perform better than the related kernel methods.

We characterize the meaning of words with language-independent numerical fingerprints, through a mathematical analysis of recurring patterns in texts. Approximating texts by Markov processes on a long-range time scale, we are able to extract topics, discover synonyms, and sketch semantic fields from a particular document of moderate length, without consulting external knowledge-base or thesaurus. Our Markov semantic model allows us to represent each topical concept by a low-dimensional vector, interpretable as algebraic invariants in succinct statistical operations on the document, targeting local environments of individual words. These language-independent semantic representations enable a robot reader to both understand short texts in a given language (automated question-answering) and match medium-length texts across different languages (automated word translation). Our semantic fingerprints quantify local meaning of words in 14 representative languages across 5 major language families, suggesting a universal and cost-effective mechanism by which human languages are processed at the semantic level. Our protocols and source codes are publicly available on https://github.com/yajun-zhou/linguae-naturalis-principia-mathematica, Comment: Main text (12 pages, 7 figures); Software Manual (ii+262 pages, 16 figures, 12 tables, available as two ancillary files). Revised according to reviewers' comments

We establish a scale separation of Kolmogorov width type between subspaces of a given Banach space under the condition that a sequence of linear maps converges much faster on one of the subspaces. The general technique is then applied to show that reproducing kernel Hilbert spaces are poor $$L^{2}$$ -approximators for the class of two-layer neural networks in high dimension, and that multi-layer networks with small path norm are poor approximators for certain Lipschitz functions, also in the $$L^{2}$$ -topology.

We study the natural function space for infinitely wide two-layer neural networks with ReLU activation (Barron space) and establish different representation formulae. In two cases, we describe the space explicitly up to isomorphism. Using a convenient representation, we study the pointwise properties of two-layer networks and show that functions whose singular set is fractal or curved (for example distance functions from smooth submanifolds) cannot be represented by infinitely wide two-layer networks with finite path-norm. We use this structure theorem to show that the only $C^1$-diffeomorphisms which Barron space are affine. Furthermore, we show that every Barron function can be decomposed as the sum of a bounded and a positively one-homogeneous function and that there exist Barron functions which decay rapidly at infinity and are globally Lebesgue-integrable. This result suggests that two-layer neural networks may be able to approximate a greater variety of functions than commonly believed.

We prove that for analytic functions in low dimension, the convergence rate of the deep neural network approximation is exponential.

Developing algorithms for solving high-dimensional partial differential equations (PDEs) has been an exceedingly difficult task for a long time, due to the notoriously difficult problem known as the "curse of dimensionality". This paper introduces a deep learning-based approach that can handle general high-dimensional parabolic PDEs. To this end, the PDEs are reformulated using backward stochastic differential equations and the gradient of the unknown solution is approximated by neural networks, very much in the spirit of deep reinforcement learning with the gradient acting as the policy function. Numerical results on examples including the nonlinear Black-Scholes equation, the Hamilton-Jacobi-Bellman equation, and the Allen-Cahn equation suggest that the proposed algorithm is quite effective in high dimensions, in terms of both accuracy and cost. This opens up new possibilities in economics, finance, operational research, and physics, by considering all participating agents, assets, resources, or particles together at the same time, instead of making ad hoc assumptions on their inter-relationships., 13 pages, 6 figures

We propose a deep learning-based method, the Deep Ritz Method, for numerically solving variational problems, particularly the ones that arise from partial differential equations. The Deep Ritz Method is naturally nonlinear, naturally adaptive and has the potential to work in rather high dimensions. The framework is quite simple and fits well with the stochastic gradient descent method used in deep learning. We illustrate the method on several problems including some eigenvalue problems.

We propose a new algorithm for solving parabolic partial differential equations (PDEs) and backward stochastic differential equations (BSDEs) in high dimension, by making an analogy between the BSDE and reinforcement learning with the gradient of the solution playing the role of the policy function, and the loss function given by the error between the prescribed terminal condition and the solution of the BSDE. The policy function is then approximated by a neural network, as is done in deep reinforcement learning. Numerical results using TensorFlow illustrate the efficiency and accuracy of the proposed algorithms for several 100-dimensional nonlinear PDEs from physics and finance such as the Allen-Cahn equation, the Hamilton-Jacobi-Bellman equation, and a nonlinear pricing model for financial derivatives., 39 pages, 15 figures

We discuss the idea of using continuous dynamical systems to model general high-dimensional nonlinear functions used in machine learning. We also discuss the connection with deep learning.

This article analyzes well-definedness and regularity of renormalized powers of Ornstein-Uhlenbeck processes and uses this analysis to establish local existence, uniqueness and regularity of strong solutions of stochastic Ginzburg-Landau equations with polynomial nonlinearities in two space dimensions and with quadratic nonlinearities in three space dimensions., Comment: 42 pages

The instability of the two-dimensional Poiseuille flow in a long channel and the subsequent transition is studied using a thermodynamic approach. The idea is to view the transition process as an initial value problem with the initial condition being Poiseuille flow plus noise, which is considered as our ensemble. Using the mean energy of the velocity fluctuation and the skin friction coefficient as the macrostate variable, we analyze the transition process triggered by the initial noises with different amplitudes. A first order transition is observed at the critical Reynolds number $Re_* \sim 5772$ in the limit of zero noise. An action function, which relates the mean energy with the noise amplitude, is defined and computed. The action function depends only on the Reynolds number, and represents the cost for the noise to trigger a transition from the laminar flow. The correlation function of the spatial structure is analyzed.

We introduce a new family of trial wave-functions based on deep neural networks to solve the many-electron Schrodinger equation. The Pauli exclusion principle is dealt with explicitly to ensure that the trial wave-functions are physical. The optimal trial wave-function is obtained through variational Monte Carlo and the computational cost scales quadratically with the number of electrons. The algorithm does not make use of any prior knowledge such as atomic orbitals. Yet it is able to represent accurately the ground-states of the tested systems, including He, H2, Be, B, LiH, and a chain of 10 hydrogen atoms. This opens up new possibilities for solving large-scale many-electron Schrodinger equation.

Recent work linking deep neural networks and dynamical systems opened up new avenues to analyze deep learning. In particular, it is observed that new insights can be obtained by recasting deep learning as an optimal control problem on difference or differential equations. However, the mathematical aspects of such a formulation have not been systematically explored. This paper introduces the mathematical formulation of the population risk minimization problem in deep learning as a mean-field optimal control problem. Mirroring the development of classical optimal control, we state and prove optimality conditions of both the Hamilton-Jacobi-Bellman type and the Pontryagin type. These mean-field results reflect the probabilistic nature of the learning problem. In addition, by appealing to the mean-field Pontryagin's maximum principle, we establish some quantitative relationships between population and empirical learning problems. This serves to establish a mathematical foundation for investigating the algorithmic and theoretical connections between optimal control and deep learning., Comment: 44 pages

In this work, we study the nonlinear instability of two-dimensional (2D) wall-bounded shear flows from the large deviation point of view. The main idea is to consider the Navier–Stokes equations perturbed by small noise in force and then examine the noise-induced transitions between the two coexisting stable solutions due to the subcritical bifurcation. When the amplitude of the noise goes to zero, the Freidlin–Wentzell (F–W) theory of large deviations defines the most probable transition path in the phase space, which is the minimizer of the F–W action functional and characterizes the development of the nonlinear instability subject to small random perturbations. Based on such a transition path we can define a critical Reynolds number for the nonlinear instability in the probabilistic sense. Then the action-based stability theory is applied to study the 2D Poiseuille flow in a short channel.

We study the optimal convergence rate for the universal estimation error. Let \(\mathcal {F}\) be the excess loss class associated with the hypothesis space and n be the size of the data set, we prove that if the Fat-shattering dimension satisfies \(\text {fat}_{\epsilon } (\mathcal {F})= O(\epsilon ^{-p})\), then the universal estimation error is of \(O(n^{-1/2})\) for \(p 2\). Among other things, this result gives a criterion for a hypothesis class to achieve the minimax optimal rate of \(O(n^{-1/2})\). We also show that if the hypothesis space is the compact supported convex Lipschitz continuous functions in \(\mathbb {R}^d\) with \(d>4\), then the rate is approximately \(O(n^{-2/d})\).

Parabolic partial differential equations (PDEs) and backward stochastic differential equations (BSDEs) are key ingredients in a number of models in physics and financial engineering. In particular, parabolic PDEs and BSDEs are fundamental tools in the state-of-the-art pricing and hedging of financial derivatives. The PDEs and BSDEs appearing in such applications are often high-dimensional and nonlinear. Since explicit solutions of such PDEs and BSDEs are typically not available, it is a very active topic of research to solve such PDEs and BSDEs approximately. In the recent article [E, W., Hutzenthaler, M., Jentzen, A., and Kruse, T. Linear scaling algorithms for solving high-dimensional nonlinear parabolic differential equations. arXiv:1607.03295 (2017)] we proposed a family of approximation methods based on Picard approximations and multilevel Monte Carlo methods and showed under suitable regularity assumptions on the exact solution for semilinear heat equations that the computational complexity is bounded by $O( d \, \epsilon^{-(4+\delta)})$ for any $\delta\in(0,\infty)$, where $d$ is the dimensionality of the problem and $\epsilon\in(0,\infty)$ is the prescribed accuracy. In this paper, we test the applicability of this algorithm on a variety of $100$-dimensional nonlinear PDEs that arise in physics and finance by means of numerical simulations presenting approximation accuracy against runtime. The simulation results for these 100-dimensional example PDEs are very satisfactory in terms of accuracy and speed. In addition, we also provide a review of other approximation methods for nonlinear PDEs and BSDEs from the literature.

We present an efficient algorithm for calculating the minimum energy path (MEP) and energy barriers between local minima on a multidimensional potential energy surface (PES). Such paths play a central role in the understanding of transition pathways between metastable states. Our method relies on the original formulation of the string method [Phys. Rev. B, 66,052301 (2002)], i.e. to evolve a smooth curve along a direction normal to the curve. The algorithm works by performing minimization steps on hyperplanes normal to the curve. Therefore the problem of finding MEP on the PES is remodeled as a set of constrained minimization problems. This provides the flexibility of using minimization algorithms faster than the steepest descent method used in the simplified string method [J. Chem. Phys., 126(16), 164103 (2007)]. At the same time, it provides a more direct analog of the finite temperature string method. The applicability of the algorithm is demonstrated using various examples.

We present for the first time an efficient iterative method to directly solve the four-component Dirac-Kohn-Sham (DKS) density functional theory. Due to the existence of the negative energy continuum in the DKS operator, the existing iterative techniques for solving the Kohn-Sham systems cannot be efficiently applied to solve the DKS systems. The key component of our method is a novel filtering step (F) which acts as a preconditioner in the framework of the locally optimal block preconditioned conjugate gradient (LOBPCG) method. The resulting method, dubbed the LOBPCG-F method, is able to compute the desired eigenvalues and eigenvectors in the positive energy band without computing any state in the negative energy band. The LOBPCG-F method introduces mild extra cost compared to the standard LOBPCG method and can be easily implemented. We demonstrate our method in the pseudopotential framework with a planewave basis set which naturally satisfies the kinetic balance prescription. Numerical results for Pt$_{2}$, Au$_{2}$, TlF, and Bi$_{2}$Se$_{3}$ indicate that the LOBPCG-F method is a robust and efficient method for investigating the relativistic effect in systems containing heavy elements., 31 pages, 5 figures

We develop a technique for generating a set of optimized local basis functions to solve models in the Kohn-Sham density functional theory for both insulating and metallic systems. The optimized local basis functions are obtained by solving a minimization problem in an admissible set determined by a large number of primitive basis functions. Using the optimized local basis set, the electron energy and the atomic force can be calculated accurately with a small number of basis functions. The Pulay force is systematically controlled and is not required to be calculated, which makes the optimized local basis set an ideal tool for ab initio molecular dynamics and structure optimization. We also propose a preconditioned Newton-GMRES method to obtain the optimized local basis functions in practice. The optimized local basis set is able to achieve high accuracy with a small number of basis functions per atom when applied to a one dimensional model problem.

A theory for noise-driven subcritical instabilities in spatially extended systems is put forward. The theory allows one to calculate the critical bifurcation parameter for a first-order phase transition in such non-equilibrium systems in the thermodynamic limit and analyse the mechanism of phase transition. Two examples with distinctive features are studied in detail to demonstrate the usefulness of the theory and the different scenarios that can occur in the thermodynamic limit of non-equilibrium systems.

We study the semi-classical limit of the Schrodinger equation in a crystal in the presence of an external potential and magnetic field. We first introduce the Bloch-Wigner transform and derive the asymptotic equations governing this transform in the semi-classical setting. For the second part, we focus on the appearance of the Berry curvature terms in the asymptotic equations. These terms play a crucial role in many important physical phenomena such as the quantum Hall effect. We give a simple derivation of these terms in different settings using asymptotic analysis.

We describe an efficient implementation of an algorithm for computing selected elements of a general sparse symmetric matrix A that can be decomposed as A = LDLT , where L is lower triangular and D is diagonal. Our implementation, which is called SelInv , is built on top of an efficient supernodal left-looking LDLT factorization of A . We discuss how computational efficiency can be gained by making use of a relative index array to handle indirect addressing. We report the performance of SelInv on a collection of sparse matrices of various sizes and nonzero structures. We also demonstrate how SelInv can be used in electronic structure calculations.

The behavior of interacting electrons in a perfect crystal under macroscopic external electric and magnetic fields is studied. Effective Maxwell equations for the macroscopic electric and magnetic fields are derived starting from time-dependent density functional theory. Effective permittivity and permeability coefficients are obtained.

An efficient parallel algorithm is presented for computing selected components of $A^{-1}$ where $A$ is a structured symmetric sparse matrix. Calculations of this type are useful for several applications, including electronic structure analysis of materials in which the diagonal elements of the Green's functions are needed. The algorithm proposed here is a direct method based on a block $LDL^T$ factorization. The selected elements of $A^{-1}$ we compute lie in the nonzero positions of $L+L^T$. We use the elimination tree associated with the block $LDL^T$ factorization to organize the parallel algorithm, and reduce the synchronization overhead by passing the data level by level along this tree using the technique of local buffers and relative indices. We demonstrate the efficiency of our parallel implementation by applying it to a discretized two dimensional Hamiltonian matrix. We analyze the performance of the parallel algorithm by examining its load balance and communication overhead, and show that our parallel implementation exhibits an excellent weak scaling on a large-scale high performance distributed-memory parallel machine.

The electronic structure of a smoothly deformed crystal is analyzed using a minimalist model in quantum many-body theory, the nonlinear tight-binding model. An extension of the classical Cauchy-Born rule for crystal lattices is established for the electronic structure under sharp stability conditions. A nonlinear elasticity model is rigorously derived. The onset of instability is briefly examined. © 2010 Wiley Periodicals, Inc.

We present a multiscale model for numerical simulations of dynamics of crystalline solids. The method combines the continuum nonlinear elasto-dynamics model, which models the stress waves and physical loading conditions, and molecular dynamics model, which provides the nonlinear constitutive relation and resolves the atomic structures near local defects. The coupling of the two models is achieved based on a general framework for multiscale modeling - the heterogeneous multiscale method (HMM). We derive an explicit coupling condition at the atomistic/continuum interface. Application to the dynamics of brittle cracks under various loading conditions is presented as test examples.

A numerical approach based on the string method is developed to study nucleation of ordered phases in first-order phase transitions. Among other things, this method allows an efficient computation of the minimum energy path (MEP) during the nucleation process. The MEP provides information about the size, shape and free energy barrier of the critical nucleus. To improve the efficiency of the string method, a special initialization process is proposed. Constraints from physical models are treated using two methods, a generalized coordinates method and a projection method. Strategies for choosing the computational domain and defining the nucleus boundary are also introduced. The validity of our approach is illustrated by two nontrivial examples from soft condensed matter physics, namely the nematic-isotropic transition of liquid crystals and the ordered-to-ordered phase transition of diblock copolymers.

We investigate noise-induced transitions in non-gradient systems when complex invariant sets emerge. Our example is the Lorenz system in three representative Rayleigh number regimes. It is found that before the homoclinic explosion bifurcation, the only transition state is the saddle point, and the transition is similar to that in gradient systems. However, when the chaotic invariant set emerges, an unstable limit cycle continues from the homoclinic trajectory. This orbit, which is embedded in a local tube-like manifold around the initial stable stationary point as a relative attractor, plays the role of the most probable exit set in the transition process. This example demonstrates how limit cycles, the next simplest invariant set beyond fixed points, can be involved in the transition process in smooth dynamical systems.

We present a new general framework for designing multiscale methods. Compared with previous work such as Brandt's systematic up-scaling, the heterogeneous multiscale method (HMM) and the ''equation-free'' approach, this new framework has the distinct feature that it does not require reinitializing the microscale model at each macro time step or each macro iteration step. In the new strategy, the macro- and micro-models evolve simultaneously using different time steps (and therefore different clocks), and they exchange data at every step. The micro-model uses its own appropriate time step. The macro-model runs at a slower pace than required by accuracy and stability considerations for the macroscale dynamics, in order for the micro-model to relax. Examples are discussed and application to modeling complex fluids is presented.

We propose an algorithm for extracting the diagonal of the inverse matrices arising from electronic structure calculation. The proposed algorithm uses a hierarchical decomposition of the computational domain. It first constructs hierarchical Schur complements of the interior points for the blocks of the domain in a bottom-up pass and then extracts the diagonal entries efficiently in a top-down pass by exploiting the hierarchical local dependence of the inverse matrices. The overall cost of our algorithm is $O(N^3/2)$ for a two dimensional problem with $N$ degrees of freedom. Numerical results in electronic structure calculation illustrate the efficiency and accuracy of the proposed algorithm.

Orbital minimization is among the most promising linear scaling algorithms for electronic structure calculation. However, to achieve linear scaling, one has to truncate the support of the orbitals and this introduces many problems, the most important of which is the occurrence of numerous local minima. In this paper, we introduce a simple modification of the orbital minimization method, by adding a localization step into the algorithm. This localization step selects the most localized representation of the subspace spanned by the orbitals obtained during the intermediate stages of the iteration process. We show that the addition of the localization step substantially reduces the chances that the iterations get trapped at local minima.

The authors consider the simplest quantum mechanics model of solids, the tight binding model, and prove that in the continuum limit, the energy of tight binding model converges to that of the continuum elasticity model obtained using Cauchy-Born rule. The technique in this paper is based mainly on spectral perturbation theory for large matrices.

We study continuum and atomistic models for the elastodynamics of crystalline solids at zero temperature. We establish sharp criterion for the regime of validity of the nonlinear elastic wave equations derived from the well-known Cauchy-Born rule.

In this paper we develop a domain decomposition method (DDM), based on the discontinuous Galerkin (DG) and the local discontinuous Galerkin (LDG) methods, for solving multiscale problems involving macro sub-domains, where a macro model is valid, and micro sub-domains, where the macro model is not valid and a more costly micro model must be used. We take two examples, one from compressible gas dynamics where the micro sub-domains are around shocks, contacts and corners of rarefaction fans, and another one from semiconductor device simulations where the micro sub-domains are around the jumps in the doping profile. The macro model is taken as the Euler equations for the gas dynamics problem and as a hydrodynamic model and a high field model for the semiconductor device problem. The micro model for both problems is taken as a kinetic equation. We pay special attention to the effective coupling between the macro sub-domains and the micro sub-domains, in which we utilize the advantage of the discontinuous Galerkin method in its compactness of the computational stencil. Numerical results demonstrate the effectiveness of our DDM-DG method in solving such multi-scale problems.

Many multiscale methods are based on the idea of extracting macroscopic behavior of solutions by solving an array of microscale models over small domains. A key ingredient in such multiscale methods is the boundary condition and the size of the computational domain over which the microscale problems are solved. This problem is systematically investigated in the present paper in the context of modeling strongly heterogeneous media. Three different boundary conditions are considered: the periodic boundary condition, Dirichlet boundary condition, and the Neumann boundary condition. Each is applied to several benchmark problems: the random checker-board problem, periodic problem with isotropic macroscale behavior, periodic problem with anisotropic macroscale behavior and periodic laminated media. In each case, convergence studies are conducted as the domain size for the microscale problem is changed. Convergence rates as well as the size of fluctuations in the computed effective coefficients are compared for the different formulations. In addition, we will discuss a mixed Dirichlet-Neumann boundary condition that is often used in porous medium modeling. We explain why that leads to unsatisfactory results and how it can be corrected. Also discussed are the different averaging methods used in extracting the effective coefficients.

We present an efficient numerical algorithm for simulating chemical kinetic systems with multiple time scales. This algorithm is an improvement of the traditional stochastic simulation algorithm (SSA), also known as Gillespie's algorithm. It is in the form of a nested SSA and uses an outer SSA to simulate the slow reactions with rates computed from realizations of inner SSAs that simulate the fast reactions. The algorithm itself is quite general and seamless, and it amounts to a small modification of the original SSA. Our analysis of such multi-scale chemical kinetic systems allows us to identify the slow variables in the system, derive effective dynamics on the slow time scale, and provide error estimates for the nested SSA. Efficiency of the nested SSA is discussed using these error estimates, and illustrated through several numerical examples.

We consider the continuum limit for models of solids that arise in density functional theory and the QM-continuum approximation of such models. Two different versions of QM- continuum approximation are proposed, depending on the level at which the Cauchy-Born rule is used, one at the level of electron density and one at the level of energy. Consistency at the interface between the smooth and the non-smooth regions is analyzed. We show that if the Cauchy-Born rule is used at the level of electron density, then the resulting QM-continuum model is free of the so-called “ghost force” at the interface. We also present dynamic models that bridge naturally the Car-Parrinello method and the QM-continuum approximation.

We present a general mathematical framework for modeling the macroscale behavior of a multiscale system using only microscale models, by formulating the effective macroscale models as dynamic models on the underlying fiber bundles. This framework allows us to carry out seamless multiscale modeling using traditional numerical techniques. At the same time, they give rise to an interesting mathematical structure and new interesting mathematical problems. We discuss several examples from homogenization problems, continuum modeling of solids based on atomistic or electronic structure models, macroscale behavior of interacting diffusion, and continuum modeling of complex fluids based on kinetic and Brownian dynamics models.

We carry out a mathematical study of the Andersen thermostat [1], which is a frequently used tool in molecular dynamics. After reformulating the continuous- and discrete-time Andersen dynamics, we prove that in both cases the Andersen dynamics is uniformly ergodic. A detailed numerical analysis is presented, establishing the rate of convergence of most commonly used numerical algorithms for the Andersen thermostat. Transport properties such as the diffusion constant are also investigated. It is proved for the Lorentz gas model where there is intrinsic diffusion that the diffusion coefficient calculated using the Andersen thermostat converges to the true diffusion coefficient in the limit of vanishing collision frequency in the Andersen thermostat. © 2007 Wiley Periodicals, Inc.

We discuss numerical methods for linear and nonlinear transport equations with multiscale velocity fields. These methods are themselves multiscaled in nature in the sense that they use macro and micro grids, multiscale test functions. We demonstrate the efficiency of these methods and apply them to two-phase flow in heterogeneous porous media.

The framework of the heterogeneous multiscale method (HMM) is used to develop numerical methods for the study of macroscale dynamics of fluids in situations, where either the constitutive relation or the boundary conditions are not explicitly available and have to be inferred from microscopic models such as molecular dynamics. Continuum hydrodynamics is used as the macroscopic model, while molecular dynamics serves as the microscopic model and is used to supply the necessary data, e.g., the stress or the boundary condition, for the macroscopic model. Scale separation is exploited so that the macroscopic variables can be evolved in macroscopic spatial/temporal scales using data that are estimated from molecular dynamics simulation on microscale spatial/temporal domains. This naturally decouples the micro and macrospatial and temporal scales whenever possible. Applications are presented for models of complex fluids, contact line dynamics, and a simple model of non-trivial fluid-solid interactions.

We analyze a class of numerical schemes proposed in [26] for stochastic differential equations with multiple time scales. Both advective and diffusive time scales are considered. Weak as well as strong convergence theorems are proven. Most of our results are optimal. They in turn allow us to provide a thorough discussion on the efficiency as well as optimal strategy for the method. c ! 2005 Wiley Periodicals, Inc.