Search

Iterative methods such as Jacobi, Gauss-Seidel, and Successive Over-Relaxation (SOR) are fundamental tools in solving large systems of linear equations across various scientific fields, particularly in the field of data science which has become increasingly relevant in the past decade. Iterative methods' use of matrix multiplication rather than matrix inverses makes them ideal for solving large systems quickly. Our research explores the factors of each method that define their respective strengths, limitations, and convergence behaviors to understand how these methods address drawbacks encountered when performing matrix operations by hand, as well as how they can be used in real world applications. After implementing each method by hand to understand how the algorithms work, we developed a Python program that assesses a user-given matrix based on each method's specific convergence criteria. The program compares the spectral radii of all three methods and chooses to execute whichever will yield the fastest convergence rate. Our research revealed the importance of mathematical modeling and understanding specific properties of the coefficient matrix. We observed that Gauss-Seidel is usually the most efficient method because it is faster than Jacobi and doesn't have as strict requirements as SOR, however SOR is ideal in terms of computation speed. We applied the knowledge we gained to create a traffic flow model of the I-485 highway in Charlotte. After creating a program that generates the matrix for this model, we were able to iteratively approximate the flow of cars through neighboring exits using data from the N.C. Department of Transportation. This information identifies which areas are the most congested and can be used to inform future infrastructure development.

Applications such as unbalanced and fully shuffled regression can be approached by optimizing regularized optimal transport (OT) distances, such as the entropic OT and Sinkhorn distances. A common approach for this optimization is to use a first-order optimizer, which requires the gradient of the OT distance. For faster convergence, one might also resort to a second-order optimizer, which additionally requires the Hessian. The computations of these derivatives are crucial for efficient and accurate optimization. However, they present significant challenges in terms of memory consumption and numerical instability, especially for large datasets and small regularization strengths. We circumvent these issues by analytically computing the gradients for OT distances and the Hessian for the entropic OT distance, which was not previously used due to intricate tensor-wise calculations and the complex dependency on parameters within the bi-level loss function. Through analytical derivation and spectral analysis, we identify and resolve the numerical instability caused by the singularity and ill-posedness of a key linear system. Consequently, we achieve scalable and stable computation of the Hessian, enabling the implementation of the stochastic gradient descent (SGD)-Newton methods. Tests on shuffled regression examples demonstrate that the second stage of the SGD-Newton method converges orders of magnitude faster than the gradient descent-only method while achieving significantly more accurate parameter estimations.

Within the likes of any highly contagious and unpredictable disease, lies a predictable and attainable growth rate that researchers can find in order to make logistical conclusions about that particular disease and its affected regions' counterparts. The foundation that researchers pull from when studying a particular disease and looking for its growth rate is the Susceptible-Infected-Removed (SIR) model, presented by a series of differential equations. The issue with the SIR model lies not in its complexity, but actually its simplicity and lack of a potentially high-finite amount of factors; the limit being bounded by the amount of data available for that particular factor. Our research involves the application of multiple regressions to pinpoint and identify our Covid lockdown periods, followed by the modification of the SIR model. This involved creating new model approximations such as the time-delayed SIR model and the reinfected SIR model in order to take into account factors such as incubation and reinfection, and get the lowest error discrepancy as possible for our infection rate. We were able to conclude that the more factors that we took into account, our error rate became lower and our results became more accurate. We could also identify outlier Metros and draw certain conclusions on performance level and the reasons behind them. We then moved on to find correlations, if any, between the infection rates and outside factors. We looked at demographic and weather data to demonstrate whether correlations appeared. We found that there are a few factors with high correlations, including graduate education and low temperatures.

We introduce a novel approach for decomposing and learning every scale of a given multiscale objective function in $\mathbb{R}^d$, where $d\ge 1$. This approach leverages a recently demonstrated implicit bias of the optimization method of gradient descent by Kong and Tao, which enables the automatic generation of data that nearly follow Gibbs distribution with an effective potential at any desired scale. One application of this automated effective potential modeling is to construct reduced-order models. For instance, a deterministic surrogate Hamiltonian model can be developed to substantially soften the stiffness that bottlenecks the simulation, while maintaining the accuracy of phase portraits at the scale of interest. Similarly, a stochastic surrogate model can be constructed at a desired scale, such that both its equilibrium and out-of-equilibrium behaviors (characterized by auto-correlation function and mean path) align with those of a damped mechanical system with the original multiscale function being its potential. The robustness and efficiency of our proposed approach in multi-dimensional scenarios have been demonstrated through a series of numerical experiments. A by-product of our development is a method for anisotropic noise estimation and calibration. More precisely, Langevin model of stochastic mechanical systems may not have isotropic noise in practice, and we provide a systematic algorithm to quantify its covariance matrix without directly measuring the noise. In this case, the system may not admit closed form expression of its invariant distribution either, but with this tool, we can design friction matrix appropriately to calibrate the system so that its invariant distribution has a closed form expression of Gibbs., Comment: 10 Figures

Inspired by the blending method developed by [P. Seleson, S. Beneddine, and S. Prudhome, \emph{A Force-Based Coupling Scheme for Peridynamics and Classical Elasticity}, (2013)] for the nonlocal-to-local coupling, we create a symmetric and consistent blended force-based Atomistic-to-Continuum (a/c) scheme for the atomistic chain in one-dimensional space. The conditions for the well-posedness of the underlying model are established by analyzing an optimal blending size and blending type to ensure the $H^1$ semi-norm stability for the blended force-based operator. We present several numerical experiments to test and confirm the theoretical findings.

Coarse-graining or model reduction is a term describing a range of approaches used to extend the time-scale of molecular simulations by reducing the number of degrees of freedom. In the context of molecular simulation, standard coarse-graining approaches approximate the potential of mean force and use this to drive an effective Markovian model. To gain insight into this process, the simple case of a quadratic energy is studied in an overdamped setting. A hierarchy of reduced models is derived and analysed, and the merits of these different coarse-graining approaches are discussed. In particular, while standard recipes for model reduction accurately capture static equilibrium statistics, it is shown that dynamical statistics such as the mean-squared displacement display systematic error, even when a system exhibits large time-scale separation. In the linear setting studied, it is demonstrated both analytically and numerically that such models can be augmented in a simple way to better capture dynamical statistics., Comment: 33 pages, 8 figures. Version accepted for publication in SIAM Multiscale Modeling and Simulation

Large time-stepping is important for efficient long-time simulations of deterministic and stochastic Hamiltonian dynamical systems. Conventional structure-preserving integrators, while being successful for generic systems, have limited tolerance to time step size due to stability and accuracy constraints. We propose to use data to innovate classical integrators so that they can be adaptive to large time-stepping and are tailored to each specific system. In particular, we introduce NySALT, Nystr\"{o}m-type inference-based schemes adaptive to large time-stepping. The NySALT has optimal parameters for each time step learnt from data by minimizing the one-step prediction error. Thus, it is tailored for each time step size and the specific system to achieve optimal performance and tolerate large time-stepping in an adaptive fashion. We prove and numerically verify the convergence of the estimators as data size increases. Furthermore, analysis and numerical tests on the deterministic and stochastic Fermi-Pasta-Ulam (FPU) models show that NySALT enlarges the maximal admissible step size of linear stability, and quadruples the time step size of the St\"{o}rmer--Verlet and the BAOAB when maintaining similar levels of accuracy., Comment: 26 pages, 7 figures

In this work we aim to develop a unified mathematical framework and a reliable computational approach to model the brittle fracture in heterogeneous materials with variability in material microstructures, and to provide statistic metrics for quantities of interest, such as the fracture toughness. To depict the material responses and naturally describe the nucleation and growth of fractures, we consider the peridynamics model. In particular, a stochastic state-based peridynamic model is developed, where the micromechanical parameters are modeled by a finite-dimensional random vector, or a combination of random variables truncating the Karhunen-Lo\`{e}ve decomposition or the principle component analysis (PCA). To solve this stochastic peridynamic problem, probabilistic collocation method (PCM) is employed to sample the random field representing the micromechanical parameters. For each sample, the deterministic peridynamic problem is discretized with an optimization-based meshfree quadrature rule. We present rigorous analysis for the proposed scheme and demonstrate its convergence for a number of benchmark problems, showing that it sustains the asymptotic compatibility spatially and achieves an algebraic or sub-exponential convergence rate in the random space as the number of collocation points grows. Finally, to validate the applicability of this approach on real-world fracture problems, we consider the problem of crystallization toughening in glass-ceramic materials, in which the material at the microstructural scale contains both amorphous glass and crystalline phases. The proposed stochastic peridynamic solver is employed to capture the crack initiation and growth for glass-ceramics with different crystal volume fractions, and the averaged fracture toughness are calculated. The numerical estimates of fracture toughness show good consistency with experimental measurements.

In this paper we present an asymptotically compatible meshfree method for solving nonlocal equations with random coefficients, describing diffusion in heterogeneous media. In particular, the random diffusivity coefficient is described by a finite-dimensional random variable or a truncated combination of random variables with the Karhunen-Lo\`{e}ve decomposition, then a probabilistic collocation method (PCM) with sparse grids is employed to sample the stochastic process. On each sample, the deterministic nonlocal diffusion problem is discretized with an optimization-based meshfree quadrature rule. We present rigorous analysis for the proposed scheme and demonstrate convergence for a number of benchmark problems, showing that it sustains the asymptotic compatibility spatially and achieves an algebraic or sub-exponential convergence rate in the random coefficients space as the number of collocation points grows. Finally, to validate the applicability of this approach we consider a randomly heterogeneous nonlocal problem with a given spatial correlation structure, demonstrating that the proposed PCM approach achieves substantial speed-up compared to conventional Monte Carlo simulations.

In a recent paper (see [7]), a quasi-nonlocal coupling method was introduced to seamlessly bridge a nonlocal diffusion model with the classical local diffusion counterpart in a one-dimensional space. The proposed coupling framework removes inconsistency on the interface, preserves the balance of fluxes, and satisfies the maximum principle of diffusion problem. However, the numerical scheme proposed in that paper does not maintain all of these properties on a discrete level. In this paper we resolve these issues by proposing a new finite difference scheme that ensures the balance of fluxes and the discrete maximum principle. We rigorously prove these results and provide the stability and convergence analyses accordingly. In addition, we provide the Courant-Friedrichs-Lewy (CFL) condition for the new scheme and test a series of benchmark examples which confirm the theoretical findings., Comment: 24 pages, 5 figures

Efficient simulation of SDEs is essential in many applications, particularly for ergodic systems that demand efficient simulation of both short-time dynamics and large-time statistics. However, locally Lipschitz SDEs often require special treatments such as implicit schemes with small time-steps to accurately simulate the ergodic measure. We introduce a framework to construct inference-based schemes adaptive to large time-steps (ISALT) from data, achieving a reduction in time by several orders of magnitudes. The key is the statistical learning of an approximation to the infinite-dimensional discrete-time flow map. We explore the use of numerical schemes (such as the Euler-Maruyama, a hybrid RK4, and an implicit scheme) to derive informed basis functions, leading to a parameter inference problem. We introduce a scalable algorithm to estimate the parameters by least squares, and we prove the convergence of the estimators as data size increases. We test the ISALT on three non-globally Lipschitz SDEs: the 1D double-well potential, a 2D multi-scale gradient system, and the 3D stochastic Lorenz equation with degenerate noise. Numerical results show that ISALT can tolerate time-step magnitudes larger than plain numerical schemes. It reaches optimal accuracy in reproducing the invariant measure when the time-step is medium-large., Comment: 20 pages, 9 figures

Local-to-Nonlocal (LtN) coupling refers to a class of methods aimed at combining nonlocal and local modeling descriptions of a given system into a unified coupled representation. This allows to consolidate the accuracy of nonlocal models with the computational expediency of their local counterparts, while often simultaneously removing additional nonlocal modeling issues such as surface effects. The number and variety of proposed LtN coupling approaches have significantly grown in recent year, yet the field of LtN coupling continues to grow and still has open challenges. This review provides an overview of the state-of-the-art of LtN coupling in the context of nonlocal diffusion and nonlocal mechanics, specifically peridynamics. We present a classification of LtN coupling methods and discuss common features and challenges. The goal of this review is not to provide a preferred way to address LtN coupling but to present a broad perspective of the field, which would serve as guidance for practitioners in the selection of appropriate LtN coupling methods based on the characteristics and needs of the problem under consideration.

The Mori-Zwanzig formalism is applied to derive an equation for the evolution of linear observables of the overdamped Langevin equation. To illustrate the resulting equation and its use in deriving approximate models, a particular benchmark example is studied both numerically and via a formal asymptotic expansion. The example considered demonstrates the important of memory effects in determining the correct temporal behaviour of such systems., Comment: 19 pages, 5 figures

In this paper, we extend the idea of "geometric reconstruction" to couple a nonlocal diffusion model directly with the classical local diffusion in one dimensional space. This new coupling framework removes interfacial inconsistency, ensures the flux balance, and satisfies energy conservation as well as the maximum principle, whereas none of existing coupling methods for nonlocal-to-local coupling satisfies all of these properties. We establish the well-posedness and provide the stability analysis of the coupling method. We investigate the difference to the local limiting problem in terms of the nonlocal interaction range. Furthermore, we propose a first order finite difference numerical discretization and perform several numerical tests to confirm the theoretical findings. In particular, we show that the resulting numerical result is free of artifacts near the boundary of the domain where a classical local boundary condition is used, together with a coupled fully nonlocal model in the interior of the domain.

We formulate the blended force-based quasicontinuum (BQCF) method for multilattices and develop rigorous error estimates in terms of the approximation parameters: atomistic region, blending region and continuum finite element mesh. Balancing the approximation parameters yields a convergent atomistic/continuum multiscale method for multilattices with point defects, including a rigorous convergence rate in terms of the computational cost. The analysis is illustrated with numerical results for a Stone--Wales defect in graphene., Comment: Fixed typos and added additional details in section 3.3 and explanations in sections 4.5 and 4.6

We developed a new self-adjoint, consistent, and stable coupling strategy for nonlocal diffusion models, inspired by the quasinonlocal atomistic-to-continuum method for crystalline solids. The proposed coupling model is coercive with respect to the energy norms induced by the nonlocal diffusion kernels as well as the $L^2$ norm, and it satisfies the maximum principle. A finite difference approximation is used to discretize the coupled system, which inherits the property from the continuous formulation. Furthermore, we design a numerical example which shows the discrepancy between the fully nonlocal and fully local diffusions, whereas the result of the coupled diffusion agrees with that of the fully nonlocal diffusion., Comment: 28 pages, 3 figures, ams.org

We introduce a general framework for the mean-field analysis of large-scale load-balancing networks with general service distributions. Specifically, we consider a parallel server network that consists of N queues and operates under the $SQ(d)$ load balancing policy, wherein jobs have independent and identical service requirements and each incoming job is routed on arrival to the shortest of $d$ queues that are sampled uniformly at random from $N$ queues. We introduce a novel state representation and, for a large class of arrival processes, including renewal and time-inhomogeneous Poisson arrivals, and mild assumptions on the service distribution, show that the mean-field limit, as $N \rightarrow \infty$, of the state can be characterized as the unique solution of a sequence of coupled partial integro-differential equations, which we refer to as the hydrodynamic PDE. We use a numerical scheme to solve the PDE to obtain approximations to the dynamics of large networks and demonstrate the efficacy of these approximations using Monte Carlo simulations. We also illustrate how the PDE can be used to gain insight into network performance., Comment: Minor change, adding acknowledgments

We demonstrate the utility of the equation free methodology developed by one of the authors (I.G.K) for the study of scalar conservation laws with disordered initial conditions. The numerical scheme is benchmarked on exact solutions in Burgers turbulence corresponding to Levy process initial data. For these initial data, the kinetics of shock clustering is described by Smoluchowski's coagulation equation with additive kernel. The equation free methodology is used to develop a particle scheme that computes self-similar solutions to the coagulation equation, including those with fat tails.

We present a comprehensive error analysis of two prototypical atomistic-to-continuum coupling methods of blending type: the energy-based and the force-based quasicontinuum methods. Our results are valid in two and three dimensions, for finite range many-body interactions (e.g., EAM type), and in the presence of lattice defects (we consider point defects and dislocations). The two key ingredients in the analysis are (i) new force and energy consistency error estimates; and (ii) a new technique for proving energy norm stability of a/c couplings that requires only the assumption that the exact atomistic solution is a stable equilibrium.

The Dyson Brownian Motion (DBM) describes the stochastic evolution of $N$ points on the line driven by an applied potential, a Coulombic repulsion and identical, independent Brownian forcing at each point. We use an explicit tamed Euler scheme to numerically solve the Dyson Brownian motion and sample the equilibrium measure for non-quadratic potentials. The Coulomb repulsion is too singular for the SDE to satisfy the hypotheses of rigorous convergence proofs for tamed Euler schemes. Nevertheless, in practice the scheme is observed to be stable for time steps of O(1/N^2) and to relax exponentially fast to the equilibrium measure with a rate constant of O(1) independent of $N$. Further, this convergence rate appears to improve with $N$ in accordance with O(1/N) relaxation of local statistics of the Dyson Brownian motion. This allows us to use the Dyson Brownian motion to sample $N\times N$ Hermitian matrices from the invariant ensembles. The computational cost of generating $M$ independent samples is O(MN^4) with a naive scheme, and O(MN^3\log N) when a fast multipole method is used to evaluate the Coulomb interaction., Comment: 11 pages; 6 Figures