Papers

We investigate the semi-discrete Active Flux numerical method on Cartesian grids applied to two-dimensional linear and nonlinear shallow water equations with Coriolis source terms. These hyperbolic systems admit non-trivial stationary solutions governed by a geostrophic equilibrium. In this setting, we analyze the stationarity-preserving properties of our semi-discrete Active Flux formulation. To do so for the linear system, we apply a discrete spatial Fourier transform to the semi-discrete method and analyze the resulting evolution matrix. We show that its discrete kernel is a non-trivial discretization of the geostrophic equilibrium, which implies that the method preserves a discrete geostrophic equilibrium state. This is an inherent feature of Active Flux and is achieved without any additional modifications. Numerical experiments confirm the method's ability to maintain geostrophic stationary states for the linear system. Additionally, we numerically investigate the performance of the Active Flux method for the nonlinear shallow water system.

We propose a time-stepping minimum-residual (MinRes) framework for transient coercive variational problems, applicable to finite-element and neural-network trial approximations. For the Backward Euler scheme, we measure the residual at each time level in the dual norm induced by the steady test-space norm, scaled by the time step. This choice yields stability estimates controlling a discrete parabolic energy quantity, up to the time-discretization defect. We develop conforming and broken-test formulations and specialize the analysis to diffusion-advection-reaction problems. For finite-dimensional test spaces, we derive fully discrete reliability estimates by supplementing the computable discrete residual with a complementary contribution that accounts for residual components the test space does not resolve. We then extend the framework to neural-network trial classes and obtain computable residual decompositions for conforming and broken polynomial test spaces. Numerical experiments confirm the expected finite-element convergence rates, examine the effect of test-space refinement for a global space-time neural approximation, and demonstrate residual-driven spatial refinement for a transient problem with a moving localized feature.

We analyze two fixed-point iterations for computing the principal square root of an M-matrix $A$. Although these iterations, with customary initialization, converge sublinearly when $A$ is a singular M-matrix, we show that, under suitable mild conditions on the initial approximation, the convergence is linear. Moreover, we provide component-wise accurate versions of these iterations, which allow us to approximate the principal square root with a component-wise relative error uniformly bounded by a small multiple of the machine precision. Numerical experiments demonstrating the effectiveness of the proposed algorithms for certain classes of problems are presented.

We investigate the fast-advection limit ($\varepsilon \to 0$) of the stochastically forced two-dimensional incompressible Navier-Stokes equations on the torus. Numerical simulations are performed at resolution $N = 128$ for four different forcing shells on both the square torus and a thin torus, with the advection parameter $\varepsilon$ varied over approximately one decade. We find that the stationary distributions become increasingly concentrated on the lowest Fourier modes as $\varepsilon \to 0$, with the energy-to-enstrophy ratio $E/\Omega$ and the fraction $R$ of enstrophy in the lowest modes both increasing monotonically. At finite $\varepsilon$ the degree of concentration depends on the forcing shell through the effective spectral value $B_1/B_0$, in agreement with the condensation bound of Sznitman and Widmayer. Linear extrapolation to $\varepsilon = 0$ suggests that the limiting values may be forcing-independent, with $B_1/B_0$ governing only the rate of convergence. The same qualitative behaviour persists on the thin torus.

In this paper, we construct axially symmetric kernels on the sphere. We start by taking an isotropic kernel and transform it into an axially symmetric kernel via a specific deformation. For a general deformation, an explicit spherical-harmonic representation is not necessarily available. For our deformation, we can obtain such a representation of the deformed kernel in terms of the expansion coefficients of the original isotropic kernel. We derive this representation by explicitly describing the effect of the deformation on the spherical harmonic basis. From this, an explicit finite sum formula for the coupling coefficients is deduced and we show that, for each fixed input mode, the coefficients decay geometrically at a rate controlled by our deformation parameter. We further note that if the native space of the isotropic kernel is norm equivalent to a Sobolev-space, then the deformed kernel has the same native space with an equivalent norm. In our numerical experiments, we show that our deformation can improve prediction when the underlying random field has an axially symmetric covariance.

We study trust-region minimization of a smooth, possibly nonconvex functional plus a convex functional composed with a bounded linear operator. A relative primal–dual gap condition controls both the error in an approximate proximal-gradient step and its linear-model decrease. Together with a computable absolute stationarity test, it yields a finite Cauchy search, convergence of the proximal stationarity measure to zero, and an $O(\varepsilon^{-2})$ bound on outer trials. The outer analysis allows the linear operator to take values in a Banach space and does not require dual attainment. When the operator takes values in a Hilbert space and the regularizer is finite and Lipschitz, the dual proximal-gradient method produces finite gaps tending to zero, provided the required proximal maps and functional values can be evaluated. We prove $O(j^{-1})$ gap bounds for both recovered and averaged primal candidates and give a sharper bound on the primal error for exactly recovered points. A semilinear elliptic control problem with unsmoothed total-variation regularization and an $L^2$ control cost illustrates the method in the full $H^1$ metric. Across five meshes, outer and state Newton counts remain constant, while interior-point iteration counts vary mildly.

In this work, we propose a hybrid iterative deep Ritz method (H-IDRM) for a class of interface problems for second-order elliptic operators. It is based on a new mixed formulation of the problem and involves solving a sequence of convex minimization problems. We employ a level-set neural network architecture, featuring a level-set representation of the interface, to accommodate the piecewise smoothness of the solution and the flux. The approach involves only volumetric representations instead of duality pairing on the interface and avoids explicit interface sampling that is inconvenient for complex interface geometries. Further, we present an analysis of the method, including the errors arising from the neural network approximation, Monte Carlo approximation, iterative scheme, and penalty parameters. Numerical experiments indicate that the H-IDRM outperforms existing neural solvers on problems with high-dimensional domains, intricate interface geometries, and lower subdomain regularity.

Short-horizon forecasts of atmospheric temperature are needed to support climate-aware digital-twin systems, but such forecasts must be produced where thermal observations are incomplete. This study evaluates a physics-informed neural network for potential-temperature forecasting, constrained by a pressure-coordinate thermodynamic advection-source equation and a diabatic-source closure fit from the preceding 12-hour period and frozen before future-time training. Using hourly ERA5 reanalysis at three pressure levels, the model is evaluated as a conditional hindcast at lead times of one, two and three hours against persistence, local-trend, and two matched neural-network baselines, one of which receives the same future meteorological forcing as the PINN, helping distinguish the physical constraint from access to future forcing. In an Oklahoma development case, mean RMSE improvement over the strongest baseline grew from 8.1% at one hour to 23.8% at three hours; under an observation-density sweep down to 5% of candidate locations, this 3-hour advantage remained 14.6–16.9%, with no evidence that lower density improves performance. Under a fixed protocol transferred to an Alabama heat event with three virtual-observation layouts, three-hour improvement ranged 19.7-24.4% with consistent origin-level wins. A parallel Montana stress test, in which fixed pressure levels intersected complex terrain, produced a three-hour degradation of roughly 17.5%, identifying a terrain-related applicability limit of the formulation. Together, these results indicate that the physics constraint's benefit grows with forecast horizon, persists under severe observation sparsity, and transfers across regions, but is bounded by the validity of a fixed vertical-coordinate representation over complex terrain, evidence relevant to physics-constrained components of climate-aware forecasting and digital-twin systems.

Multicomponent reacting-flow discretizations must conserve total species mass, preserve positive mass fractions, and dissipate a convex mixing potential. These properties are controlled at different stages of a solver, however, and a certificate established at one stage need not survive the next. We study this compatibility problem for mixture-averaged diffusion, non-orthogonal finite-volume reconstruction, and positivity-preserving Patankar time integration. The analysis is isothermal, uses prescribed density, and uses synthetic detailed-balance reaction networks for the chemistry tests. At the face level, the entropy-variable jump is shown to equal a log-mean Gibbs Hessian applied exactly to the mass-fraction jump. Requiring conservation and entropy neutrality of the mass-correction term on the simplex tangent space then yields a unique componentwise-positive face composition and an exact two-point entropy certificate. The arithmetic-weight face criterion is also satisfied by production neutral-air/carbon data and a quasi-neutral ambipolar-air control, but a face-level certificate does not automatically extend to subsequent discretization layers. A weighted least-squares non-orthogonal correction can, in the tested prototype, convert a dissipative two-point operator into an anti-dissipative multipoint operator. Likewise, Patankar rescaling can preserve positivity and conservation while increasing the free energy; an explicit finite-step rational counterexample establishes this failure independently of ensemble statistics. A fully discrete two-cell example further shows that a positive face entropy margin does not guarantee positivity of the fully discrete update. Thus conservation, positivity, and free-energy dissipation are layer-specific properties and require distinct discrete certificates.

Variational data assimilation computes a Bayesian update through optimization on a fixed posterior geometry, whose conditioning is determined by the prior covariance and the observation operator. In strongly anisotropic, weakly observed, or highly informative regimes, this geometry can become severely ill conditioned, leading to substantial sensitivity to solver design and preconditioning. This work investigates an alternative computational route based on flow-based transport, using the ensemble score filter as a training-free realization. Rather than repeatedly optimizing over the full posterior geometry, the flow-based update transports the predictive distribution through a sequence of intermediate distributions whose covariance and curvature are progressively regularized before reaching the target posterior. In the linear-Gaussian setting, we characterize this mechanism through the evolution of the diffused posterior covariance and the associated time-dependent curvature, revealing a continuation from a near-isotropic reference geometry to the ill-conditioned posterior geometry targeted by variational assimilation. Numerical stress tests show that this geometric regularization makes the flow-based method substantially less sensitive to severe ill-conditioning than the variational method, while a heterogeneous $10^4$-dimensional Lorenz-96 benchmark demonstrates improved robustness and lower delivered computational cost. These results identify posterior-geometry regularization as a key mechanism distinguishing transport-based from optimization-based data assimilation.

Topology optimization (TO) remains computationally intensive due to repeated finite element analysis (FEA) evaluations required at each iteration. While neural network-based surrogates offer potential acceleration, existing approaches often suffer from gradient inconsistency between predicted objectives and sensitivities, leading to optimization instability. This work presents KATOsuper, an objective-agnostic framework that couples neural-reparameterized topology optimization with a Sensitivity-Consistent Fourier Neural Operator (SC-FNO). The framework employs the forward_split architecture, which derives deployed sensitivities via automatic differentiation through the predicted objective field and thereby preserves consistency between the predicted objective and the gradient used for optimization. The case studies include three 2D benchmark problems and three 3D structures considering compliance or stress minimization. A physics-informed multi-channel input encoding with Fourier position embedding enables resolution-invariant learning, supporting zero-shot extrapolation beyond the training resolution, with useful performance at moderate scaling factors and topology-preserving exploration at up to 64x without retraining. The framework extends to 3D through KATO3D, featuring novel KANConv3D blocks with learnable B-spline activations. KATOsuper demonstrates 15–110x deployment-time speedup over MATLAB baselines while maintaining competitive optimality, with the clearest gains observed in complex 3D and stress-optimization cases. The insight that sensitivity direction matters more than magnitude enables robust optimization even with approximate physics evaluation, extensible to other differentiable physics-driven design objectives.

We systematically investigate breather solutions in the modified Korteweg-de Vries (mKdV) equation via numerical simulations. We show that the breather solutions are stable under a variety of perturbations, including amplitude changes, modifications of internal parameters, and perturbations of the nonlinearity of the equation. Our results are not only consistent with the known analytical studies on the stability of breathers in integrable systems, extending them further, but also provide numerical evidence that such breather-type structures can persist and remain stable over the simulated time scales in some non-integrable settings with symmetric potentials.

We study the reconstruction of a directional albedo field from confocal non-line-of-sight measurements. For mirror-symmetric vector fields in Sobolev spaces, radial preprocessing reduces the data to spherical means of the divergence with centers on the relay wall. Full data determine this divergence uniquely, with divergence-free fields forming the entire ambiguity. An explicit Fourier–sine formula recovers the irrotational Helmholtz component and reproduces the data. The Fourier–sine transforms of the model data form a weighted Hilbert space that we characterize exactly, with norm equal to the Sobolev norm of the reconstructed field. This weighted range remains well defined even when the preprocessed data fail to belong to standard Sobolev spaces. For Schwartz fields in the model class, such Sobolev regularity holds exactly when the divergence has zero depth integral. Measurements on an open subset of the relay wall, for all radii $0<r<R_{\max}$, uniquely determine the divergence in the union of the corresponding balls. An FFT-based algorithm with Stolt interpolation implements the reconstruction in $O(N^3\log N)$ operations on an $N^3$ grid. Tests on explicitly defined fields assess reconstruction accuracy and invariance under divergence-free perturbations. We present imaging examples of potential and vector-field reconstruction from model-generated, rendered, and measured transients.

Residual-based adaptive GMsFEM for high-contrast elliptic problems repeatedly evaluates local weighted $H^{-1}$ indicators on every coarse neighborhood, making indicator evaluation a recurring cost in repeated-query settings. We introduce a non-intrusive Gaussian-process (GP) surrogate for the indicator scores used in Dörfler marking. The operational predictor is the GP posterior mean, algebraically equivalent to a kernel ridge regression (KRR) estimator under the stated convention; it uses compressed local solution and spectral features without changing the multiscale solve, local spectral construction, or basis enrichment. A nonuniform perturbed-marking result quantifies how pointwise score errors affect the exact indicator mass captured by surrogate-selected neighborhoods, while a conditional bounded-discrepancy KRR pathway identifies sufficient assumptions for such score bounds. In controlled held-out in-distribution tests, the surrogate-guided method gives error-versus-DoF trends comparable with classical $H^{-1}$-residual offline adaptivity and evaluates the online indicator component 2.0-2.1 times faster, excluding offline data generation and GP training.

We establish a fully discrete weak-to-strong (W2S) multi-entropy stability theory for arbitrarily high-order finite-volume and discontinuous Galerkin (DG) approximations of the ideal compressible magnetohydrodynamics (MHD) equations on general polytopal meshes, whereby a single numerical update simultaneously satisfies discrete entropy inequalities for any prescribed finite family of convex Harten entropy pairs. Because physical entropies are defined only for positive density and pressure, the central analytical difficulty lies in reconciling discrete entropy stability with positivity preservation, the magnetic divergence constraint, and the nonconservative Godunov–Powell coupling. Building upon the provably positivity-preserving MHD framework of Wu and Shu, we develop a weak multi-entropy analysis of the common cell-average evolution of finite-volume and DG methods. Through convex decomposition and relative entropy, we establish two sufficient criteria for this weak stability general polytopal meshes. Compatible magnetic corrections and locally divergence-free approximations cooperatively offset the nonconservative Godunov–Powell coupling and provide the discrete cancellations essential to both positivity-preservation and entropy stability. A W2S lifting then unifies positivity and entropy limiting into a single cellwise scaling limiter, achieving strong multi-entropy stability while strictly preserving cell averages and the locally divergence-free magnetic field. High-order temporal accuracy follows from variable-step strong-stability-preserving multistep methods. Numerical experiments confirm the theoretical findings and demonstrate computational robustness. This framework provides a rigorous, fully discrete foundation for unifying physical admissibility, magnetic divergence control, and multi-entropy stability in high-order MHD approximations.

Row-sparse projection provides a useful tool in machine learning (ML) when it comes to, for example, feature selection, aiming to choose most relevant features for various ML objectives. One way to seek a high quality row-sparse projection is to combine an ML objective, such as the ones for PCA, LDA, and OCCA, with the matrix $(2,1)$-norm regularization which is nonsmooth. Such combinations result in challenging optimization problems on the Stiefel manifold that need to be solved efficiently. In this paper, a unifying NEPv framework is established to efficiently deal with optimization on the Stiefel manifold with the $(2,1)$-norm regularization. The effect of the $(2,1)$-norm regularization is also investigated. The wide applicability of the framework is demonstrated through the combinations of common learning objectives in today's data science applications with the $(2,1)$-norm regularization. Numerical experiments are presented to illustrate the use of the NEPv approach and to gain insights as to what a proper regularizing parameter should have in real-world applications.

We present a distributed framework for large-scale Bayesian inverse problems governed by the eikonal equation, with a specific focus on seismic traveltime tomography. Traditional deterministic approaches often fail to provide the uncertainty quantification (UQ) necessary for ill-posed problems, while conventional Bayesian sampling methods such as Markov chain Monte Carlo (MCMC) suffer from the curse of dimensionality and slow convergence in high-dimensional model spaces. The proposed framework addresses these challenges through a three-tier computational strategy. First, we utilize the Fast Marching Method (FMM) to solve the eikonal equation, ensuring high numerical accuracy. Second, we reformulate the global tomographic objective into a decentralized consensus form, allowing the inversion to be decomposed into independent subproblems solved in parallel via the Alternating Direction Method of Multipliers (ADMM). This architecture eliminates the need for the explicit construction of large-scale sensitivity matrices, significantly reducing the memory footprint for 3D surveys. Finally, we integrate Stein Variational Gradient Descent (SVGD) within the ADMM workers to perform approximate posterior sampling. By evolving a set of model particles along a functional gradient direction that balances data-fitting forces with a repulsive kernel-based diversity force, we obtain an ensemble from which posterior summaries are computed. We derive a data-space Gauss-Newton update using the Woodbury matrix identity to further accelerate the particle evolution in large-scale 3D problems. Numerical experiments on complex 2D and 3D models demonstrate that the algorithm achieves stable convergence, produces high-fidelity velocity reconstructions, and provides posterior uncertainty maps.

In sparse random matrices, spectral outliers (eigenvalues and singular values located away from the bulk) emerge due to degree fluctuations: high degrees inflate the operator norm, while low column degrees reduce the least singular value. As proved by Feige and Ofek (2005) and Le, Levina, and Vershynin (2017), degree regularization enforces concentration at the expected norm scale. However, precise bounds incorporating the cutoffs remain unexplored and challenging since regularization introduces dependencies among entries. For the first time in the literature, we provide variance- and cutoff-dependent bounds for extreme singular values and eigenvalues of regularized inhomogeneous random matrices. In the absence of regularization, our lower bound for the least singular value matches the same leading constant obtained by Brailovskaya and van Handel (2024). Moreover, our error term vanishes under the milder condition $d/\log N\to\infty$, as opposed to their stronger requirement $d/(\log N)^4\to\infty$. A key ingredient is to extend spectral radius bounds for nonbacktracking matrices to the dependent setting. We build on approaches for independent cases established by Benaych-Georges, Bordenave, and Knowles (2020), as well as Dumitriu and Zhu (2024), and carefully handle edges traversed only once. Our proof framework separates deterministic spectral comparisons from probabilistic estimates: once Loewner inequalities and columnwise variance controls are established, the remaining probabilistic analysis boils down to verifying the graph moment conditions formulated in this paper. We hope this framework can be extended to handle general random matrices with more complex dependencies.

These notes develop the mathematical foundations and construction of a Fourier-Bessel wavelet family inspired by the disk harmonics of Shaqfa et al.[9]. We begin with the relevant properties of Bessel and modified Bessel functions and introduce the wavelet properties required for the construction. We then derive the Fourier-Bessel disk harmonics as solutions to the Helmholtz equation on the unit disk subject to a Neumann boundary condition. Building on this basis, we construct a wavelet family by applying a Gaussian spatial envelope and introducing a zero-mean correction for the zeroth angular order. We derive the corresponding normalisation constants for $L^2$-based applications and discuss $L^1$-based normalisation for frequency-domain peak consistency. Finally, we derive a closed-form Fourier-domain representation of the resulting wavelets. The main motivation is the approximately linear spacing, which converges to $\pi$ between consecutive radial eigenvalues. Rather than replacing the conventional dyadic organisation of wavelet families, this construction lays out the foundation to explore whether a more uniform radial frequency allocation can be useful for applications in which broad and balanced frequency coverage is desirable.

This paper presents a unified error analysis of randomized neural networks (RaNNs) for elliptic interface problems. We first derive an integral representation of Barron functions using the standard $\tanh$ activation, which yields a $C^2$-approximation bound valid for any hidden-parameter distribution supported on a cube with strictly positive density. We then establish a quantitative $\mathcal{O}(n^{-1})$ generalization bound for the composite loss that enforces the partial differential equation, boundary, and interface conditions simultaneously, where $n$ denotes the number of training samples. Combining these results with the optimization guarantees for RaNNs, we obtain an a priori $L^2$ error estimate for the neural network solution in the presence of interface discontinuities. Numerical experiments are presented to validate the theoretical results.

We introduce a new formulation of the Longest Edge Bisection (LEB) of tetrahedra entirely in sextuple space R6, where tetrahedra are represented by the squares of their edge lengths. This representation renders the LEB refinement equations fully linear and eliminates the need for coordinate-based data structures. A central difficulty in three-dimensional LEB arises when a tetrahedron possesses multiple longest edges, making the refinement rule intrinsically multivalued. We formalize this phenomenon through the notion of Multiform Longest Edge Bisection (MLEB), which systematically explores all admissible longest-edge choices. To encode this multivalued behavior, we introduce the concept of bisection patterns, defined as sequences of sextuple permutations governing the refinement process. We prove that the set of sextuples sharing a common LEB pattern forms a convex region in R6. For structurally significant families of tetrahedra, including the R1+ family and the Liu-Joe family, we show that the infinite refinement tree collapses into a finite directed graph with eight states. Remarkably, both families are governed by the same graph, differing only in their initial state. This directed-graph formulation provides a unified combinatorial description of the refinement process and offers an efficient computational framework for deep iterative LEB analysis.

In this work, we present a linear Robin-type domain decomposition scheme for solving the mixed formulation of Richards' equation (mLRDD-scheme), governing flow in variably saturated porous media. Assuming a highly heterogeneous porous medium consisting of multiple blocks/layers of distinct materials, we apply non-overlapping domain decomposition to isolate individual materials, yielding near-homogeneous conditions within each subdomain. Richards' equation is discretized in time by backward Euler and linearized using the L-scheme. A Robin-type interface condition, consistent with the mixed formulation, is then derived to couple the subdomains. For spatial discretization, mixed finite elements are employed, ensuring local mass conservation and providing the flux variables consistent with the Robin condition. The convergence of the proposed scheme is rigorously proved. Heterogeneous, multidomain numerical examples in 2D/3D are presented to demonstrate the efficiency and robustness of the scheme, along with numerical validation of the theoretical convergence results.

Gaussian random fields are a versatile tool used in the fields of stochastic PDEs, uncertainty quantification, and geostatistical simulation. One way to obtain them is to use a truncated Karhunen-Loève expansion (KLE). Computing the expansion requires the leading eigenpairs of an $N \times N$ covariance matrix, where $N$ is the total number of grid cells. These are usually computed with a Krylov eigensolver, which relies on the covariance operator only within matrix-vector products. Stored densely, the matrix takes $\mathcal{O}(N^{2})$ memory and each product $\mathcal{O}(N^{2})$ time. For a stationary kernel on an equispaced grid, the covariance matrix becomes (block-) Toeplitz and the product evaluates in $\mathcal{O}(N\log N)$ time using FFT-based circulant embedding, without the need to assemble the dense matrix. In a matched single-threaded comparison, the median speedup of the eigensolve grows from $18 \times$ at $N = 4096$ to $183 \times$ at $N = 2^{15}$. This makes it possible to compute discretized fields that would otherwise be infeasible to compute in the standard formulation. We show that the same construction carries over to non-separable kernels as well as $d$ dimensions, using block-Toeplitz matrices. We extend it to piecewise-constant fields on arbitrary domains given as subsets of a tensor grid. Computational savings grow with problem size, and storage drops from $\mathcal{O}(N^{2})$ to $\mathcal{O}(2^{d}N)$.

Drifting provides a direct route to one-step generative models, but applying it directly to stochastic transition modeling requires multiple samples of the next state conditioned on the same current state. Standard trajectory data, however, typically provide only one realized next state for each observed current state and therefore do not provide an empirical approximation of the corresponding conditional distribution over possible next states. We introduce block-triangular joint drifting, which instead applies a projected drift field to the empirically accessible joint distribution of consecutive states. Importantly, the block-triangular architecture preserves the current-state marginal while making its second component a direct sampler of the conditional distribution of possible next states. The resulting surrogate generates stochastic trajectories with one model evaluation per time step, without auxiliary generative steps between time steps. Numerical experiments demonstrate accurate marginal and trajectory-dependent statistics and favorable accuracy-cost tradeoffs compared with deterministic, diffusion-, flow-, and distillation-based generative surrogate models.

Principal component analysis (PCA) is a fundamental tool to reduce the dimensionality of the data in many applications. PCA finds a few signal directions that contain most of the variability of the data by computing the eigenvectors of the sample covariance matrix. In this work, we focus on the spiked covariance model, in which the data vectors are defined by a few orthogonal signals plus an isotropic Gaussian noise, and our goal is to estimate one or more of the leading signals. Our main theoretical finding is that the subspace spanned by several leading eigenvectors of the sample covariance matrix contains significant information about the desired signals long before the individual eigenvectors converge to the population principal components. To prove this, we derive a posteriori bounds for the angle between the subspace spanned by the desired population signals and the subspace obtained from the sample using perturbation theory for singular vectors. This leads to a new algorithm, SuperPCA (SUbsPace subsamplER PCA), which capitalizes on an approximate eigenspace of the sample covariance matrix to find the leading signals far more efficiently and accurately than classical PCA in the high-dimensional, multi-signal setting. SuperPCA exploits only a small number of subsampled coordinates of the data, which can lead to tremendous savings in data acquisition cost, especially when the signals are approximately sparse. For the same number of measurements, SuperPCA can offer a factor $10$ improvement in accuracy compared to the classical PCA method.

We develop a new endgame for numerically computing singular solutions to complex analytic systems by means of a homotopy path whose endpoint is the solution to be computed. When the solution is singular, prediction-correction path tracking may fail at the endpoint. Instead of tracking the path all the way to its end, endgames sample nonsingular points on the incoming path and approximate the endpoint by interpolating the samples with a truncated Puiseux series. By collecting sample points along a logarithmic spiral in the complex plane, the proposed spiral endgame balances between the superior numerical conditioning of a circular sample (known as the Cauchy endgame) and the update efficiency of a linear sample (known as the power series endgame), these being tied together as two extremes of the new approach. The numerical conditioning of the methods are analyzed and compared, and computational examples illustrate the effectiveness of the new spiral sampling method compared with the circular sampling method of the Cauchy endgame.

Compliant arterial wall mechanics, contributing to the Windkessel effect in the carotid artery (CA), have significant impact on hemodynamic patterns and surface shear indicators for cardiovascular diseases, e.g. atherosclerosis. In this study, we extend the linear elastic Fluid-Structure Interaction (FSI) framework to account for nonlinear strain-dependent wall behavior. The model incorporates a Young's modulus generalized from tensile tests on silicone phantoms to capture the nonlinear stress-strain relation of arterial tissue. The final computational model, using a resistance-type boundary condition and in vitro measured stress-strain relation is validated against both, in vitro assessed silicon CA phantom as well as published clinical data on the flow splitting to the daughter branches in CA bifurcations. To better reflect physiological conditions, our model is subsequently extended to incorporate prestress of patient-specific geometries, and clinically measured stress-strain relations, followed by validation against clinical CA data. The present study demonstrates the feasibility of strain-dependent Young's (elastic) modulus as a means to enhance the capacity of the linear elastic framework to accurately represent the physiologically nonlinear mechanics of arterial walls, striking a balance between implementation effort and physiological fidelity. This approach yields realistic strains, volumetric inflation behavior and nonlinear pressure-volume relationships. Furthermore, the study reveals the importance of proper outlet boundary conditions and the shortcomings of resistance-type boundary condition in patient-specific modeling, leading to non-physiological pressure profiles and non-physiological temporal flow splitting.

The increasing complexity in design and manufacturing has driven the need for advanced techniques for fast modeling problems with large-dimensional parameter spaces. Avoiding high-fidelity finite element models while achieving fast and accurate simulations in such contexts is challenging. Parametric model order reduction (pMOR) has drawn significant attention in recent years. Nevertheless, the curse of dimensionality in parameter spaces has severely limited its effectiveness. The active subspace (AS) approach has been successfully applied to pMOR for systems with many parameters. However, the balance between accuracy and compactness of the reduced model remains problematic for such systems with high-dimensional parameter spaces. It often results in models that are either not small enough or not accurate enough. In this paper, we propose an iterative active subspace (IAS) approach for parametric model order reduction, which, to some extent, addresses the trade-off between accuracy and reduced model size and achieves substantial computational gains compared to the original active subspace method.

In this paper, an abstract framework is introduced for the derivation of reliable a posteriori error estimators for variational inequalities of the second kind. Under certain conditions (local) efficiency of the error estimators can also be derived. The framework is based on the estimation of a residuum which results from an equivalent mixed formulation and from the replacement of the Lagrange multiplier by its (possibly post-processed) discrete counterpart. The approach is applied to an idealized frictional contact problem and a Mosolov type model problem in an $hp$-finite element discretization setup. Several numerical experiments demonstrate the broad applicability of the theoretical findings.

We present a new semi-implicit structure-preserving (SP) finite volume discretization for a unified first-order hyperbolic model of continuum mechanics. This framework provides a common mathematical description for fluids and solids and incorporates transport, viscous effects, heat conduction, and elastic deformations within a single system of hyperbolic partial differential equations. The coexistence of several physical mechanisms gives rise to multiple characteristic wave speeds, resulting in severe time-step restrictions for fully explicit discretizations. To overcome this difficulty, we develop a four-split scheme in which the governing equations are decomposed into convective, temperature, mechanical, and pressure subsystems. The convective subsystem is the only one that is advanced explicitly in time, while the remaining subsystems are treated implicitly in a sequential manner. As a consequence, the time step restriction of the proposed method depends only on the material velocity and is independent of the acoustic, shear, and thermal wave speeds that characterize the model. The scheme is designed to preserve key structural properties of the underlying equations on unstructured grids. In particular, a compatible vertex-staggered spatial discretization on triangles allows the curl-free involutions associated with the distortion field and thermal impulse to be respected whenever the relaxation source terms in the governing PDE are linear or absent. At the same time, the scheme is asymptotic preserving as it is consistent with the low Mach number limit and with the stiff relaxation limits of the GPR model that recover the classical Navier-Stokes-Fourier equations in fluid mechanics. A series of numerical experiments covering both fluids and solids demonstrates the accuracy, robustness, and multi-scale capabilities of the proposed approach.