Search

The quantum thermal average plays a central role in describing the thermodynamic properties of a quantum system. From the computational perspective, the quantum thermal average can be computed by the path integral molecular dynamics (PIMD), but the knowledge on the quantitative convergence of such approximations is lacking. We propose an alternative computational framework named the continuous loop path integral molecular dynamics (CL-PIMD), which replaces the ring polymer beads by a continuous loop in the spirit of the Feynman--Kac formula. By truncating the number of normal modes to a finite integer $N\in\mathbb N$, we quantify the discrepancy of the statistical average of the truncated CL-PIMD from the true quantum thermal average, and prove that the truncated CL-PIMD has uniform-in-$N$ geometric ergodicity. These results show that the CL-PIMD provides an accurate approximation to the quantum thermal average, and serves as a mathematical justification of the PIMD methodology.

The voltage-conductance kinetic equation for an ensemble of neurons has been studied by many scientists and mathematicians, while its rigorous analysis is still at a premature stage. In this work, we obtain for the first time the exponential convergence to the steady state of this kinetic model in the linear setting. Our proof is based on a probabilistic reformulation, which allows us to investigate microscopic trajectories and bypass the difficulties raised by the special velocity field and boundary conditions in the macroscopic equation. We construct an associated stochastic process, for which proving the minorization condition becomes tractable, and the exponential ergodicity is then proved using Harris' theorem., 23 pages, 3 figures

We investigate the tumor boundary instability induced by nutrient consumption and supply based on a Hele-Shaw model derived from taking the incompressible limit of a cell density model. We analyze the boundary stability/instability in two scenarios: 1) the front of the traveling wave; 2) the radially symmetric boundary. In each scenario, we investigate the boundary behaviors under two different nutrient supply regimes, in vitro, and in vivo. Our main conclusion is that for either scenario, the in vitro regime always stabilizes the tumor's boundary regardless of the nutrient consumption rate. However, boundary instability may occur when the tumor cells aggressively consume nutrients, and the nutrient supply is governed by the in vivo regime., Comment: 33 pages, 10 figures

In this work, we consider the Fokker-Planck equation of the Nonlinear Noisy Leaky Integrate-and-Fire (NNLIF) model for neuron networks. Due to the firing events of neurons at the microscopic level, this Fokker-Planck equation contains dynamic boundary conditions involving specific internal points. To efficiently solve this problem and explore the properties of the unknown, we construct a flexible numerical scheme for the Fokker-Planck equation in the framework of spectral methods that can accurately handle the dynamic boundary condition. This numerical scheme is stable with suitable choices of test function spaces, and asymptotic preserving, and it is easily extendable to variant models with multiple time scales. We also present extensive numerical examples to verify the scheme properties, including order of convergence and time efficiency, and explore unique properties of the model, including blow-up phenomena for the NNLIF model and learning and discriminative properties for the NNLIF model with learning rules.

In neuroscience, the time elapsed since the last discharge has been used to predict the probability of the next discharge. Such predictions can be improved taking into account the last two discharge times, and possibly more. Such multi-time processes arise in many other areas and there is no universal limitation on the number of times to be used. This observation leads us to study the infinite-times renewal equation as a simple model to understand the meaning and properties of such partial differential equations depending on an infinite number of variables.We define two notions of solutions, prove existence and uniqueness of solutions, possibly measures. We also prove the long time convergence, with exponential rate, to the steady state in different, strong or weak, topologies depending on assumptions on the coefficients.

In this paper, we use the Bayesian inversion approach to study the data assimilation problem for a family of tumor growth models described by porous-medium type equations. The models contain uncertain parameters and are indexed by a physical parameter $m$, which characterizes the constitutive relation between density and pressure. Based on these models, we employ the Bayesian inversion framework to infer parametric and nonparametric unknowns that affect tumor growth from noisy observations of tumor cell density. We establish the well-posedness and the stability theories for the Bayesian inversion problem and further prove the convergence of the posterior distribution in the so-called incompressible limit, $m \rightarrow \infty$. Since the posterior distribution across the index regime $m\in[2,\infty)$ can thus be treated in a unified manner, such theoretical results also guide the design of the numerical inference for the unknown. We propose a generic computational framework for such inverse problems, which consists of a typical sampling algorithm and an asymptotic preserving solver for the forward problem. With extensive numerical tests, we demonstrate that the proposed method achieves satisfactory accuracy in the Bayesian inference of the tumor growth models, which is uniform with respect to the constitutive relation., Comment: 29 pages, 14 figures

The noisy leaky integrate-and-fire (NLIF) model describes the voltage configurations of neuron networks with an interacting many-particles system at a microscopic level. When simulating neuron networks of large sizes, computing a coarse-grained mean-field Fokker-Planck equation solving the voltage densities of the networks at a macroscopic level practically serves as a feasible alternative in its high efficiency and credible accuracy. However, the macroscopic model fails to yield valid results of the networks when simulating considerably synchronous networks with active firing events. In this paper, we propose a multi-scale solver for the NLIF networks, which inherits the low cost of the macroscopic solver and the high reliability of the microscopic solver. For each temporal step, the multi-scale solver uses the macroscopic solver when the firing rate of the simulated network is low, while it switches to the microscopic solver when the firing rate tends to blow up. Moreover, the macroscopic and microscopic solvers are integrated with a high-precision switching algorithm to ensure the accuracy of the multi-scale solver. The validity of the multi-scale solver is analyzed from two perspectives: firstly, we provide practically sufficient conditions that guarantee the mean-field approximation of the macroscopic model and present rigorous numerical analysis on simulation errors when coupling the two solvers; secondly, the numerical performance of the multi-scale solver is validated through simulating several large neuron networks, including networks with either instantaneous or periodic input currents which prompt active firing events over a period of time., Comment: 25 Pages, 18 Figures

It has been noticed that when the waiting time distribution exhibits a transition from an intermediate time power law decay to a long-time exponential decay in the continuous time random walk model, a transition from anomalous diffusion to normal diffusion can be observed at the population level. However, the mechanism behind the transition of waiting time distribution is rarely studied. In this paper, we provide one possible mechanism to explain the origin of such transition. A jump model terminated by a state-dependent Poisson clock is studied by a formal asymptotic analysis for the time evolutionary equation of its probability density function. The waiting time behavior under a more relaxed setting can be rigorously characterized by probability tools. Both approaches show the transition phenomenon of the waiting time $T$, which is further verified numerically by particle simulations. Our results indicate that a small drift and strong noise in the state equation and a stiff response in the Poisson rate are crucial to the transitional phenomena.

In this article, we introduce the frozen Gaussian sampling (FGS) algorithm to solve the scalar wave equation in the high-frequency regime. The FGS algorithm is a Monte Carlo sampling strategy based on the frozen Gaussian approximation, which greatly reduces the computation workload in the wave propagation and reconstruction. In this work, we propose feasible and detailed procedures to implement the FGS algorithm to approximate scalar wave equations with Gaussian initial conditions and WKB initial conditions respectively. For both initial data cases, we rigorously analyze the error of applying this algorithm to wave equations of dimensionality $d \geq 3$. In Gaussian initial data cases, we prove that the sampling error due to the Monte Carlo method is independent of the typical wave number. We also derive a quantitative bound of the sampling error in WKB initial data cases. Finally, we validate the performance of the FGS and the theoretical estimates about the sampling error through various numerical examples, which include using the FGS to solve wave equations with both Gaussian and WKB initial data of dimensionality $d = 1, 2$, and $3$.

The nonlinear noisy leaky integrate-and-fire (NNLIF) model is a popular mean-field description of a large number of interacting neurons, which has attracted mathematicians to study from various aspects. A core property of this model is the finite time blow-up of the firing rate, which scientifically corresponds to the synchronization of a neuron network, and mathematically prevents the existence of a global classical solution. In this work, we propose a new generalized solution based on reformulating the PDE model with a specific change of variable in time. A firing rate dependent timescale is introduced, in which the transformed equation can be shown to be globally well-posed for any connectivity parameter even in the event of the blow-up. The generalized solution is then defined via the backward change of timescale, and it may have a jump when the firing rate blows up. We establish properties of the generalized solution including the characterization of blow-ups and the global well-posedness in the original timescale. The generalized solution provides a new perspective to understand the dynamics when the firing rate blows up as well as the continuation of the solution after a blow-up., 43 pages, 3 figures

The random batch method provides an efficient algorithm for computing statistical properties of a canonical ensemble of interacting particles. In this work, we study the error estimates of the fully discrete random batch method, especially in terms of approximating the invariant distribution. Using a triangle inequality framework, we show that the long-time error of the method is $O(\sqrt{\tau} + e^{-\lambda t})$, where $\tau$ is the time step and $\lambda$ is the convergence rate which does not depend on the time step $\tau$ or the number of particles $N$. Our results also apply to the McKean-Vlasov process, which is the mean-field limit of the interacting particle system as the number of particles $N\rightarrow\infty$.

The voltage-conductance kinetic model for the collective behavior of neurons has been studied by scientists and mathematicians for two decades, but the rigorous analysis of its solution structure has been only partially obtained in spite of plenty of numerical evidence in various scenarios. In this work, we consider a simplified voltage-conductance model in which the velocity field in the voltage variable is in a separable form. The long time behavior of the simplified model is fully investigated leading to the following dichotomy: either the density function converges to the global equilibrium, or the firing rate diverges as time goes to infinity. Besides, the fast conductance asymptotic limit is justified and analyzed, where the solution to the limit model either blows up in finite time, or globally exists leading to time periodic solutions. An important implication of these results is that the non-separable velocity field, or physically the leaky mechanism, is a key element for the emergence of periodic solutions in the original model based on the available numerical evidence.

The voltage-conductance equation determines the probability distribution of a stochastic process describing a fluctuation-driven neuronal network arising in the visual cortex. Its structure and degeneracy share many common features with the kinetic Fokker-Planck equation, which has attracted much attention recently. We prove an L $\infty$ bound on the steady state solution and the long term convergence of the evolution problem towards this stationary state. Despite the hypoellipticity property, the difficulty is to treat the boundary conditions that change type along the boundary. This leads us to use specific weights in the Alikakos iterations and adapt the relative entropy method.

We study the geometric ergodicity and the long-time behavior of the Random Batch Method for interacting particle systems, which exhibits superior numerical performance in recent large-scale scientific computing experiments. We show that for both the interacting particle system (IPS) and the random batch interacting particle system (RB–IPS), the distribution laws converge to their respective invariant distributions exponentially, and the convergence rate does not depend on the number of particles [Formula: see text], the time step [Formula: see text] for batch divisions or the batch size [Formula: see text]. Moreover, the Wasserstein-1 distance between the invariant distributions of the IPS and the RB–IPS is bounded by [Formula: see text], showing that the RB–IPS can be used to sample the invariant distribution of the IPS accurately with greatly reduced computational cost.

We study a two-speed model with variant drift fields, which generalizes the Goldstein-Taylor model but the two advection fields are not necessarily in proportion. Due to the lack of a local equilibrium structure of the steady state, prevailing hypocoercivity methods could not be directly applied to such a system. To prove the exponential convergence to equilibrium for this model, we use a recent generalization of the Gearhart-Pr\"uss theorem [Sci. China Math. 64 (2021), no. 3, 507-518] to gain semigroup bounds, where the relative entropy estimate plays a central role in gaining desired resolvent estimates via a compactness argument., Comment: 11 pages

The stochastic growth-fragmentation model describes the temporal evolution of a structured cell population through a discrete-time and continuous-state Markov chain. The simulations of this stochastic process and its invariant measure are of interest. In this paper, we propose a numerical scheme for both the simulation of the process and the computation of the invariant measure, and show that under appropriate assumptions, the numerical chain converges to the continuous growth-fragmentation chain with an explicit error bound. With a triangle inequality argument, we are also able to quantitatively estimate the distance between the invariant measures of these two Markov chains.

In this paper, we investigate fast algorithms in the small fraction order regime to approximate the Caputo derivative $^C_0D_t^\alpha u(t)$ when $\alpha$ is small. We focus on two fast algorithms, i.e. FIR and FIDR, both relying on the sum-of-exponential approximation to reduce the cost of evaluating the history part. FIR is the numerical scheme originally proposed in [16], and FIDR is an alternative scheme proposed in [26], and we show that the latter is superior when $\alpha$ is small. With quantitative estimates, we prove that given a certain error threshold, the computational cost of evaluating the history part of the Caputo derivative can be decreased as $\alpha$ gets small. Hence, only minimal cost for the fast evaluation is required in the small $\alpha$ regime, which matches prevailing protocols in engineering practice. We also present improved stability and error analysis of FIDR for solving linear fractional diffusion equations, which achieves clear dependence of the error bound on the fraction order $\alpha$. Finally, we carry out systematic numerical studies for the performances of both FIR and FIDR schemes, where we explore the trade-off between accuracy and efficiency when $\alpha$ is small., Comment: 42 pages, 5 figures, submitted to Journal of Scientific Computing

A decentralized blockchain is a distributed ledger that is often used as a platform for exchanging goods and services. This ledger is maintained by a network of nodes that obeys a set of rules, called a consensus protocol, which helps to resolve inconsistencies among local copies of a blockchain. In this paper, we build a mathematical framework for the consensus protocol designer, specifying (a) the measurement of a resource which nodes strategically invest in and compete for to win the right to build new blocks in the blockchain; and (b) a payoff function for such efforts. Thus, the equilibrium of an associated stochastic differential game can be implemented by selecting nodes in proportion to this specified resource and penalizing dishonest nodes by its loss. This associated, induced game can be further analyzed using mean field games. The problem can be broken down into two coupled PDEs, where an individual node's optimal control path is solved using a Hamilton-Jacobi-Bellman equation, and where the evolution of states distribution is characterized by a Fokker-Planck equation. We develop numerical methods to compute the mean field equilibrium for both steady states at the infinite time horizon and evolutionary dynamics. As an example, we show how the mean field equilibrium can be applied to the Bitcoin blockchain mechanism design. We demonstrate that a blockchain can be viewed as a mechanism that operates in a decentralized setup and propagates properties of the mean field equilibrium over time, such as the underlying security of the blockchain., Comment: 28 pages, 13 figures

In this paper, we develop a Monte Carlo algorithm named the Frozen Gaussian Sampling (FGS) to solve the semiclassical Schr��dinger equation based on the frozen Gaussian approximation. Due to the highly oscillatory structure of the wave function, traditional mesh-based algorithms suffer from "the curse of dimensionality", which gives rise to more severe computational burden when the semiclassical parameter \(\ep\) is small. The Frozen Gaussian sampling outperforms the existing algorithms in that it is mesh-free in computing the physical observables and is suitable for high dimensional problems. In this work, we provide detailed procedures to implement the FGS for both Gaussian and WKB initial data cases, where the sampling strategies on the phase space balance the need of variance reduction and sampling convenience. Moreover, we rigorously prove that, to reach a certain accuracy, the number of samples needed for the FGS is independent of the scaling parameter \(\ep\). Furthermore, the complexity of the FGS algorithm is of a sublinear scaling with respect to the microscopic degrees of freedom and, in particular, is insensitive to the dimension number. The performance of the FGS is validated through several typical numerical experiments, including simulating scattering by the barrier potential, formation of the caustics and computing the high-dimensional physical observables without mesh., 41 pages, 7 figures

In this article, we develop and analyse a new spectral method to solve the semi-classical Schr��dinger equation based on the Gaussian wave-packet transform (GWPT) and Hagedorn's semi-classical wave-packets (HWP). The GWPT equivalently recasts the highly oscillatory wave equation as a much less oscillatory one (the $w$ equation) coupled with a set of ordinary differential equations governing the dynamics of the so-called GWPT parameters. The Hamiltonian of the $ w $ equation consists of a quadratic part and a small non-quadratic perturbation, which is of order $ \mathcal{O}(\sqrt{\varepsilon }) $, where $ \varepsilon\ll 1 $ is the rescaled Planck's constant. By expanding the solution of the $ w $ equation as a superposition of Hagedorn's wave-packets, we construct a spectral method while the $ \mathcal{O}(\sqrt{\varepsilon}) $ perturbation part is treated by the Galerkin approximation. This numerical implementation of the GWPT avoids imposing artificial boundary conditions and facilitates rigorous numerical analysis. For arbitrary dimensional cases, we establish how the error of solving the semi-classical Schr��dinger equation with the GWPT is determined by the errors of solving the $ w $ equation and the GWPT parameters. We prove that this scheme has the spectral convergence with respect to the number of Hagedorn's wave-packets in one dimension. Extensive numerical tests are provided to demonstrate the properties of the proposed method., A revised version. An appendix is added to establish estimate (3.15) and (3.20)

In this paper, we propose a tumor growth model to incorporate and investigate the spatial effects of autophagy. The cells are classified into two phases: normal cells and autophagic cells, whose dynamics are also coupled with the nutrients. First, we construct a reaction-(cross-)diffusion system describing the evolution of cell densities, where the drift is determined by the negative gradient of the joint pressure, and the reaction terms manifest the unique mechanism of autophagy. Next, in the incompressible limit, such a cell density model naturally connects to a free boundary system, describing the geometric motion of the tumor region. Analyzing the free boundary model in a special case, we show that the ratio of the two phases of cells exponentially converges to a ``well-mixed" limit. Within this ``well-mixed" limit, we obtain an analytical solution of the free boundary system which indicates the exponential growth of the tumor size in the presence of autophagy in contrast to the linear growth without it. Numerical simulations are also provided to illustrate the analytical properties and to explore more scenarios.

In this paper, we consider the field model for complex ionic fluids with an energy variational structure, and analyze the well-posedness to this model with regularized kernels. Furthermore, we deduce the estimate of the maximal density function to quantify the finite size effect. On the numerical side, we adopt a finite volume scheme to the field model, which satisfies the following properties: positivity-preserving, mass conservation and energy dissipation. Besides, series of numerical experiments are provided to demonstrate the properties of the steady state and the finite size effect by showing the equilibrium profiles with different values of the parameter in the kernel., 25 pages, 14 figures

At the continuous level, we consider two types of tumor growth models: the cell density model, which is based on the fluid mechanical construction, is more favorable for scientific interpretation and numerical simulations; and the free boundary model, as the incompressible limit of the former, is more tractable when investigating the boundary propagation. In this work, we aim to investigate the boundary propagation speeds in those models based on asymptotic analysis of the free boundary model and efficient numerical simulations of the cell density model. We derive, for the first time, some analytical solutions for the free boundary model with pressures jumps across the tumor boundary in multi-dimensions with finite tumor sizes. We further show that in the large radius limit, the analytical solutions to the free boundary model in one and multiple spatial dimensions converge to traveling wave solutions. The convergence rate in the propagation speeds are algebraic in multi-dimensions as opposed to the exponential convergence in 1D. We also propose an accurate front capturing numerical scheme for the cell density model, and extensive numerical tests are provided to verify the analytical findings.

There is a growing interest in investigating numerical approximations of the water wave equation in recent years, whereas the lack of rigorous analysis of its time discretization inhibits the design of more efficient algorithms. In this work, we focus on a nonlocal hyperbolic model, which essentially inherits the features of the water wave equation, and is simplified from the latter. For the constant coefficient case, we carry out systematical stability studies of the fully discrete approximation of such systems with the Fourier spectral approximation in space and general Runge-Kutta method in time. In particular, we discover the optimal time step constraints, in the form of a modified CFL condition, when certain explicit Runge-Kutta method are applied. Besides, the convergence of the semi-discrete approximation of variable coefficient case is shown, which naturally connects to the water wave equation. Extensive numerical tests have been performed to verify the stability conditions and simulations of the simplified hyperbolic model in the high frequency regime and the water wave equation are also provided.

We propose a semi-discrete scheme for 2D Keller-Segel equations based on a symmetrization reformation, which is equivalent to the convex splitting method and is free of any nonlinear solver. We show that, this new scheme is unconditionally stable as long as the initial condition does not exceed certain threshold, and it asymptotically preserves the quasi-static limit in the transient regime. Furthermore, we prove that the fully discrete scheme is conservative and positivity preserving, which makes it ideal for simulations. The analogical schemes for the radial symmetric cases and the subcritical degenerate cases are also presented and analyzed. With extensive numerical tests, we verify the claimed properties of the methods and demonstrate their superiority in various challenging applications.