Associated to a finite measure on the real line with finite moments are recurrence coefficients in a three-term formula for orthogonal polynomials with respect to this measure. These recurrence coefficients are frequently inputs to modern computational tools that facilitate evaluation and manipulation of polynomials with respect to the measure, and such tasks are foundational in numerical approximation and quadrature. Although the recurrence coefficients for classical measures are known explicitly, those for nonclassical measures must typically be numerically computed. We survey and review existing approaches for computing these recurrence coefficients for univariate orthogonal polynomial families and propose a novel "predictor-corrector" algorithm for a general class of continuous measures. We combine the predictor-corrector scheme with a stabilized Lanczos procedure for a new hybrid algorithm that computes recurrence coefficients for a fairly wide class of measures that can have both continuous and discrete parts. We evaluate the new algorithms against existing methods in terms of accuracy and efficiency.

We consider a least-squares/relaxation finite element method for the numerical solution of the prescribed Jacobian equation. We look for its solution via a least-squares approach. We introduce a relaxation algorithm that decouples this least-squares problem into a sequence of local nonlinear problems and variational linear problems. We develop dedicated solvers for the algebraic problems based on Newton's method and we solve the differential problems using mixed low-order finite elements. Various numerical experiments demonstrate the accuracy, efficiency and the robustness of the proposed method, compared for instance to augmented Lagrangian approaches.

We analyze a class of meshfree semi-Lagrangian methods for solving advection problems on smooth, closed surfaces with solenoidal velocity field. In particular, we prove the existence of an embedding equation whose corresponding semi-Lagrangian methods yield the ones in the literature for solving problems on surfaces. Our analysis allows us to apply standard bulk domain convergence theories to the surface counterparts. In addition, we provide detailed descriptions for implementing the proposed methods to run on point clouds. After verifying the convergence rates against the theory, we show that the proposed method is a robust building block for more complicated problems, such as advection problems with non-solenoidal velocity field, inviscid Burgers' equations and systems of reaction advection diffusion equations for pattern formation.

In this manuscript a POD-Galerkin based Reduced Order Model for unsteady Fluid-Structure Interaction problems is presented. The model is based on a partitioned algorithm, with semi-implicit treatment of the coupling conditions. A Chorin-Temam projection scheme is applied to the incompressible Navier-Stokes problem, and a Robin coupling condition is used for the coupling between the fluid and the solid. The coupled problem is based on an Arbitrary Lagrangian Eulerian formulation, and the Proper Orthogonal Decomposition procedure is used for the generation of the reduced basis. We extend existing works on a segregated Reduced Order Model for Fluid-Structure Interaction to unsteady problems that couple an incompressible, Newtonian fluid with a linear elastic solid, in two spatial dimensions. We consider three test cases to assess the overall capabilities of the method: an unsteady, non-parametrized problem, a problem that presents a geometrical parametrization of the solid domain, and finally, a problem where a parametrization of the solid's shear modulus is taken into account.

The regularized optimal mass transport (rOMT) problem adds a diffusion term to the continuity equation in the original dynamic formulation of the optimal mass transport (OMT) problem proposed by Benamou and Brenier. We show that the rOMT model serves as a powerful tool in computational fluid dynamics for visualizing fluid flows in the glymphatic system. In the present work, we describe how to modify the previous numerical method for efficient implementation, resulting in a significant reduction in computational runtime. Numerical results applied to synthetic and real-data are provided.

This paper is devoted to the numerical solution of constrained energy minimization problems arising in computational physics and chemistry such as the Gross-Pitaevskii and Kohn-Sham models. In particular, we introduce Riemannian Newton methods on the infinite-dimensional Stiefel and Grassmann manifolds. We study the geometry of these two manifolds, its impact on the Newton algorithms, and present expressions of the Riemannian Hessians in the infinite-dimensional setting, which are suitable for variational spatial discretizations. A series of numerical experiments illustrates the performance of the methods and demonstrates their supremacy compared to other well-established schemes such as the self-consistent field iteration and gradient descent schemes.

A new relaxed variant of interior point method for low-rank semidefinite programming problems is proposed in this paper. The method is a step outside of the usual interior point framework. In anticipation to converging to a low-rank primal solution, a special nearly low-rank form of all primal iterates is imposed. To accommodate such a (restrictive) structure, the first order optimality conditions have to be relaxed and are therefore approximated by solving an auxiliary least-squares problem. The relaxed interior point framework opens numerous possibilities how primal and dual approximated Newton directions can be computed. In particular, it admits the application of both the first- and the second-order methods in this context. The convergence of the method is established. A prototype implementation is discussed and encouraging preliminary computational results are reported for solving the SDP-reformulation of matrix-completion problems.

The Loewner framework is one of the most successful data-driven model order reduction techniques. If is the cardinality of a given data set, the so-called Loewner and shifted Loewner matrices and can be defined by solely relying on information encoded in the considered data set and they play a crucial role in the computation of the sought rational model approximation.In particular, the singular value decomposition of a linear combination of and provides the tools needed to construct accurate models which fulfill important approximation properties with respect to the original data set. However, for highly-sampled data sets, the dense nature of and leads to numerical difficulties, namely the failure to allocate these matrices in certain memory-limited environments or excessive computational costs. Even though they do not possess any sparsity pattern, the Loewner and shifted Loewner matrices are extremely structured and, in this paper, we show how to fully exploit their Cauchy-like structure to reduce the cost of computing accurate rational models while avoiding the explicit allocation of and . In particular, the use of the format allows us to remarkably lower both the computational cost and the memory requirements of the Loewner framework obtaining a novel scheme whose costs scale with .

We develop new multilevel Monte Carlo (MLMC) methods to estimate the expectation of the smallest eigenvalue of a stochastic convection-diffusion operator with random coefficients. The MLMC method is based on a sequence of finite element (FE) discretizations of the eigenvalue problem on a hierarchy of increasingly finer meshes. For the discretized, algebraic eigenproblems we use both the Rayleigh quotient (RQ) iteration and implicitly restarted Arnoldi (IRA), providing an analysis of the cost in each case. By studying the variance on each level and adapting classical FE error bounds to the stochastic setting, we are able to bound the total error of our MLMC estimator and provide a complexity analysis. As expected, the complexity bound for our MLMC estimator is superior to plain Monte Carlo. To improve the efficiency of the MLMC further, we exploit the hierarchy of meshes and use coarser approximations as starting values for the eigensolvers on finer ones. To improve the stability of the MLMC method for convection-dominated problems, we employ two additional strategies. First, we consider the streamline upwind Petrov-Galerkin formulation of the discrete eigenvalue problem, which allows us to start the MLMC method on coarser meshes than is possible with standard FEs. Second, we apply a homotopy method to add stability to the eigensolver for each sample. Finally, we present a multilevel quasi-Monte Carlo method that replaces Monte Carlo with a quasi-Monte Carlo (QMC) rule on each level. Due to the faster convergence of QMC, this improves the overall complexity. We provide detailed numerical results comparing our different strategies to demonstrate the practical feasibility of the MLMC method in different use cases. The results support our complexity analysis and further demonstrate the superiority over plain Monte Carlo in all cases.

We present an algorithm for fast generation of quasi-uniform and variable-spacing nodes on domains whose boundaries are represented as computer-aided design (CAD) models, more specifically non-uniform rational B-splines (NURBS). This new algorithm enables the solution of partial differential equations within the volumes enclosed by these CAD models using (collocation-based) meshless numerical discretizations. Our hierarchical algorithm first generates quasi-uniform node sets directly on the NURBS surfaces representing the domain boundary, then uses the NURBS representation in conjunction with the surface nodes to generate nodes within the volume enclosed by the NURBS surface. We provide evidence for the quality of these node sets by analyzing them in terms of local regularity and separation distances. Finally, we demonstrate that these node sets are well-suited (both in terms of accuracy and numerical stability) for meshless radial basis function generated finite differences discretizations of the Poisson, Navier-Cauchy, and heat equations. Our algorithm constitutes an important step in bridging the field of node generation for meshless discretizations with isogeometric analysis.

We introduce an -version discontinuous Galerkin finite element method (DGFEM) for the linear Boltzmann transport problem. A key feature of this new method is that, while offering arbitrary order convergence rates, it may be implemented in an almost identical form to standard multigroup discrete ordinates methods, meaning that solutions can be computed efficiently with high accuracy and in parallel within existing software. This method provides a unified discretisation of the space, angle, and energy domains of the underlying integro-differential equation and naturally incorporates both local mesh and local polynomial degree variation within each of these computational domains. Moreover, general polytopic elements can be handled by the method, enabling efficient discretisations of problems posed on complicated spatial geometries. We study the stability and -version a priori error analysis of the proposed method, by deriving suitable -approximation estimates together with a novel inf-sup bound. Numerical experiments highlighting the performance of the method for both polyenergetic and monoenergetic problems are presented.

This paper studies two hybrid discontinuous Galerkin (HDG) discretizations for the velocity-density formulation of the compressible Stokes equations with respect to several desired structural properties, namely provable convergence, the preservation of non-negativity and mass constraints for the density, and gradient-robustness. The later property dramatically enhances the accuracy in well-balanced situations, such as the hydrostatic balance where the pressure gradient balances the gravity force. One of the studied schemes employs an -conforming velocity ansatz space which ensures all mentioned properties, while a fully discontinuous method is shown to satisfy all properties but the gradient-robustness. Also higher-order schemes for both variants are presented and compared in three numerical benchmark problems. The final example shows the importance also for non-hydrostatic well-balanced states for the compressible Navier-Stokes equations.

We propose two implicit numerical schemes for the low-rank time integration of stiff nonlinear partial differential equations. Our approach uses the preconditioned Riemannian trust-region method of Absil, Baker, and Gallivan, 2007. We demonstrate the efficiency of our method for solving the Allen-Cahn and the Fisher-KPP equations on the manifold of fixed-rank matrices. Our approach allows us to avoid the restriction on the time step typical of methods that use the fixed-point iteration to solve the inner nonlinear equations. Finally, we demonstrate the efficiency of the preconditioner on the same variational problems presented in Sutti and Vandereycken, 2021.

We present a higher-order space-time adaptive method for the numerical solution of the Richards equation that describes a flow motion through variably saturated media. The discretization is based on the space-time discontinuous Galerkin method, which provides high stability and accuracy and can naturally handle varying meshes. We derive reliable and efficient a posteriori error estimates in the residual-based norm. The estimates use well-balanced spatial and temporal flux reconstructions which are constructed locally over space-time elements or space-time patches. The accuracy of the estimates is verified by numerical experiments. Moreover, we develop the -adaptive method and demonstrate its efficiency and usefulness on a practically relevant example.

We extend the Adaptive Antoulas-Anderson (AAA) algorithm to develop a data-driven modeling framework for linear systems with quadratic output (LQO). Such systems are characterized by two transfer functions: one corresponding to the linear part of the output and another one to the quadratic part. We first establish the joint barycentric representations and the interpolation theory for the two transfer functions of LQO systems. This analysis leads to the proposed AAA-LQO algorithm. We show that by interpolating the transfer function values on a subset of samples together with imposing a least-squares minimization on the rest, we construct reliable data-driven LQO models. Two numerical test cases illustrate the efficiency of the proposed method.

This paper presents a comprehensive study of non-linear systems of multi-order fractional differential equations from aspects of theory and numerical approximation. In this regard, we first establish the well-posedness of the underlying problem by investigations concerning the existence, uniqueness, and influence of perturbed data on the behavior of the solutions as well as smoothness of the solutions under some assumptions on the given data. Next, from the numerical perspective, we develop and analyze a well-conditioned and high-order numerical approach based on an operational spectral Galerkin method. The main advantage of our strategy is that it characterizes the approximate solution via some recurrence formulas, instead of solving a complex non-linear block algebraic system that requires high computational costs. Moreover, the familiar spectral accuracy is justified in a weighted -norm, and some practical test problems are provided to approve efficiency of the proposed method.

We introduce two new lowest order methods, a mixed method, and a hybrid discontinuous Galerkin method, for the approximation of incompressible flows. Both methods use divergence-conforming linear Brezzi-Douglas-Marini space for approximating the velocity and the lowest order Raviart-Thomas space for approximating the vorticity. Our methods are based on the physically correct viscous stress tensor of the fluid, involving the symmetric gradient of velocity (rather than the gradient), provide exactly divergence-free discrete velocity solutions, and optimal error estimates that are also pressure robust. We explain how the methods are constructed using the minimal number of coupling degrees of freedom per facet. The stability analysis of both methods are based on a Korn-like inequality for vector finite elements with continuous normal component. Numerical examples illustrate the theoretical findings and offer comparisons of condition numbers between the two new methods.

We use the behavior of the norm of the solutions of linear hyperbolic equations with discontinuous coefficient matrices as a surrogate to infer stability of discontinuous Galerkin spectral element methods (DGSEM). Although the norm is not bounded in terms of the initial data for homogeneous and dissipative boundary conditions for such systems, the norm is easier to work with than a norm that discounts growth due to the discontinuities. We show that the DGSEM with an upwind numerical flux that satisfies the Rankine-Hugoniot (or conservation) condition has the same energy bound as the partial differential equation does in the norm, plus an added dissipation that depends on how much the approximate solution fails to satisfy the Rankine-Hugoniot jump.

In this work we study various continuous finite element discretization for two dimensional hyperbolic partial differential equations, varying the polynomial space (Lagrangian on equispaced, Lagrangian on quadrature points () and Bernstein), the stabilization techniques (streamline-upwind Petrov-Galerkin, continuous interior penalty, orthogonal subscale stabilization) and the time discretization (Runge-Kutta (RK), strong stability preserving RK and deferred correction). This is an extension of the one dimensional study by Michel et al. (J Sci Comput 89(2):31, 2021. 10.1007/s10915-021-01632-7), whose results do not hold in multi-dimensional frameworks. The study ranks these schemes based on efficiency (most of them are mass-matrix free), stability and dispersion error, providing the best CFL and stabilization coefficients. The challenges in two-dimensions are related to the Fourier analysis. Here, we perform it on two types of periodic triangular meshes varying the angle of the advection, and we combine all the results for a general stability analysis. Furthermore, we introduce additional high order viscosity to stabilize the discontinuities, in order to show how to use these methods for tests of practical interest. All the theoretical results are thoroughly validated numerically both on linear and non-linear problems, and error-CPU time curves are provided. Our final conclusions suggest that elements combined with SSPRK and OSS stabilization is the most promising combination.

We propose a novel use of generative adversarial networks (GANs) (i) to make predictions in time (PredGAN) and (ii) to assimilate measurements (DA-PredGAN). In the latter case, we take advantage of the natural adjoint-like properties of generative models and the ability to simulate forwards and backwards in time. GANs have received much attention recently, after achieving excellent results for their generation of realistic-looking images. We wish to explore how this property translates to new applications in computational modelling and to exploit the adjoint-like properties for efficient data assimilation. We apply these methods to a compartmental model in epidemiology that is able to model space and time variations, and that mimics the spread of COVID-19 in an idealised town. To do this, the GAN is set within a reduced-order model, which uses a low-dimensional space for the spatial distribution of the simulation states. Then the GAN learns the evolution of the low-dimensional states over time. The results show that the proposed methods can accurately predict the evolution of the high-fidelity numerical simulation, and can efficiently assimilate observed data and determine the corresponding model parameters.

We propose a novel Model Order Reduction framework that is able to handle solutions of hyperbolic problems characterized by multiple travelling discontinuities. By means of an optimization based approach, we introduce suitable calibration maps that allow us to transform the original solution manifold into a lower dimensional one. The novelty of the methodology is represented by the fact that the optimization process does not require the knowledge of the discontinuities location. The optimization can be carried out simply by choosing some reference control points, thus avoiding the use of some implicit shock tracking techniques, which would translate into an increased computational effort during the offline phase. In the online phase, we rely on a non-intrusive approach, where the coefficients of the projection of the reduced order solution onto the reduced space are recovered by means of an Artificial Neural Network. To validate the methodology, we present numerical results for the 1D Sod shock tube problem, for the 2D double Mach reflection problem, also in the parametric case, and for the triple point problem.