Search

In this paper, we investigate the safety mechanisms of instruction fine-tuned large language models (LLMs). We discover that re-weighting MLP neurons can significantly compromise a model's safety, especially for MLPs in end-of-sentence inferences. We hypothesize that LLMs evaluate the harmfulness of prompts during end-of-sentence inferences, and MLP layers plays a critical role in this process. Based on this hypothesis, we develop 2 novel white-box jailbreak methods: a prompt-specific method and a prompt-general method. The prompt-specific method targets individual prompts and optimizes the attack on the fly, while the prompt-general method is pre-trained offline and can generalize to unseen harmful prompts. Our methods demonstrate robust performance across 7 popular open-source LLMs, size ranging from 2B to 72B. Furthermore, our study provides insights into vulnerabilities of instruction-tuned LLM's safety and deepens the understanding of the internal mechanisms of LLMs.

The integrate and fire equation is a classical model for neural assemblies which can exhibit finite time blow-up. A major open problem is to understand how to continue solutions after blow-up. Here we study an approach based on random discharge models and a change of time which generates a classical global solution to the expense of a strong absorption rate 1/$\epsilon$. We prove that in the limit $\epsilon$ $\rightarrow$ 0 + , a global solution is recovered where the integrate and fire equation is reformulated with a singular measure. This describes the dynamics after blow-up and also gives information on the blow-up phenomena itself.The major difficulty is to handle nonlinear terms. To circumvent it, we establish two new estimates, a kind of equi-integrability of the discharge measure and a L 2 estimate of the density. The use of the new timescale turns out to be fundamental for those estimates.

We propose a novel discrete Poisson equation approach to estimate the statistical error of a broad class of numerical integrators for the underdamped Langevin dynamics. The statistical error refers to the mean square error of the estimator to the exact ensemble average with a finite number of iterations. With the proposed error analysis framework, we show that when the potential function $U(x)$ is strongly convex in $\mathbb R^d$ and the numerical integrator has strong order $p$, the statistical error is $O(h^{2p}+\frac1{Nh})$, where $h$ is the time step and $N$ is the number of iterations. Besides, this approach can be adopted to analyze integrators with stochastic gradients, and quantitative estimates can be derived as well. Our approach only requires the geometric ergodicity of the continuous-time underdamped Langevin dynamics, and relaxes the constraint on the time step.

We consider a special type of fast reaction-diffusion systems in which the coefficients of the reaction terms of the two substances are much larger than those of the diffusion terms while the diffusive motion to the substrate is negligible. Specifically speaking, the rate constants of the reaction terms are $O(1/\epsilon)$ while the diffusion coefficients are $O(1)$ where the parameter $\epsilon$ is small. When the rate constants of the reaction terms become highly large, i.e. $\epsilon$ tends to 0, the singular limit behavior of such a fast reaction-diffusion system is inscribed by the Stefan problem with latent heat, which brings great challenges in numerical simulations. In this paper, we adopt a semi-implicit scheme, which is first-order accurate in time and can accurately approximate the interface propagation even when the reaction becomes extremely fast, that is to say, the parameter $\epsilon$ is sufficiently small. The scheme satisfies the positivity, bound preserving properties and has $L^2$ stability and the linearized stability results of the system. For better performance on numerical simulations, we then construct a semi-implicit Runge-Kutta scheme which is second-order accurate in time. Numerous numerical tests are carried out to demonstrate the properties, such as the order of accuracy, positivity and bound preserving, the capturing of the sharp interface with various $\epsilon$ and to simulate the dynamics of the substances and the substrate, and to explore the heat transfer process, such as solid melting or liquid solidification in two dimensions.

In this paper, we consider a Hele-Shaw model that describes tumor growth subject to nutrient supply. This model was recently studied in \cite{feng2022tumor} via asymptotic analysis. Our contributions are twofold: Firstly, we provide a rigorous derivation of this Hele-Shaw model by taking the incompressible limit of the porous medium reaction-diffusion equation, which solidifies the mathematical foundations of the model. Secondly, from a bifurcation theory perspective, we prove the existence of non-symmetric traveling wave solutions to the model, which reflect the intrinsic boundary instability in tumor growth dynamics., Comment: 23 pages, 2 figures

We introduce a novel reformulation of the mean-field system for pulse-coupled oscillators. It is based on writing a closed equation for the inverse distribution function associated to the probability density of oscillators with a given phase in a suitable time scale. This new framework allows to show a hidden contraction/expansion of certain distances leading to a full clarification of the long-time behavior, existence of steady states, rates of convergence, and finite time blow-up of classical solutions for a large class of monotone phase response functions. In the process, we get insights about the origin of obstructions to global-in-time existence and uniform in time estimates on the firing rate of the oscillators.

We develop an enthalpy-based modeling and computational framework to quantify uncertainty in Stefan problems with an injection boundary. Inspired by airfoil icing studies, we consider a system featuring an injection boundary inducing domain changes and a free boundary separating phases, resulting in two types of moving boundaries. Our proposed enthalpy-based formulation seamlessly integrates thermal diffusion across the domain with energy fluxes at the boundaries, addressing a modified injection condition for boundary movement. Uncertainty then stems from random variations in the injection boundary. The primary focus of our Uncertainty Quantification (UQ) centers on investigating the effects of uncertainty on free boundary propagation. Through mapping to a reference domain, we derive an enthalpy-based numerical scheme tailored to the transformed coordinate system, facilitating a simple and efficient simulation. Numerical and UQ studies in one and two dimensions validate the proposed model and the extended enthalpy method. They offer intriguing insights into ice accretion and other multiphysics processes involving phase transitions.

Non-linear effects in long-haul, high-speed optical fiber systems significantly hinder channel capacity. While the Digital Backward Propagation algorithm (DBP) with adaptive filter (ADF) can mitigate these effects, it suffers from an overwhelming computational complexity. Recent solutions have incorporated deep neural networks in a data-driven strategy to alleviate this complexity in the DBP model. However, these models are often limited to a specific symbol rate and channel number, necessitating retraining for different settings, their performance declines significantly under high-speed and high-power conditions. We introduce Meta-DSP, a novel data-driven nonlinear compensation model based on meta-learning that processes multi-modal data across diverse transmission rates, power levels, and channel numbers. This not only enhances signal quality but also substantially reduces the complexity of the nonlinear processing algorithm. Our model delivers a 0.7 dB increase in the Q-factor over Electronic Dispersion Compensation (EDC), and compared to DBP, it curtails computational complexity by a factor of ten while retaining comparable performance. From the perspective of the entire signal processing system, the core idea of Meta-DSP can be employed in any segment of the overall communication system to enhance the model's scalability and generalization performance. Our research substantiates Meta-DSP's proficiency in addressing the critical parameters defining optical communication networks.

In this paper, we study the semiclassical behavior of Bloch electrons in the presence of Weyl nodes, which are singular points in the band structure of certain materials. We carry out asymptotic analysis and present a rigorous derivation of the semiclassical asymptotic expansion of the current of Bloch electrons with the presence of Weyl nodes. The analysis shows that the current contains two parts, one independent of the Weyl nodes and the other a contribution from the singular points. This work provides a theoretical foundation towards a rigorous justification of recent scientific discoveries in Weyl semimetals. The main innovation of this paper is a new strategy to deal with the singular points with quantitative estimates, which may have broader applications in multiscale models with singularities.

Prompt Engineering (PE) has emerged as a critical technique for guiding Large Language Models (LLMs) in solving intricate tasks. Its importance is highlighted by its potential to significantly enhance the efficiency and effectiveness of human-machine interaction. As tasks grow increasingly complex, recent advanced PE methods have extended beyond the limitations of single-round interactions to embrace multi-round interactions, which allows for a deeper and more nuanced engagement with LLMs. In this paper, we propose an optimal control framework tailored for multi-round interactions with LLMs. This framework provides a unified mathematical structure that not only systematizes the existing PE methods but also sets the stage for rigorous analytical improvements. Furthermore, we extend this framework to include PE via ensemble methods and multi-agent collaboration, thereby enlarging the scope of applicability. By adopting an optimal control perspective, we offer fresh insights into existing PE methods and highlight theoretical challenges that warrant future research. Besides, our work lays a foundation for the development of more effective and interpretable PE methods.

This paper considers mean field games with optimal stopping time (OSMFGs) where agents make optimal exit decisions, the coupled obstacle and Fokker-Planck equations in such models pose challenges versus classic MFGs. This paper proposes a generalized fictitious play algorithm that computes OSMFG mixed equilibria by iteratively solving pure strategy systems, i.e. approximating mixed strategies through averaging pure strategies according to a certain updating rule. The generalized fictitious play allows for a broad family of learning rates and the convergence to the mixed strategy equilibrium can be rigorously justified. The algorithm also incorporates efficient finite difference schemes of the pure strategy system, and numerical experiments demonstrate the effectiveness of the proposed method in robustly and efficiently computing mixed equilibria for OSMFGs.

We develop a stochastic Schr\"odinger equation (SSE) framework to simulate real-time dynamics of Anderson-Holstein (AH) impurities coupled to a continuous fermionic bath. The bath degrees of freedom are incorporated through fluctuating terms determined by exact system-bath correlations, which is derived in an ab initio manner. We show that such an SSE treatment provides a middle ground between numerically expansive microscopic simulations and macroscopic master equations. Computationally, the SSE model enables efficient numerical methods for propagating stochastic trajectories. We demonstrate that this approach not only naturally provides microscopically-detailed information unavailable from reduced models, but also captures effects beyond master equations, thus serves as a promising tool to study open quantum dynamics emerging in physics and chemistry., Comment: 16 pages, 8 figures

The quantum thermal average is a central topic in quantum physics and can be represented by the path integrals. For the computational perspective, the path integral representation (PIR) needs to be approximated in a finite-dimensional space, and the convergence of such approximation is termed as the convergence of the PIR. In this paper, we establish the Trotter product formula in the trace form, which connects the quantum thermal average and the Boltzmann distribution of a continuous loop in a rigorous way. We prove the qualitative convergence of the standard PIR, and obtain the explicit convergence rates of the continuous loop PIR. These results showcase various approaches to approximate the quantum thermal average, which provide theoretical guarantee for the path integral approaches of quantum thermal equilibrium systems, such as the path integral molecular dynamics.

Motivated by the incompressible limit of a cell density model, we propose a free boundary tumor growth model where the pressure satisfies an obstacle problem on an evolving domain $\Omega(t)$, and the coincidence set $\Lambda(t)$ captures the emerging necrotic core. We contribute to the analytical characterization of the solution structure in the following two aspects. By deriving a semi-analytical solution and studying its dynamical behavior, we obtain quantitative transitional properties of the solution separating phases in the development of necrotic cores and establish its long time limit with the traveling wave solutions. Also, we prove the existence of traveling wave solutions incorporating non-zero outer densities outside the tumor bulk, provided that the size of the outer density is below a threshold.

The quantum thermal average plays a central role in describing the thermodynamic properties of a quantum system. Path integral molecular dynamics (PIMD) is a prevailing approach for computing quantum thermal averages by approximating the quantum partition function as a classical isomorphism on an augmented space, enabling efficient classical sampling, but the theoretical knowledge of the ergodicity of the sampling is lacking. Parallel to the standard PIMD with $N$ ring polymer beads, we also study the Matsubara mode PIMD, where the ring polymer is replaced by a continuous loop composed of $N$ Matsubara modes. Utilizing the generalized $\Gamma$ calculus, we prove that both the Matsubara mode PIMD and the standard PIMD have uniform-in-$N$ ergodicity, i.e., the convergence rate towards the invariant distribution does not depend on the number of modes or beads $N$.

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: 30 pages, 17 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

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., Comment: 23 pages, 3 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.

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.

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 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$.

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

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., Comment: 43 pages, 3 figures

In this article, we propose a Frozen Gaussian Sampling (FGS) algorithm for simulating nonadiabatic quantum dynamics at metal surfaces with a continuous spectrum. This method consists of a Monte-Carlo algorithm for sampling the initial wave packets on the phase space and a surface-hopping type stochastic time propagation scheme for the wave packets. We prove that to reach a certain accuracy threshold, the sample size required is independent of both the semiclassical parameter $\varepsilon$ and the number of metal orbitals $N$, which makes it one of the most promising methods to study the nonadiabatic dynamics. The algorithm and its convergence properties are also validated numerically. Furthermore, we carry out numerical experiments including exploring the nuclei dynamics, electron transfer and finite-temperature effects, and demonstrate that our method captures the physics which can not be captured by classical surface hopping trajectories., Comment: 41 pages, 10 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 $N$, the time step $\tau$ for batch divisions or the batch size $p$. Moreover, the Wasserstein distance between the invariant distributions of the IPS and the RB-IPS is bounded by $O(\sqrt{\tau})$, 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

In this paper, we develop a Monte Carlo algorithm named the Frozen Gaussian Sampling (FGS) to solve the semiclassical Schr\"odinger 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., Comment: 41 pages, 7 figures

In this work, we are concerned with a Fokker-Planck equation related to the nonlinear noisy leaky integrate-and-fire model for biological neural networks which are structured by the synaptic weights and equipped with the Hebbian learning rule. The equation contains a small parameter $\varepsilon$ separating the time scales of learning and reacting behavior of the neural system, and an asymptotic limit model can be derived by letting $\varepsilon\to 0$, where the microscopic quasi-static states and the macroscopic evolution equation are coupled through the total firing rate. To handle the endowed flux-shift structure and the multi-scale dynamics in a unified framework, we propose a numerical scheme for this equation that is mass conservative, unconditionally positivity preserving, and asymptotic preserving. We provide extensive numerical tests to verify the schemes' properties and carry out a set of numerical experiments to investigate the model's learning ability, and explore the solution's behavior when the neural network is excitatory., Comment: 24 pages, 8 figures. arXiv admin note: text overlap with arXiv:1706.05796 by other authors

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

In this article, we develop and analyse a new spectral method to solve the semi-classical Schr\"odinger 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\"odinger 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., Comment: A revised version. An appendix is added to establish estimate (3.15) and (3.20)

In the mean field integrate-and-fire model, the dynamics of a typical neuron within a large network is modeled as a diffusion-jump stochastic process whose jump takes place once the voltage reaches a threshold. In this work, the main goal is to establish the convergence relationship between the regularized process and the original one where in the regularized process, the jump mechanism is replaced by a Poisson dynamic, and jump intensity within the classically forbidden domain goes to infinity as the regularization parameter vanishes. On the macroscopic level, the Fokker-Planck equation for the process with random discharges (i.e. Poisson jumps) are defined on the whole space, while the equation for the limit process is on the half space. However, with the iteration scheme, the difficulty due to the domain differences has been greatly mitigated and the convergence for the stochastic process and the firing rates can be established. Moreover, we find a polynomial-order convergence for the distribution by a re-normalization argument in probability theory. Finally, by numerical experiments, we quantitatively explore the rate and the asymptotic behavior of the convergence for both linear and nonlinear models., Comment: This is a new version of the existed paper and I should not submit a new one. Please see arXiv:2009.04679

An efficient sampling method, the pmmLang+RBM, is proposed to compute the quantum thermal average in the interacting quantum particle system. Benefiting from the random batch method (RBM), the pmmLang+RBM reduces the complexity due to the interaction forces per timestep from $O(NP^2)$ to $O(NP)$, where $N$ is the number of beads and $P$ is the number of particles. Although the RBM introduces a random perturbation of the interaction forces at each timestep, the long time effects of the random perturbations along the sampling process only result in a small bias in the empirical measure of the pmmLang+RBM from the target distribution, which also implies a small error in the thermal average calculation. We numerically study the convergence of the pmmLang+RBM, and quantitatively investigate the dependence of the error in computing the thermal average on the parameters including the batch size, the timestep, etc. We also propose an extension of the pmmLang+RBM, which is based on the splitting Monte Carlo method and is applicable when the interacting potential contains a singular part.

In the mean field integrate-and-fire model, the dynamics of a typical neuron within a large network is modeled as a diffusion-jump stochastic process whose jump takes place once the voltage reaches a threshold. In this work, the main goal is to establish the convergence relationship between the regularized process and the original one where in the regularized process, the jump mechanism is replaced by a Poisson dynamic, and jump intensity within the classically forbidden domain goes to infinity as the regularization parameter vanishes. On the macroscopic level, the Fokker-Planck equation for the process with random discharges (i.e. Poisson jumps) are defined on the whole space, while the equation for the limit process is on the half space. However, with the iteration scheme, the difficulty due to the domain differences has been greatly mitigated and the convergence for the stochastic process and the firing rates can be established. Moreover, we find a polynomial-order convergence for the distribution by a re-normalization argument in probability theory. Finally, by numerical experiments, we quantitatively explore the rate and the asymptotic behavior of the convergence for both linear and nonlinear models.

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 a nonlinear and nonlocal parabolic model for multi-species ionic fluids and introduce a semi-implicit finite volume scheme, which is second order accurate in space, first order in time and satisfies the following properties: positivity preserving, mass conservation and energy dissipation. Besides, our scheme involves a fast algorithm on the convolution terms with singular but integrable kernels, which otherwise impedes the accuracy and efficiency of the whole scheme. Error estimates on the fast convolution algorithm are shown next. Numerous numerical tests are provided to demonstrate the properties, such as unconditional stability, order of convergence, energy dissipation and the complexity of the fast convolution algorithm. Furthermore, extensive numerical experiments are carried out to explore the modeling effects in specific examples, such as, the steric repulsion, the concentration of ions at the boundary and the blowup phenomenon of the Keller-Segel equations., Comment: 27 pages, 17 figures

With the path integral approach, the thermal average in a multi-electronic-state quantum systems can be approximated by the ring polymer representation on an extended configuration space, where the additional degrees of freedom are associated with the surface index of each bead. The primary goal of this work is to propose a more efficient sampling algorithm for the calculation of such thermal averages. We reformulate the extended ring polymer approximation according to the configurations of the surface indexes, and by introducing a proper reference measure, the reformulation is recast as a ratio of two expectations of function expansions. By quantitatively estimating the sub-estimators, and minimizing the total variance of the sampled average, we propose a multi-level Monte Carlo path integral molecular dynamics method (MLMC-PIMD) to achieve an optimal balance of computational cost and accuracy.

In this work, the primary goal is to establish rigorous connection between the Fokker-Planck equation of neural networks with its microscopic model: the diffusion-jump stochastic process that captures the mean field behavior of collections of neurons in the integrate-and-fire model. The proof is based on a novel iteration scheme: with an auxiliary random variable counting the firing events, both the density function of the stochastic process and the solution of the PDE problem admit series representations, and thus the difficulty in verifying the link between the density function and the PDE solution in each sub problem has been greatly mitigated. The iteration approach provides a generic frame in integrating the probability approach with PDE techniques, with which we prove that the density function of the diffusion-jump stochastic process is indeed the classical solution of the Fokker-Planck equation with a unique flux-shift structure.

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., Comment: 25 pages, 14 figures

In this work, we are concerned with the Fokker-Planck equations associated with the Nonlinear Noisy Leaky Integrate-and-Fire model for neuron networks. Due to the jump mechanism at the microscopic level, such Fokker-Planck equations are endowed with an unconventional structure: transporting the boundary flux to a specific interior point. While the equations exhibit diversified solutions from various numerical observations, the properties of solutions are not yet completely understood, and by far there has been no rigorous numerical analysis work concerning such models. We propose a conservative and conditionally positivity preserving scheme for these Fokker-Planck equations, and we show that in the linear case, the semi-discrete scheme satisfies the discrete relative entropy estimate, which essentially matches the only known long time asymptotic solution property. We also provide extensive numerical tests to verify the scheme properties, and carry out several sets of numerical experiments, including finite-time blowup, convergence to equilibrium and capturing time-period solutions of the variant models.

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.

In this article, we propose a novel method for sampling potential functions based on noisy observation data of a finite number of observables in quantum canonical ensembles, which leads to the accurate sampling of a wide class of test observables. The method is based on the Bayesian inversion framework, which provides a platform for analyzing the posterior distribution and naturally leads to an efficient numerical sampling algorithm. We highlight that, the stability estimate is obtained by treating the potential functions as intermediate variables in the following way: the discrepancy between two sets of observation data of training observables can bound the distance between corresponding posterior distributions of potential functions, while the latter naturally leads to a bound of the discrepancies between corresponding thermal averages of test observables. Besides, the training observables can be more flexible than finite samples of the local density function, which are mostly used in previous researches. The method also applies to the multi-level quantum systems in the non-adiabatic regime. In addition, we provide extensive numerical tests to verify the accuracy and efficiency of the proposed algorithm.

Micromagnetics simulations require accurate approximation of the magnetization dynamics described by the Landau-Lifshitz-Gilbert equation, which is nonlinear, nonlocal, and has a non-convex constraint, posing interesting challenges in developing numerical methods. In this paper, we propose two second-order semi-implicit projection methods based on the second-order backward differentiation formula and the second-order interpolation formula using the information at previous two temporal steps. Unconditional unique solvability of both methods is proved, with their second-order accuracy verified through numerical examples in both 1D and 3D. The efficiency of both methods is compared to that of another two popular methods. In addition, we test the robustness of both methods for the first benchmark problem with a ferromagnetic thin film material from National Institute of Standards and Technology.

In this work, we study the semi-classical limit of the Schr\"odinger equation with random inputs, and show that the semi-classical Schr\"odinger equation produces $O(\varepsilon)$ oscillations in the random variable space. With the Gaussian wave packet transform, the original Schr\"odinger equation is mapped to an ODE system for the wave packet parameters coupled with a PDE for the quantity $w$ in rescaled variables. Further, we show that the $w$ equation does not produce $\varepsilon$ dependent oscillations, and thus it is more amenable for numerical simulations. We propose multi-level sampling strategy in implementing the Gaussian wave packet transform, where in the most costly part, i.e. simulating the $w$ equation, it is sufficient to use $\varepsilon$ independent samples. We also provide extensive numerical tests as well as meaningful numerical experiments to justify the properties of the numerical algorithm, and hopefully shed light on possible future directions.

We investigate the continuum limit that the number of beads goes to infinity in the ring polymer representation of thermal averages. Studying the continuum limit of the trajectory sampling equation sheds light on possible preconditioning techniques for sampling ring polymer configurations with large number of beads. We propose two preconditioned Langevin sampling dynamics, which are shown to have improved stability and sampling accuracy. We present a careful mode analysis of the preconditioned dynamics and show their connections to the normal mode, the staging coordinate and the Matsubara mode representation for ring polymers. In the case where the potential is quadratic, we show that the continuum limit of the preconditioned mass modified Langevin dynamics converges to its equilibrium exponentially fast, which suggests that the finite-dimensional counterpart has a dimension-independent convergence rate. In addition, the preconditioning techniques can be naturally applied to the multi-level quantum systems in the nonadiabatic regime, which are compatible with various numerical approaches.