Large Systems Of Equations Research Articles

A Multigrid Solver for PDE-Constrained Optimization with Uncertain Inputs

In this manuscript, we present a collective multigrid algorithm to solve efficiently the large saddle-point systems of equations that typically arise in PDE-constrained optimization under uncertainty, and develop a novel convergence analysis of collective smoothers and collective two-level methods. The multigrid algorithm is based on a collective smoother that at each iteration sweeps over the nodes of the computational mesh, and solves a reduced saddle-point system whose size is proportional to the number N of samples used to discretized the probability space. We show that this reduced system can be solved with optimal O(N) complexity. The multigrid method is tested both as a stationary method and as a preconditioner for GMRES on three problems: a linear-quadratic problem, possibly with a local or a boundary control, for which the multigrid method is used to solve directly the linear optimality system; a nonsmooth problem with box constraints and L1\\documentclass[12pt]{minimal} \\usepackage{amsmath} \\usepackage{wasysym} \\usepackage{amsfonts} \\usepackage{amssymb} \\usepackage{amsbsy} \\usepackage{mathrsfs} \\usepackage{upgreek} \\setlength{\\oddsidemargin}{-69pt} \\begin{document}$$L^1$$\\end{document}-norm penalization on the control, in which the multigrid scheme is used as an inner solver within a semismooth Newton iteration; a risk-averse problem with the smoothed CVaR risk measure where the multigrid method is called within a preconditioned Newton iteration. In all cases, the multigrid algorithm exhibits excellent performances and robustness with respect to the parameters of interest.

A Method for Specifying Complete Signature Randomization and an Algebraic Algorithm Based on It

To eliminate the limitations of signature randomization in known algebraic algorithms with a hidden group, the security of which is based on the computational complexity of solving large systems of power equations, a method for ensuring complete randomization is proposed. Based on this method, a new algorithm of the indicated type was developed, using a four-dimensional finite non-commutative associative algebra as an algebraic basis. We obtained estimates of the security of algorithms to direct attacks as well as from attacks based on known signatures, which confirm the effectiveness of the proposed signature randomization method. Due to the relatively small size and signature of the public and private keys, the developed algorithm is of interest as a potential practical post-quantum digital signature scheme.

Quasi-Interpolation on Chebyshev Grids with Boundary Corrections

Quasi-interpolation is a powerful tool for approximating functions using radial basis functions (RBFs) such as Gaussian kernels. This avoids solving large systems of equations as in RBF interpolation. However, quasi-interpolation with Gaussian kernels on compact intervals can have significant errors near the boundaries. This paper proposes a quasi-interpolation method with Gaussian kernels using Chebyshev points and boundary corrections to improve the approximation near the boundaries. The boundary corrections use a linear approximation of the function beyond the interval to estimate the truncation error and add correction terms. Numerical studies on test functions show that the proposed method reduces errors significantly near boundaries compared to quasi-interpolation without corrections, for both equally spaced and Chebyshev points. The convergence and accuracy with the boundary corrections are generally better with Chebyshev points compared to equally spaced points. The proposed method provides an efficient way to perform quasi-interpolation on compact intervals while controlling the boundary errors. This study introduces a novel approach to quasi-interpolation modification, which significantly enhances convergence rates and minimizes errors at boundary points, thereby advancing the methods for boundary approximation.

Jacobian-free Locally Linearized Runge–Kutta method of Dormand and Prince for large systems of differential equations

This paper introduces Jacobian-free Locally Linearized Runge–Kutta formulas of Dormand and Prince for integrating large systems of initial value problems. With these new Jacobian-free formulas, an adaptive time-stepping variable-order scheme with regulated Krylov dimension is constructed with new developments in the computation of the phi-function action over vectors through Jacobian-free Krylov–Padé approximations, and in the procedures for controlling the approximation errors and for estimating the dimension of the Krylov subspaces. At each integration step, the implemented scheme performs only one Krylov subspace decomposition and a few computations of low dimensional exponential matrices, which contrasts with the implementation of other exponential-type integrators of high order. Regarding to existing schemes, the performed numerical study demonstrates a higher effectiveness of the constructed scheme in the integration of a variety of test equations.

Variable Time-stepping Exponential Integrators for Chemical Reactors with Analytical Jacobians

Chemical combustion problems are known to be stiff and therefore difficult to efficiently integrate in time when numerically simulated. Implicit methods, such as backwards differentiation formula (BDF), are widely considered to be the state-of-the-art methods owing their capability of taking relatively large time-steps while maintaining accurate combustion characteristics. Exponential time integration methods have recently demonstrated the ability to accurately and efficiently solve large scale systems of ordinary differential equations. This study introduces a novel adaptive time stepping exponential integrator named EPI3V. Its performance is measured on spatially homogeneous isobaric reactive mixtures involving three hydrocarbon fuel mechanisms. The full combustion process is simulated using gas compositions with sufficient temperature to obtain auto-ignition. Simulations are run until the steady state is obtained, then a comparison of the computational efficiency and accuracy between a BDF and EPI3V method is made. The novel EPI3V method exhibits comparable computational efficiency to a well-established implementation of the variable time-stepping BDF implicit methods for two of the mechanisms investigated. In certain situations it even demonstrates a slight advantage over the implicit solver. However, in one specific case, the EPI3V shows relative performance degradation compared to the implicit method, but it still converges for this case. These results indicate that exponential time integration methods may be applicable to a larger variety of combustion problems.

Simple statistical tests selection based parallel computating method ensures the guaranteed global extremum identification

The article proposes a parallel computing oriented method for solving the global minimum finding problems, in which continuous objective functions satisfy the Hölder condition, and the control parameters domain limited by continuous functions is characterized by a positive Lebesgue measure. A typical example of such a task is the discrepancy minimizing problem between the left and right parts of some large system of equilibrium equations (this is a usual situation when describing real process using Markov chain). The method is based on simple statistical tests, thanks to which, at each iteration, growing sets of potential global minima and sets of decrements necessary for estimating the values of the Hölder constants are formed. The article theoretically substantiates and empirically proves the guaranteed convergence of the authors’ method to the real global minimum, which occurs at an exponential rate. For the continuous iterations number, analytical upper estimates of the spacing between the potential global minima and real global minima are formalized, as well as an estimate of the probability of overcoming this spacing is formalized. The decrements sequence approximation, estimation of a priori unknown Hölder constants, estimation of the average number of iterations of the method and probabilistic characteristics of the final solution are analytically justified. In addition to the theoretical proof, the adequacy of the authors’ method has been confirmed empirically. It turned out that both the quality characteristics of the initial results calculated by the authors’ method and the time to obtain them are practically independent of the size of the search area. This expected result is a significant advantage of the authors’ method over analogues.

Tensor reduction of loop integrals

The computational cost associated with reducing tensor integrals to scalar integrals using the Passarino-Veltman method is dominated by the diagonalisation of large systems of equations. These systems of equations are sized according to the number of independent tensor elements that can be constructed using the metric and external momenta. In this article, we present a closed-form solution of this diagonalisation problem in arbitrary tensor integrals. We employ a basis of tensors whose building blocks are the external momentum vectors and a metric tensor transverse to the space of external momenta. The scalar integral coefficients of the basis tensors are obtained by mapping the basis elements to the elements of an orthogonal dual basis. This mapping is succinctly expressed through a formula that resembles the ordering of operators in Wick’s theorem.Finally, we provide examples demonstrating the application of our tensor reduction formula to Feynman diagrams in QCD 2 → 2 scattering processes, specifically up to three loops.

An Adaptive Fast-Multipole-Accelerated Hybrid Boundary Integral Equation Method for Accurate Diffusion Curves

In theory, diffusion curves promise complex color gradations for infinite-resolution vector graphics. In practice, existing realizations suffer from poor scaling, discretization artifacts, or insufficient support for rich boundary conditions. Previous applications of the boundary element method to diffusion curves have relied on polygonal approximations, which either forfeit the high-order smoothness of Bézier curves, or, when the polygonal approximation is extremely detailed, result in large and costly systems of equations that must be solved. In this paper, we utilize the boundary integral equation method to accurately and efficiently solve the underlying partial differential equation. Given a desired resolution and viewport, we then interpolate this solution and use the boundary element method to render it. We couple this hybrid approach with the fast multipole method on a non-uniform quadtree for efficient computation. Furthermore, we introduce an adaptive strategy to enable truly scalable infinite-resolution diffusion curves.

Reproducibility of biophysical in silico neuron states and spikes from event-based partial histories.

Biophysically detailed simulations of neuronal activity often rely on solving large systems of differential equations; in some models, these systems have tens of thousands of states per cell. Numerically solving these equations is computationally intensive and requires making assumptions about the initial cell states. Additional realism from incorporating more biological detail is achieved at the cost of increasingly more states, more computational resources, and more modeling assumptions. We show that for both a point and morphologically-detailed cell model, the presence and timing of future action potentials is probabilistically well-characterized by the relative timings of a moderate number of recent events alone. Knowledge of initial conditions or full synaptic input history is not required. While model time constants, etc. impact the specifics, we demonstrate that for both individual spikes and sustained cellular activity, the uncertainty in spike response decreases as the number of known input events increases, to the point of approximate determinism. Further, we show cellular model states are reconstructable from ongoing synaptic events, despite unknown initial conditions. We propose that a strictly event-based modeling framework is capable of representing the complexity of cellular dynamics of the differential-equations models with significantly less per-cell state variables, thus offering a pathway toward utilizing modern data-driven modeling to scale up to larger network models while preserving individual cellular biophysics.

Planning ride-pooling services with detour restrictions for spatially heterogeneous demand: A multi-zone queuing network approach

This study presents a multi-zone queuing network model for steady-state ride-pooling operations that serve heterogeneous demand, and then builds upon this model to optimize the design of ride-pooling services. Spatial heterogeneity is addressed by partitioning the study region into a set of relatively homogeneous zones, and a set of criteria are imposed to avoid significant detours among matched passengers. A generalized multi-zone queuing network model is then developed to describe how vehicles’ states transition within each zone and across neighboring zones, and how passengers are served by idle or partially occupied vehicles. A large system of equations is constructed based on the queuing network model to analytically evaluate steady-state system performance. Then, we formulate a constrained nonlinear program to optimize the design of ride-pooling services, such as zone-level vehicle deployment, vehicle routing paths, and vehicle rebalancing operations. A customized solution approach is also proposed to decompose and solve the optimization problem. The proposed model and solution approach are applied to a hypothetical case and a real-world Chicago case study, so as to demonstrate their applicability and to draw insights. Agent-based simulations are also used to corroborate results from the proposed analytical model. These numerical examples not only reveal interesting insights on how ride-pooling services serve heterogeneous demand, but also highlight the importance of addressing demand heterogeneity when designing ride-pooling services.

A very high order Flux Reconstruction (FR) method for the numerical simulation of 1-D compositional fluid flow model in petroleum reservoirs

Compositional reservoir simulation is an extremely important tool for modeling the fluid flow in complex petroleum reservoirs, as in cases where the reservoir fluid is volatile and composed of several pseudo-components with distinct characteristics, or reservoirs that require the use of enhanced oil recovery process. For these cases, simpler black-oil models are not suitable. The compositional model involves the solution of a large system of partial differential equations, that comprises the mass conservation, the Darcy’s law and the fugacity constraints. However, the high number of equations and constraints render the compositional problems extremely complex and computationally demanding solutions. To improve accuracy and reduce the high computational costs, one can use higher than first or even second-order methods to approximate the advective flux terms in the hyperbolic conservation laws that describe the multicomponent transport in the reservoir. Those very high-order schemes can replace low-order approaches, such as the First Order Upwind (FOU) method, which is traditionally used in commercial petroleum reservoir simulators (CMG, 2019; ECLIPSE, 2022), obtaining more accurate solutions with reduced computational cost. In this work, we consider a three-phase and isothermal fluid flow of water, oil and gas and assume that there is no mass transfer between the water and the hydrocarbon phases. Physical dispersion and capillary effects are neglected. To solve the system of Partial Differential Equations (PDEs), we used the classical IMPEC (Implicit Pressure Explicit Composition) approach and, to model the complex phase behavior, we used the Peng–Robinson equation of state (PR-EOS). The classical Two Point Flux Approximation (TPFA) method is used to spatially discretize the diffusion terms in the pressure equation. For the first time in literature the very high order Flux Reconstruction (FR) method is adapted for the numerical simulation of the compositional model in petroleum reservoirs. The FR method is a numerical scheme employed to discretize the hyperbolic conservation laws using general grids. To avoid spurious oscillations in the vicinity of solution discontinuities, the h-MLP (hierarchical Multi-dimensional Limiting Process) is applied in the reconstruction stage. For the time integration, we have used the third-order Runge–Kutta approach. To appraise the robustness of our formulation, we have solved some one-dimensional benchmark problems found in literature and compared the FR solution with the FOU method and with the second-order Monotonic Upstream-Centered Scheme for Conservation Laws (MUSCL) method. From our results, the FR method has shown improved accuracy, efficiency and capability of detecting sharper shocks when compared to the other methods that we have considered for the numerical modeling of the compositional flow in the present paper.

Equilibrium and surviving species in a large Lotka-Volterra system of differential equations.

Lotka-Volterra (LV) equations play a key role in the mathematical modeling of various ecological, biological and chemical systems. When the number of species (or, depending on the viewpoint, chemical components) becomes large, basic but fundamental questions such as computing the number of surviving species still lack theoretical answers. In this paper, we consider a large system of LV equations where the interactions between the various species are a realization of a random matrix. We provide conditions to have a unique equilibrium and present a heuristics to compute the number of surviving species. This heuristics combines arguments from Random Matrix Theory, mathematical optimization (LCP), and standard extreme value theory. Numerical simulations, together with an empirical study where the strength of interactions evolves with time, illustrate the accuracy and scope of the results.

A solution of the Crow-Kimura evolution model on fluctuating fitness landscape

The article discusses the Crow-Kimura model in the context of random transitions between different fitness landscapes. The duration of epochs, during which the fitness landscape is constant over time, is modeled by an exponential distribution. To obtain an exact solution, a system of functional equations is required. However, to approximate the model, we consider the cases of slow or fast transitions and calculate the first-order corrections using either the transition rate or its inverse. Specifically, we focus on the case of slow transitions and find that the average fitness is equal to the average fitness for evolution on static fitness landscapes, but with the addition of a load term. We also investigate the model for a small number of genes and identify the exact transition points to the transient phase.

Построение С1-гладких кусочно-квадратичных функций при решении краевых задач уравнений 4-го порядка на треугольной сетке

In this paper, we present one approach to constructing continuously differentiable piecewise-quadratic functions on a triangular grid, based on smoothing a piecewise-linear function in the vicinity of edges and triangulation nodes. In the works [8], [7], the issues of approximation of the functional (1) in triangular grids and the convergence of the variational method for solving the boundary value problem of the equation (2) were studied. However, in the numerical solution there are difficulties associated with the construction of continuously differentiable piecewise polynomial functions on triangulations. In particular, their construction requires solving large systems of equations at each step of the variational method. When trying to get by with only continuous piecewise polynomial functions, we got a negative result (divergence of approximate solutions was found) [9]. In this paper, we circumvent the difficulties that arise — we indicate a method for constructing piecewise polynomial functions that have continuous partial derivatives. With the help of this class of functions, we have obtained formulas for approximating the functional (1) and tested the method on the example of a biharmonic equation.

Multi-objective optimization through a novel Bayesian approach for industrial manufacturing of Polyvinyl Acetate

ABSTRACT Manufacturing long-chain branched polymers enhances the quality of products at the cost of production time, translating into high manufacturing costs. In this manuscript, we study the cost-versus-quality trade-off while optimizing the operating conditions for large-scale industrial production of Polyvinyl acetate (PVAc). PVAc polymerization is emulated using a large system of stiff differential equations resulting in time-expensive function evaluation. Thus, we aim to achieve the Pareto optimality using Multi-Objective Bayesian Optimization (MOBO) that introduces q-Expected Hyper-Volume Improvement (qEHVI) as the novel acquisition function. Comparison with Non-dominated Sorting Genetic Algorithm-II (NSGA-II), applied to the high-fidelity model, resulted in a similar Pareto front with only 1% of the high-fidelity calls required by NSGA-II. The results indicate the efficiency of MOBO for PVAc optimization and present a real-world application of a generic method that can be implemented to solve time-expensive multi-objective optimization in manufacturing processes

A new, fast method for solving finite-element equations iteratively based on Gauss–Seidel

Solving large equation systems is the most time-consuming part of finite-element modelling – iterative techniques are favoured for models with numerous degrees of freedom where direct techniques have high storage requirements. Classical iterative techniques such as Gauss–Seidel (GS) are robust due to guaranteed convergence and algorithmic simplicity. Performance of iterative techniques chiefly depends on system scale and stiffness matrix properties – which are influenced by structural configuration. However, it is possible to adjust an iterative algorithm such that its speed is greatly enhanced for a certain class of structural configurations. This paper presents an adjusted GS solver, ‘constrained Gauss–Seidel’ (CGS), which has been formulated to solve typical multi-storey structures with an enhanced speed. The innovation in CGS originates from the adoption of a diaphragmatic relaxation mechanism that results in dividing equations into two groups to optimally deal with the different unknown types. In this paper, the concept and algorithm of the newly developed CGS method are elucidated. Then, 16 practical examples are solved to assess the solving speed of CGS against other iterative methods – GS and modified Gauss–Seidel (MGS, MGS*). The convergence speed of CGS attained 33 times, 3.7 times and 2 times those of GS, MGS and MGS*, respectively.

Analysis of Magnetohydrodynamic Micropolar Nanofluid Flow due to Radially Stretchable Rotating Disk Employing Spectral Method

The present analysis is aimed at examining MHD micropolar nanofluid flow past a radially stretchable rotating disk with the Cattaneo-Christov non-Fourier heat and non-Fick mass flux model. To begin with, the model is developed in the form of nonlinear partial differential equations (PDEs) for momentum, microrotation, thermal, and concentration with their boundary conditions. Employing suitable similarity transformation, the boundary layer micropolar nanofluid flows governing these PDEs are transformed into large systems of dimensionless coupled nonlinear ordinary differential equations (ODEs). These dimensionless ODEs are solved numerically by means of the spectral local linearization method (SLLM). The consequences of more noticeable involved parameters on different flow fields and engineering quantities of interest are thoroughly inspected, and the results are presented via graph plots and tables. The obtained results confirm that SLLM is a stable, accurate, convergent, and computationally very efficient method to solve a large coupled system of equations. The radial velocity grows while the tangential velocity, temperature, and concentration distributions turn down as the value of the radial stretching parameter improves, and hence, in practical applications, radial stretching of the disk is helpful to advance the cooling process of the rotating disk. The occurrence of microrotation viscosity in microrotation parameters ( A 1 − A 6 ) declines the radial velocity profile, and the kinetic energy of the fluid is reduced to some extent far away from the surface of the disk. The novelty of the study is the consideration of microscopic effects occurring from the micropolar fluid elements such as micromotion and couple stress, the effects of non-Fourier’s heat and non-Fick’s mass flux, and the effect of radial stretching disk on micropolar nanofluid flow, heat, and mass transfer.

Multiscale solver for multi-component reaction–diffusion systems in heterogeneous media

Coupled nonlinear system of reaction–diffusion equations describing multi-component (species) interactions with heterogeneous coefficients is considered. Finite volume method based approximation for the space is used to construct semi-discrete form for the computation of numerical solutions. Two techniques for time approximations, namely, a fully implicit (FI) and a semi-implicit (SI) schemes are examined. The fully implicit scheme is constructed using Newton’s method and leads to the coupled system of equations on each nonlinear and time iterations which is computationally rather expensive. In order to minimize the latter hurdle, an efficient and fast multiscale solver is proposed for reaction–diffusion systems in heterogeneous media. To construct fast solver, we apply a semi-implicit scheme that leads to an uncoupled system for each individual component. Problems in heterogeneous domains require a very fine grid for accurate solutions of large systems of equations at each time step iteration. Here we present a multiscale model reduction technique to reduce the size of the discrete system. Multiscale solver is based on the uncoupled operator of the problem and constructed by the use of Generalized Multiscale Finite Element Method (GMsFEM). In GMsFEM we use a diffusion part of the operator and construct multiscale basis functions by solving spectral problems in each local domain associated with the coarse grid nodes. We collect multiscale basis functions to construct a projection/prolongation matrix and generate reduced order model on the coarse grid for fast solution. Moreover, the prolongation operator is used to reconstruct a fine-scale solution and accurate approximation of the reaction part of the problem which then leads to a very accurate and computationally effective multiscale solver. We provide numerical results for two species competition test problems in two-dimensional domain with heterogeneous inclusions. Our computed solutions predict strong competition between the two species and demonstrate that diffusion can dictate the dominance of one species over the other in heterogeneous environments. We investigate the influence of number of the multiscale basis functions to the method accuracy and ability to work with different values of the diffusion coefficients.

Quantum algorithms for geologic fracture networks

Solving large systems of equations is a challenge for modeling natural phenomena, such as simulating subsurface flow. To avoid systems that are intractable on current computers, it is often necessary to neglect information at small scales, an approach known as coarse-graining. For many practical applications, such as flow in porous, homogenous materials, coarse-graining offers a sufficiently-accurate approximation of the solution. Unfortunately, fractured systems cannot be accurately coarse-grained, as critical network topology exists at the smallest scales, including topology that can push the network across a percolation threshold. Therefore, new techniques are necessary to accurately model important fracture systems. Quantum algorithms for solving linear systems offer a theoretically-exponential improvement over their classical counterparts, and in this work we introduce two quantum algorithms for fractured flow. The first algorithm, designed for future quantum computers which operate without error, has enormous potential, but we demonstrate that current hardware is too noisy for adequate performance. The second algorithm, designed to be noise resilient, already performs well for problems of small to medium size (order 10–1000 nodes), which we demonstrate experimentally and explain theoretically. We expect further improvements by leveraging quantum error mitigation and preconditioning.

Modeling of the Activators Regenerated by Electron Transfer Copolymerization of Styrene-Acrylonitrile with Prediction of the Bivariate Molecular Weight Distribution─Copolymer Composition Distribution Using Parallel Computing and Probability Generating Functions

A mathematical model of the activators regenerated by electron transfer atom transfer radical polymerization copolymerization of styrene-acrylonitrile (SAN) is developed for the first time that calculates the bivariate molecular weight distribution (MWD)-copolymer composition distribution (CCD) of the copolymer. In addition, the model calculates the overall MWD, the overall CCD, average molecular weights and composition, and conversion. The model is implemented in the open-source programming language Julia, a high-level language designed for scientific computing, which is easy to use and that performs very fast calculations. The probability generating function technique is used to process the infinite population balances of the system. The mathematical model, composed of a very large system of equations, is solved by applying parallel computing to speed up its execution. In this way, the distributions of polymer properties are calculated in a very efficient manner. The evolution of the joint MWD–CCD and other molecular properties along the reaction path is analyzed under different operating conditions. The model allows obtaining very detailed insights into the copolymer microstructure during the polymerization reaction.