In this work, a family of symmetric interpolation points are generated on the four-dimensional simplex (i.e., the pentatope). Interpolation points are crucial toward the furtherance of space–time finite element methods, which themselves find applications in shock fitting, fluid structure interactions, and even rotating detonation engines. The points generated herein are optimized in order to minimize the Lebesgue constant. The process of generating these points closely follows that outlined by Warburton (J Eng Math 56:247–262, 2006). Here, Warburton generated optimal interpolation points on the triangle and tetrahedron by formulating explicit geometric warping and blending functions and applying these functions to equidistant nodal distributions. The locations of the resulting points were Lebesgue-optimized. In our work, we extend this procedure to four dimensions and construct interpolation points on the pentatope up to order 10. The Lebesgue constants of our nodal sets are calculated and are shown to outperform those of equidistant nodal distributions, as well as many of the point distributions of Isaac (SIAM J Sci Comput 42(6):4046–4062, 2020).

Among numerous first-order algorithms, the Fast Iterative Shrinkage-Thresholding Algorithm, known as FISTA, is renowned for its convergence speed of O ( 1 k 2 ) in terms of the objective function value, where k denotes the number of iterations. Additionally, various improvements to FISTA have been proposed in the literature. Among them, the convergence rate of a parameterized FISTA (PFISTA) proposed by Liang, Luo and Schonlieb [SIAM J Sci Comput. 2022:44(3):A1069–A1091] is obviously faster than that of the original FISTA, which can reach o ( 1 k 2 ) convergence rate. Building upon the idea of parameterization, we introduce the fully parameterized APFISTA for complete parameterization of the inertia term, aiming to enhance the current situation of partially parameterized inertia terms. We establish the convergence speed of its objective function value and the sequence generated by the algorithm. Finally, we apply the proposed method to solve the least squares problem and linear inverse problems. The obtained numerical results illustrate the practical behavior and theoretical analysis of the proposed approach.

In this study, we present an hp-multigrid preconditioner for a divergence-conforming HDG scheme for the generalized Stokes and the Navier–Stokes equations using an augmented Lagrangian formulation. Our method relies on conforming simplicial meshes in two- and three-dimensions. The hp-multigrid algorithm is a multiplicative auxiliary space preconditioner that employs the lowest-order space as the auxiliary space, and we develop a geometric multigrid method as the auxiliary space solver. For the generalized Stokes problem, the crucial ingredient of the geometric multigrid method is the equivalence between the condensed lowest-order divergence-conforming HDG scheme and a Crouzeix–Raviart discretization with a pressure-robust treatment as introduced in Linke and Merdon (Comput Methods Appl Mech Engrg 311:304–326, 2022), which allows for the direct application of geometric multigrid theory on the Crouzeix–Raviart discretization. The numerical experiments demonstrate the robustness of the proposed hp-multigrid preconditioner with respect to mesh size and augmented Lagrangian parameter, with iteration counts insensitivity to polynomial order increase. Inspired by the works by Benzi and Olshanskii (SIAM J Sci Comput 28:2095–2113, 2006) and Farrell et al. (SIAM J Sci Comput 41:A3073–A3096, 2019), we further test the proposed preconditioner on the divergence-conforming HDG scheme for the Navier–Stokes equations. Numerical experiments show a mild increase in the iteration counts of the preconditioned GMRes solver with the rise in Reynolds number up to 103\\documentclass[12pt]{minimal} \\usepackage{amsmath} \\usepackage{wasysym} \\usepackage{amsfonts} \\usepackage{amssymb} \\usepackage{amsbsy} \\usepackage{mathrsfs} \\usepackage{upgreek} \\setlength{\\oddsidemargin}{-69pt} \\begin{document}$$10^3$$\\end{document}.

The solution of time-dependent hyperbolic conservation laws on cut cell meshes causes the small cell problem: standard schemes are not stable on the arbitrarily small cut cells if an explicit time stepping scheme is used and the time step size is chosen based on the size of the background cells. In May and Berger (J Sci Comput 71: 919–943, 2017), the mixed explicit-implicit approach in general and MUSCL-Trap (monotonic upwind scheme for conservation laws and trapezoidal scheme) in particular have been introduced to solve this problem by using implicit time stepping on the cut cells. Theoretical and numerical results have indicated that this might lead to a loss in accuracy when switching between the explicit and implicit time stepping. In this contribution, we examine this in more detail and will prove in one dimension that the specific combination MUSCL-Trap of an explicit second-order and an implicit second-order scheme results in a fully second-order mixed scheme. As this result is unlikely to hold in two dimensions, we also introduce two new versions of mixed explicit-implicit schemes based on exchanging the explicit scheme. We present numerical tests in two dimensions where we compare the new versions with the original MUSCL-Trap scheme.

Balanced truncation is a well-established model order reduction method which has been applied to a variety of problems. Recently, a connection between linear Gaussian Bayesian inference problems and the system-theoretic concept of balanced truncation has been drawn (Qian et al in Sci Comput 91:29, 2022). Although this connection is new, the application of balanced truncation to data assimilation is not a novel idea: it has already been used in four-dimensional variational data assimilation (4D-Var). This paper discusses the application of balanced truncation to linear Gaussian Bayesian inference, and, in particular, the 4D-Var method, thereby strengthening the link between systems theory and data assimilation further. Similarities between both types of data assimilation problems enable a generalisation of the state-of-the-art approach to the use of arbitrary prior covariances as reachability Gramians. Furthermore, we propose an enhanced approach using time-limited balanced truncation that allows to balance Bayesian inference for unstable systems and in addition improves the numerical results for short observation periods.

We present an extended range of stable flux reconstruction (FR) methods on triangles through the development and application of the summation-by-parts framework in two-dimensions. This extended range of stable schemes is then shown to contain the single parameter schemes of Castonguay et al. (J Sci Comput 51:224–256, 2011) on triangles, and our definition enables wider stability bounds to be developed for those single parameter families. Stable upwinded spectral difference (SD) schemes on triangular elements have previously been found using Fourier analysis. We used our extended range of FR schemes to investigate the linear stability of SD methods on triangles, and it was found that a only first order SD scheme could be recovered within this set of FR methods.

The aim of this work is to apply a semi-implicit (SI) strategy in an implicit-explicit (IMEX) Runge–Kutta (RK) setting introduced in Boscarino et al. (J Sci Comput 68:975–1001, 2016) to a sequence of 1D time-dependent partial differential equations (PDEs) with high order spatial derivatives. This strategy gives a great flexibility to treat these equations, and allows the construction of simple linearly implicit schemes without any Newton’s iteration. Furthermore, the SI IMEX-RK schemes so designed does not need any severe time step restriction that usually one has using explicit methods for the stability, i.e. Delta t = {mathcal {O}}(Delta t^k) for the kth (k ge 2) order PDEs. For the space discretization, this strategy is combined with finite difference schemes. We illustrate the effectiveness of the schemes with many applications to dissipative, dispersive and biharmonic-type equations. Numerical experiments show that the proposed schemes are stable and can achieve optimal orders of accuracy.

We develop flux globalization based well-balanced central-upwind schemes for hydrodynamic equations with general free energy. The proposed schemes are well-balanced in the sense that they are capable of exactly preserving quite complicated steady-state solutions and also exactly capturing traveling waves, even when vacuum regions are present. In order to accurately track interfaces of the vacuum regions and near vacuum parts of the solution, we use the technique introduced in Chertock et al. (J Sci Comput 90:Paper No. 9, 2022) and design a hybrid approach: inside the no vacuum regions, we use the flux globalization based well-balanced central-upwind scheme, while elsewhere we implement the central-upwind scheme similar to the one proposed in Bollermann et al. (J Sci Comput 56:267–290, 2013) in the context of wet/dry fronts in the shallow water equations. The advantages of the proposed schemes are demonstrated on a number of challenging numerical examples.

Within the last years pressure robust methods for the discretization of incompressible fluids have been developed. These methods allow the use of standard finite elements for the solution of the problem while simultaneously removing a spurious pressure influence in the approximation error of the velocity of the fluid, or the displacement of an incompressible solid. To this end, reconstruction operators are utilized mapping discretely divergence free functions to divergence free functions. This work shows that the modifications proposed for Stokes equation by Linke (Comput Methods Appl Mech Eng 268:782–800, 2014) also yield gradient robust methods for nearly incompressible elastic materials without the need to resort to discontinuous finite elements methods as proposed in Fu et al. (J Sci Comput 86(3):39–30, 2021).

A new family of non-hydrostatic layer-averaged models for the non-stationary Euler equations is presented in this work, with improved dispersion relations. They are a generalisation of the layer-averaged models introduced in Fernández-Nieto et al. (Commun Math Sci 16(05):1169–1202, 2018), named LDNH models, where the vertical profile of the horizontal velocity is layerwise constant. This assumption implies that solutions of LDNH can be seen as a first order Galerkin approximation of Euler system. Nevertheless, it is not a fully (x, z) Galerkin discretisation of Euler system, but just in the vertical direction (z). Thus, the resulting model only depends on the horizontal space variable (x), and therefore specific and efficient numerical methods can be applied (see Escalante-Sanchez et al. in J Sci Comput 89(55):1–35, 2021). This work focuses on particular weak solutions where the horizontal velocity is layerwise linear on z and possibly discontinuous across layer interfaces. This approach allows the system to be a second-order approximation in the vertical direction of Euler system. Several closure relations of the layer-averaged system with non-hydrostatic pressure are presented. The resulting models are named LIN-NH_k models, with k=0,1,2. Parameter k indicates the degree of the vertical velocity profile considered in the approximation of the vertical momentum equation. All the introduced models satisfy a dissipative energy balance. Finally, an analysis and a comparison of the dispersive properties of each model are carried out. We show that Models LIN-NH_1 and LIN-NH_2 provide a better dispersion relation, group velocity and shoaling than LDNH models.

An hp-mixed discontinuous Galerkin method for the biharmonic eigenvalue problem

We propose a new algorithmic framework for constrained compressed sensing models that admit nonconvex sparsity-inducing regularizers including the log-penalty function as objectives, and nonconvex loss functions such as the Cauchy loss function and the Tukey biweight loss function in the constraint. Our framework employs iteratively reweighted $$\ell _1$$ and $$\ell _2$$ schemes to construct subproblems that can be efficiently solved by well-developed solvers for basis pursuit denoising such as SPGL1 by van den Berg and Friedlander (SIAM J Sci Comput 31:890-912, 2008). We propose a new termination criterion for the subproblem solvers that allows them to return an infeasible solution, with a suitably constructed feasible point satisfying a descent condition. The feasible point construction step is the key for establishing the well-definedness of our proposed algorithm, and we also prove that any accumulation point of this sequence of feasible points is a stationary point of the constrained compressed sensing model, under suitable assumptions. Finally, we compare numerically our algorithm (with subproblems solved by SPGL1 or the alternating direction method of multipliers) against the SCP $$_\textrm{ls}$$ in Yu et al. (SIAM J Optim 31: 2024-2054, 2021) on solving constrained compressed sensing models with the log-penalty function as the objective and the Cauchy loss function in the constraint, for badly scaled measurement matrices. Our computational results show that our approaches return solutions with better recovery errors, and are always faster.

In this work, we investigate the diffusive optical tomography (DOT) problem in the case that limited boundary measurements are available. Motivated by the direct sampling method (DSM) proposed in Chow et al. (SIAM J Sci Comput 37(4):A1658–A1684, 2015), we develop a deep direct sampling method (DDSM) to recover the inhomogeneous inclusions buried in a homogeneous background. In this method, we design a convolutional neural network to approximate the index functional that mimics the underling mathematical structure. The benefits of the proposed DDSM include fast and easy implementation, capability of incorporating multiple measurements to attain high-quality reconstruction, and advanced robustness against the noise. Numerical experiments show that the reconstruction accuracy is improved without degrading the efficiency, demonstrating its potential for solving the real-world DOT problems.

In this paper, numerical experiments are carried out to investigate the impact of penalty parameters in the numerical traces on the resonance errors of high-order multiscale discontinuous Galerkin (DG) methods (Dong et al. in J Sci Comput 66: 321–345, 2016; Dong and Wang in J Comput Appl Math 380: 1–11, 2020) for a one-dimensional stationary Schrödinger equation. Previous work showed that penalty parameters were required to be positive in error analysis, but the methods with zero penalty parameters worked fine in numerical simulations on coarse meshes. In this work, by performing extensive numerical experiments, we discover that zero penalty parameters lead to resonance errors in the multiscale DG methods, and taking positive penalty parameters can effectively reduce resonance errors and make the matrix in the global linear system have better condition numbers.

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 (Cubature) 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. https://doi.org/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 Cubature elements combined with SSPRK and OSS stabilization is the most promising combination.

In this paper, some enhanced error estimates are derived for the augmented subspace methods which are designed for solving eigenvalue problems. For the first time, we strictly prove that the augmented subspace methods have the second order convergence rate between the two iteration steps, which is better than the existing theoretical results in Lin and Xie (Math Comp 84:71–88, 2015), Xie et al. (SIAM J Numer Anal 57(6):2519–2550, 2019), Xu et al. (SIAM J Sci Comput 42(5):A2655–A2677, 2020) and more consistent with the results of actual numerical test. These sharper estimates explicitly depict the dependence of convergence rate on the coarse spaces, which provides new advantages for the augmented subspace methods. Some numerical examples are finally presented to validate these new estimate results and the efficiency of our algorithms.

Multifidelity methods are widely used for estimating quantities of interest (QoI) in computational science by employing numerical simulations of differing costs and accuracies. Many methods approximate numerical-valued statistics that represent only limited information, e.g., scalar statistics, about the QoI. Further quantification of uncertainty, e.g., for risk assessment, failure probabilities, or confidence intervals, requires estimation of the full distributions. In this paper, we generalize the ideas in (Xu et al. in SIAM J Sci Comput 44(1):A150–A175, 2022) to develop a multifidelity method that approximates the full distribution of scalar-valued QoI. The main advantage of our approach compared to alternative methods is that we require no particular relationships among the high and lower-fidelity models (e.g. model hierarchy), and we do not assume any knowledge of model statistics including correlations and other cross-model statistics before the procedure starts. Under suitable assumptions in the framework above, we achieve provable 1-Wasserstein metric convergence of an algorithmically constructed distributional emulator via an exploration–exploitation strategy. We also prove that crucial policy actions taken by our algorithm are budget-asymptotically optimal. Numerical experiments are provided to support our theoretical analysis.

This paper considers a way of generalizing the t-SVD of third-order tensors (regarded as tubal matrices) to tensors of arbitrary order \(N\) [which can be similarly regarded as tubal tensors of order \((N-1)\)]. Such a generalization is different from the t-SVD for tensors of order greater than three (Martin et al. in SIAM J Sci Comput 35(1):A474–A490, 2013). The decomposition is called Hot-SVD, since it can be recognized as a tensor–tensor product version of the celebrated higher order SVD (HOSVD). The existence of Hot-SVD is proved. To this end, the “small-t” transpose for third-order tensors is introduced. This transpose is crucial in the verification of Hot-SVD, since it serves as a bridge between tubal tensors and their unfolding tubal matrices. We establish some properties of Hot-SVD, analogous to those of HOSVD. The truncated Hot-SVD and sequentially truncated Hot-SVD are then introduced, with \(\sqrt{N}\)-error bounds established for an \((N+1)\)-th-order tensor. We provide numerical examples to validate Hot-SVD, truncated Hot-SVD, and sequentially truncated Hot-SVD.

In this work we study the solution of the Riemann problem for the barotropic version of the conservative symmetric hyperbolic and thermodynamically compatible (SHTC) two-phase flow model introduced in Romenski et al. (J Sci Comput 42(1):68, 2009, Quart Appl Math 65(2):259–279, 2007). All characteristic fields are carefully studied and explicit expressions are derived for the Riemann invariants and the Rankine–Hugoniot conditions. Due to the presence of multiple characteristics in the system under consideration, non-standard wave phenomena can occur. Therefore we briefly review admissibility conditions for discontinuities and then discuss possible wave interactions. In particular we will show that overlapping rarefaction waves are possible and moreover we may have shocks that lie inside a rarefaction wave. In contrast to nonconservative two phase flow models, such as the Baer–Nunziato system, we can use the advantage of the conservative form of the model under consideration. Furthermore, we show the relation between the considered conservative SHTC system and the corresponding barotropic version of the nonconservative Baer–Nunziato model. Additionally, we derive the reduced four equation Kapila system for the case of instantaneous relaxation, which is the common limit system of both, the conservative SHTC model and the non-conservative Baer–Nunziato model. Finally, we compare exact solutions of the Riemann problem with numerical results obtained for the conservative two-phase flow model under consideration, for the non-conservative Baer–Nunziato system and for the Kapila limit. The examples underline the previous analysis of the different wave phenomena, as well as differences and similarities of the three systems.

In this work we study pivoting strategies for the preconditioner presented in Bru (SIAM J Sci Comput 30(5):2302–2318, 2008) which computes the LU factorization of a matrix A. This preconditioner is based on the Inverse Sherman Morrison (ISM) decomposition [Preconditioning sparse nonsymmetric linear systems with the Sherman–Morrison formula. Bru (SIAM J Sci Comput 25(2):701–715, 2003), that using recursion formulas derived from the Sherman-Morrison formula, obtains the direct and inverse LU factors of a matrix. We present a modification of the ISM decomposition that allows for pivoting, and so the computation of preconditioners for any nonsingular matrix. While the ISM algorithm at a given step computes only a new pair of vectors, the new pivoting algorithm in the k-th step also modifies all the remaining vectors from k+1 to n. Thus, it can be seen as a right looking version of the ISM decomposition. The results of numerical experiments with ill-conditioned and highly indefinite matrices arising from different applications show the robustness of the new algorithm, since it is able to solve problems that are not possible to solve otherwise.