Stable Solver Research Articles

Sparse Resultant-Based Minimal Solvers in Computer Vision and Their Connection with the Action Matrix

Many computer vision applications require robust and efficient estimation of camera geometry from a minimal number of input data measurements. Minimal problems are usually formulated as complex systems of sparse polynomial equations. The systems usually are overdetermined and consist of polynomials with algebraically constrained coefficients. Most state-of-the-art efficient polynomial solvers are based on the action matrix method that has been automated and highly optimized in recent years. On the other hand, the alternative theory of sparse resultants based on the Newton polytopes has not been used so often for generating efficient solvers, primarily because the polytopes do not respect the constraints amongst the coefficients. In an attempt to tackle this challenge, here we propose a simple iterative scheme to test various subsets of the Newton polytopes and search for the most efficient solver. Moreover, we propose to use an extra polynomial with a special form to further improve the solver efficiency via Schur complement computation. We show that for some camera geometry problems our resultant-based method leads to smaller and more stable solvers than the state-of-the-art Gröbner basis-based solvers, while being significantly smaller than the state-of-the-art resultant-based methods. The proposed method can be fully automated and incorporated into existing tools for the automatic generation of efficient polynomial solvers. It provides a competitive alternative to popular Gröbner basis-based methods for minimal problems in computer vision. Additionally, we study the conditions under which the minimal solvers generated by the state-of-the-art action matrix-based methods and the proposed extra polynomial resultant-based method, are equivalent. Specifically, we consider a step-by-step comparison between the approaches based on the action matrix and the sparse resultant, followed by a set of substitutions, which would lead to equivalent minimal solvers.

Lattice real-time simulations with learned optimal kernels

We present a simulation strategy for the real-time dynamics of quantum fields, inspired by reinforcement learning. It builds on the complex Langevin approach, which it amends with system-specific prior information, a necessary prerequisite to overcome this exceptionally severe sign problem. The optimization process underlying our machine-learning approach is made possible by deploying inherently stable solvers of the complex Langevin stochastic process and a novel optimality criterion derived from insight into so-called boundary terms. This conceptual and technical progress allows us to both significantly extend the range of real-time simulations in 1+1d scalar field theory beyond the state of the art and to avoid discretization artifacts that plagued previous real-time field theory simulations. Limitations of and promising future directions are discussed. Published by the American Physical Society 2024

Deflated domain decomposition method for structural problems

The paper presents a fast and stable solver algorithm for structural problems. The point is the distance between the eigenvector of the constrained stiffness matrix and the unconstrained matrix. The coarse motions are close to the kernel of the unconstrained matrix. We use lower-frequency deformation modes to construct an iterative solver algorithm through domain decomposition expressing near-rigid-body motions, deflation algorithms, and two-level algorithms. We remove the coarse space from the solution space and hand over the iteration space to the fine space. Our solver is parallelized, and the solver thus has two sets of domain decomposition. One decomposition generates the coarse space, and the other is for parallelization. The basic framework of the solver is the parallel conjugate gradient (CG) method on the fine space. We use the CG method for the basic framework instead of the (simplest) domain decomposition method. We conducted benchmark tests using elastic static analysis for thin plate models. A comparison with the standard CG solver results shows the new solver’s high-speed performance and remarkable stability.

On a sparse and stable solver on graded meshes for solving high-dimensional parabolic pricing PDEs

The purpose of this work is to investigate the multi-dimensional Black-Scholes partial differential equation with variable coefficients numerically. The problem is of practical importance due to option pricing at the presence of multi assets. Since by increasing the dimension, the curse of dimensionality restricts the computations, the proposed solver will be constructed based on sparse arrays. Toward this goal, fourth-order finite difference approximations on graded meshes are introduced and then employed through semi-discretization. Then a sixth-order Runge-Kutta solver is employed for finding the resolution of the derived set of ordinary differential equations. The stability of the proposed scheme is furnished in detail as well. Numerical testings are given to uphold the accuracy and efficacy of the proposed procedure.

A class of linearization-based collocation methods for initial value and boundary value engineering problems

Nonlinear initial value and boundary value problems are generally classified as two different categories of problems and the methods utilized to solve them are also different. This paper proposed a uniform linearization-collocation framework through which the initial value and boundary value problems can be solved in a similar way. An optional pre-conditioning procedure is provided to guarantee the computational accuracy. Rigorous analysis about the convergence rate is given. Different from the commonly seen convergence rate analysis based on one-dimensional system and mean-value theorem, this analysis is based on multidimensional system and polynomial approximation. A perturbed orbit propagation problem, a radiation heat transfer problem and an orbit optimization problem are solved to verify the validity and efficiency of the proposed methods. The proposed methods not only outperform the highly optimized MATLAB built-in integrator ODE45 in orbit propagation, but also outperform the bvp4c and bvp5c in solving the boundary value problem. Most importantly, the proposed methods are capable of solving some problems where the bvp4c and bvp5c are invalid. The quasi-linearization-collocation method proved to be a highly accurate, efficient and stable solver for ordinary differential equations with various types of constraints.

High‐order isogeometric modified method of characteristics for two‐dimensional coupled Burgers' equations

AbstractThis paper presents a novel isogeometric modified method of characteristics for the numerical solution of the two‐dimensional nonlinear coupled Burgers' equations. The method combines the modified method of characteristics and the high‐order NURBS () elements to discretize the governing equations. The Lagrangian interpretation in this isogeometric analysis greatly reduces the time truncation errors in the Eulerian methods. A third‐order explicit Runge–Kutta scheme is used for the discretization in time. We present a detailed description of the algorithm used for the calculation of departure points and the interpolation stage. Our focus is on constructing highly accurate and stable solvers for the two‐dimensional nonlinear coupled Burgers' equations at high Reynolds numbers. A variety of benchmark tests and numerical examples are provided to show the effectiveness, accuracy, and performance of the proposed modified method of characteristics by virtue of potential advantages of isogeometric analysis. The method developed is anticipated to provide new research directions to the practical calculation of incompressible flows and to studies of their physical behavior.

Entropy Stable Discontinuous Galerkin Finite Element Method with Multi-Dimensional Slope Limitation for Euler Equations

Abstract We present an entropy stable Discontinuous Galerkin (DG) finite element method to approximate systems of 2-dimensional symmetrizable conservation laws on unstructured grids. The scheme is constructed using a combination of entropy conservative fluxes and entropy-stable numerical dissipation operators. The method is designed to work on structured as well as on unstructured meshes. As solutions of hyperbolic conservation laws can develop discontinuities (shocks) in finite time, we include a multidimensional slope limitation step to suppress spurious oscillations in the vicinity of shocks. The numerical scheme has two steps: the first step is a finite element calculation which includes calculations of fluxes across the edges of the elements using 1-D entropy stable solver. The second step is a procedure of stabilization through a truly multi-dimensional slope limiter. We compared the Entropy Stable Scheme (ESS) versus Roe’s solvers associated with entropy corrections and Osher’s solver. The method is illustrated by computing solution of the two stationary problems: a regular shock reflection problem and a 2-D flow around a double ellipse at high Mach number.

Stable solvers for real-time Complex Langevin

This study explores the potential of modern implicit solvers for stochastic partial differential equations in the simulation of real-time complex Langevin dynamics. Not only do these methods offer asymptotic stability, rendering the issue of runaway solution moot, but they also allow us to simulate at comparatively large Langevin time steps, leading to lower computational cost. We compare different ways of regularizing the underlying path integral and estimate the errors introduced due to the finite Langevin time steps. Based on that insight, we implement benchmark (non-)thermal simulations of the quantum anharmonic oscillator on the canonical Schwinger-Keldysh contour of short real-time extent.

A parallel-in-time approach for wave-type PDEs

Numerical solutions to wave-type PDEs utilizing method-of-lines require the ODE solver's stability domain to include a large stretch of the imaginary axis surrounding the origin. We show here that extrapolation based solvers of Gragg-Bulirsch-Stoer (GBS) type can meet this requirement. Extrapolation methods utilize several independent time stepping sequences, making them highly suited for parallel execution. Traditional extrapolation schemes use all time stepping sequences to maximize the method's order of accuracy. The present method instead maintains a desired order of accuracy while employing additional time stepping sequences to shape the resulting stability domain. We optimize the extrapolation coefficients to maximize the stability domain's imaginary axis coverage. This yields a family of explicit schemes that approaches maximal time step size for wave propagation problems. On a computer with several cores we achieve both high order and fast time to solution compared with traditional ODE integrators.

Efficient modeling of three-dimensional convection–diffusion problems in stationary flows

A large number of transport processes in physics, chemistry, and engineering are described by a convection–diffusion equation. This equation is notoriously difficult to solve due to the presence of convection-related first-order gradient differential operators. We describe a new and efficient numerical method for solving the convection–diffusion equation for laminar flows within channels of arbitrary cross section. It is based on reducing the convection–diffusion equation to a set of pure diffusion equations with a complex-valued potential for which fast and numerically stable solvers are readily available. Additionally, we use an eigenvector projection method that allows us to compute snapshots of full concentration distributions over millions of finite elements within a few seconds using a conventional state-of-the-art desktop computer. Our results will be important for all applications where diffusion and convection are both important for correctly describing material transport.

Influence of the convective energy formulation for melting problems with enthalpy methods

The simulation of melting phenomena in the frame work of continuum theory can be handled through many different approaches, among which the fixed-grid methods are the most common ones. Even though such methods have been in use for several decades, only scarce publications exist addressing their performance when solving moving boundary problems. Recently, five of the most used energy formulations with a strong coupling between temperature and enthalpy have been studied in an extensive benchmark both qualitatively and quantitatively. The present paper extends this work by including the role of the formulation of the convective term in the energy equation. Two cases are considered: a first one where all the enthalpy is convected, while in the second one it is the sensible enthalpy only. Similarly to what was done in the aforementioned study, this analysis also includes investigations concerning the modeling of the phase transition, through the associated temperature range for the mushy zone, and several numerical parameters, namely the time step, the mesh coarsening, the CFL condition and the tolerance for the energy equation. The numerical results are compared to a well known quasi-2D melting experiment from the literature in order to determine clearly wrong results. In addition, the liquid fraction and its variance are analyzed for every solver. Except for the solver using an apparent heat capacity method, all of them give reasonable results for a broad range of parameters. The sensible enthalpy convection formulation is in general more stable than the all enthalpy convection formulation. For the chosen parameters and the most stable solver - which is practically not affected by the tolerance, the CFL or the maximum time step - it appears that parameters having the strongest influence on the liquid fraction are the width of the mushy zone followed by the mesh and the convective formulation.

Machine-learning assisted steady-state profile predictions using global optimization techniques

Predicting plasma profiles with a stiff turbulent transport model is important for experimental analysis and development of operation scenarios. Due to the sensitivity of turbulent fluxes to profile gradients, robust predictions are still arduous with a stiff model incorporated in a conventional transport code. With global optimization techniques employed, the new steady-state transport code, global optimization version of the transport equation stable solver, has been developed to overcome these difficulties. It enables us to attain smooth profiles of diffusivity and temperature even though jagged profiles thereof are inclined to emerge in simulations with a stiff model. A neural-network-based surrogate model of a transport model is developed to compensate slow computation inherent to global optimization. Hyperparameter optimization realizes the surrogate model with very good accuracy.

On the numerical stability in dynamical distributed simulations

This work takes aim at studying numerical stability in distributed simulations through dynamic stability and stability criteria for explicit solvers. This is done by studying outer stability limits, for example stability conditions when handling unstable subsystems or marginally stable solvers. To conclude global stability of a distributed system simulation both dynamic stability and solver stability must hold, and this work combines these stability criteria into one unified criterion for distributed linear dynamical systems. Some examples are given in order to both highlight numerical stability issues and to prove stability in different case studies. The derived stability criterion is also extended to include distributed systems containing nonlinear dynamics.

A Numerically Stable Solver for Positive Semidefinite Quadratic Programs Based on Nonnegative Least Squares

This paper proposes a new algorithm for solving convex quadratic programs (QP) subject to linear inequality and equality constraints. The method extends an approach recently proposed by the author for solving strictly convex QP's using nonnegative least squares, by making it numerically more robust and able to handle also the nonstrictly convex case, equality constraints, and warm starting from an initial guess. Robustness is achieved by introducing an outer proximal-point iteration scheme that regularizes the problem without altering the solution, and by adaptively weighting the least squares problems encountered while solving the problem. The performance of the resulting QP solver in terms of speed and robustness in the single precision arithmetic is assessed against other optimization algorithms tailored to embedded model predictive control applications.

A frequency domain finite element solver for acoustic simulations of 3D rooms with microperforated panel absorbers

Microperforated panel (MPP) absorbers, which provide broadband sound absorption without the use of fibrous materials, have favorable material properties that support recyclability, flexibility of design, hygiene demands, and cleaning. Many earlier studies have specifically examined the development of absorbers themselves. However, to use the absorption performance of MPP absorbers sufficiently in room acoustic applications, it is beneficial to develop accurate prediction methods of sound fields in rooms with MPP absorbers. Such methods are expected to be useful for room acoustics design and absorber design tools. This study constructs a frequency domain finite element (FE) solver for acoustic simulations of a practical sized room with MPP absorbers. Then the accuracy and effectiveness are evaluated. In the FE solver, spatial domains are discretized by fourth-order accurate FEs in terms of dispersion error, and MPP absorbers are modeled using first-order hexahedral limp MPP elements that can deal with sound propagation in the backing structure of absorbers. First, the accuracy of present FE solver is demonstrated using impedance tube problems in comparison with conventional second-order accurate FEs. Results show higher convergence of solutions for the present FE solver. Then, exploration of an iterative solver for efficient multi-frequency analyses reveals that the recently developed CSQMOR is a faster and more stable solver. Finally, comparison with a conventional surface impedance model based on a locally reacting assumption confirms the effectiveness of present FE solver by presenting the importance of dealing with the incident angle dependence of reactance of a rigid-backed air cavity in the modeling of single-leaf MPP absorbers.

Application of a wall function to simulate turbulent flows in foil bearings at high rotational speeds

For turbulent simulations in hydrodynamic applications, a fine mesh close to wall boundaries is required to correctly predict the friction and heat fluxes. This results in a remarkable increase of computational cost due to the high aspect ratio cells and reduced allowable time step. To reduce the computational cost, a compressible wall function is introduced and validated with the representative flows. Moreover, high aspect ratio cells lead to simulation instabilities. These spurious oscillations are smoothed through the addition of a fourth-order artificial dissipation term. To verify that the solution accuracy is not affected, the method of manufactured solution is applied. These two additions result in a fast and stable solver for turbulent simulations of foil thrust bearings.

Why Do Field-Based Methods Fail to Model Plasmonics?

The paper studies plasmonics modeling issues and examines the reasons behind the failure of the field-based methods relying on Pade approximations widely used in the analysis of photonic devices based on dielectric materials. Through a study of evanescent, radiation, guided, and surface modes of a plasmonic structure where the failure appears clearly, we demonstrate the physical explanation of this failure and suggest some remedies. We developed a Bidirectional Beam Propagation Method (BiBPM) by adopting a Blocked Schur (BS) algorithm to introduce an unconditionally stable method for plasmonic structures with strong discontinuities. Central to BiBPMs is the accurate calculation of the square root operators that is very widely performed using Pade approximations. However, recent reports demonstrate convergence of Pade that is too slow to lend itself a stable solver in plasmonics [1] . Moreover, Pade approximations completely fail in handling such a strong discontinuity between dielectric and plasmonic waveguides, where a very-wide spectrum of modes could be excited. Alternatively, we propose calculating these operators by the twice faster BS algorithm. Beyond the computational speed, our suggested approach overbears the Pade-based BiBPMs instability and accuracy problems, thanks to the proper physical treatment of surface and evanescent waves: the notorious sources of instability. Through the plasmonic discontinuity problems, the superiority of BS approach has been determined numerically and explained physically.

Method of Lines Transpose: An Efficient Unconditionally Stable Solver for Wave Propagation

Building upon recent results obtained in Causley and Christlieb (SIAM J Numer Anal 52(1):220–235, 2014), Causley et al. (Math Comput 83(290):2763–2786, 2014, Method of lines transpose: high order L-stable O(N) schemes for parabolic equations using successive convolution, 2015), we describe an efficient second-order, unconditionally stable scheme for solving the wave equation, based on the method of lines transpose (MOL $$^T$$ ), and the resulting semi-discrete (i.e. continuous in space) boundary value problem. In Causley and Christlieb (SIAM J Numer Anal 52(1):220–235, 2014), unconditionally stable schemes of high order were derived, and in Causley et al. (Method of lines transpose: high order L-stable O(N) schemes for parabolic equations using successive convolution, 2015) a high order, fast $$\mathcal {O}(N)$$ spatial solver was derived, which is matrix-free and is based on dimensional-splitting. In this work, are interested in building a wave solver, and our main concern is the development of boundary conditions. We demonstrate all desired boundary conditions for a wave solver, including outflow boundary conditions, in 1D and 2D. The scheme works in a logically Cartesian fashion, and the boundary points are embedded into the regular mesh, without incurring stability restrictions, so that boundary conditions are imposed without any reduction in the order of accuracy. We demonstrate how the embedded boundary approach works in the cases of Dirichlet and Neumann boundary conditions. Further, we develop outflow and periodic boundary conditions for the MOL $$^T$$ formulation. Our solver is designed to couple with particle codes, and so special attention is also paid to the implementation of point sources, and soft sources which can be used to launch waves into waveguides.

A particle-in-cell method for the simulation of plasmas based on an unconditionally stable field solver

We propose a new particle-in-cell (PIC) method for the simulation of plasmas based on a recently developed, unconditionally stable solver for the wave equation. This method is not subject to a CFL restriction, limiting the ratio of the time step size to the spatial step size, typical of explicit methods, while maintaining computational cost and code complexity comparable to such explicit schemes. We describe the implementation in one and two dimensions for both electrostatic and electromagnetic cases, and present the results of several standard test problems, showing good agreement with theory with time step sizes much larger than allowed by typical CFL restrictions.

Optimization-based additive decomposition of weakly coercive problems with applications

We present an abstract mathematical framework for an optimization-based additive decomposition of a large class of variational problems into a collection of concurrent subproblems. The framework replaces a given monolithic problem by an equivalent constrained optimization formulation in which the subproblems define the optimization constraints and the objective is to minimize the mismatch between their solutions. The significance of this reformulation stems from the fact that one can solve the resulting optimality system by an iterative process involving only solutions of the subproblems. Consequently, assuming that stable numerical methods and efficient solvers are available for every subproblem, our reformulation leads to robust and efficient numerical algorithms for a given monolithic problem by breaking it into subproblems that can be handled more easily. An application of the framework to the Oseen equations illustrates its potential.