2026-09-15 | | Total: 47
We introduce a new traffic flow model. The main idea is to describe acceleration effects in the Aw-Rascle-Zhang model through an additional driver-dependent equation, thereby allowing individual driver characteristics to evolve dynamically rather than being prescribed solely by the traffic density. The resulting system provides a hyperbolic description of traffic dynamics and admits an explicit analysis of its characteristic structure, Riemann problem, and entropy properties. To investigate the stability of the proposed model, we perform a Chapman-Enskog expansion up to second order and derive the corresponding diffusive and higher-order corrections. This analysis reveals how the stability properties depend on the underlying velocity and pressure functions. Several examples are presented to illustrate the theoretical results.
We study the temporal error growth of the Strang splitting method for the periodic cubic nonlinear Schrödinger (NLS) equation. The leading error is governed by a forced linearized equation, whose growth depends sharply on dimension and the sign of the nonlinearity. In the 1D defocusing case, we prove a uniform quadratic upper bound $C(1+T^2)τ^2$ using the global Birkhoff transformation and the resulting degenerate structure of the linearized flow. In the higher-dimensional defocusing case, the exponential growth is constructed using arbitrarily small unstable standing waves. Moreover, to prove the higher-dimensional defocusing upper bound, we use the periodic Strichartz estimates of Killip and Vişan to show that $\int_0^T\|u(t)\|_{L^\infty}^2\,dt\le C(u_0)(1+T)$. This yields an exponential rate independent of the final time, despite possible growth of higher Sobolev norms. The exponential lower bound in the focusing case is obtained from the unstable linearized dynamics around a plane wave, with the Akhmediev breather providing the underlying mechanism. In addition, for 1D focusing case with small initial data, the appliction of the local Birkhoff transformation ensures that the error grows at most quadratically in time. Various numerical experiments confirm the sharp linear, quadratic, and exponential error growth rates.
The finite difference (FD) method is commonly used to approximate derivatives of smooth functions, with accuracy typically determined by stencil size and derivative order. However, certain centered stencils exhibit unexpectedly higher accuracy, a phenomenon known as superconvergence, which has been observed in practice but lacks rigorous explanation. We present a mathematical framework for superconvergence in centered FD approximations based on Taylor expansions of the truncation error and the resulting linear system for the FD coefficients. By analyzing symmetry properties of these coefficients and their interaction with the parity of the derivative order, we identify conditions under which higher-order error terms cancel. We show that superconvergence occurs for odd-order derivatives with even centered stencils and for even-order derivatives with odd centered stencils, while no superconvergence occurs for even derivatives with even centered stencils. Numerical experiments in MATLAB confirm the predicted convergence rates.
We present a new perspective on integral operator kernels given by a nonlinear function of the distance. This perspective is useful since it permits one to quasilinearize the nonlinear transformations acting on the tensor Haar expansion of the distance, d(x,y). This leads to a new representation comprised of a principal and residual term with desirable features; namely, enhanced local regularity of the principal component (and in certain situations rapid off-diagonal decay leading to O(N) storage) and a residual expansion of improved approximation as a consequence of enhanced global regularity. Numerical experiments on the potential kernel are included
Nested cubature rules, in which each refinement retains all previously used nodes and thus reuses earlier function evaluations, are natural in multilevel and adaptive integration. For every fixed $s>d/2$, we show that the equal-weight QMC integration rate on $\mathbb S^d$ is compatible with such nested point sets. In the subcritical range $d/2<s<d$, cumulative unions of geometrically growing QMC blocks yield nested QMC design sequences whose successive cardinality ratios converge to any prescribed $ρ>1$. At and above the critical index $s=d$, where the block-averaging estimate no longer yields the optimal rate, we prove an equal-weight completion theorem based on low-frequency discrepancy cancellation. A block-sensitive estimate sharpens the iteration and yields nested QMC designs for all $s\ge d$ with $N_{j+1}\lesssim(j+1)^{2s/d-1}N_j$. Every prescribed finite point set also admits an optimal-rate completion.
In this work, we propose an unfitted finite element scheme to approximate the solution of the heat equation on moving domains. We use the φ-FEM paradigm in which the computational domain is described implicitly by a level-set function φ. This function is incorporated into the variational formulation in order to enforce the boundary conditions, which allows the use of unfitted meshes in space. Such a strategy avoids the need for remeshing at each time step and makes it possible to handle complex geometrical evolutions. Moreover, the φ-FEM approach has the advantage of being simple to implement within standard finite element libraries. We introduce a fully discrete scheme that combines the φ-FEM spatial discretization with the lowest-order discontinuous Galerkin method in time. Under regularity assumptions on the level-set function and a mild restriction linking the time step to the mesh size, we establish an optimal a priori error estimate in the L2(0,T;H1) norm. Finally, we present several numerical experiments that confirm this convergence rate, exhibit a second-order convergence in the L{\infty}(0,T;L2) norm, and illustrate the robustness and accuracy of the proposed method.
Asynchronous iterative methods are attractive for large-scale parallel computing because they reduce synchronization and communication overhead. Existing convergence-rate analyses, however, have primarily focused on shared memory implementations, whereas distributed memory systems introduce more general and potentially inconsistent communication delays. In this work, we revisit asynchronous Jacobi and randomized Gauss--Seidel (RGS) methods for symmetric positive definite linear systems from a unified perspective. We first introduce a general asynchronous model that encompasses both shared memory and distributed memory settings and expresses the two methods through a common coordinate-update framework. We then establish linear convergence in expectation under an explicit stability condition. The resulting convergence bound depends algebraically on the delay through the quantity $\sqrt{ρτ}+ρτ$, where $τ$ is the maximum communication delay and $ρ$ reflects the communication pattern of the underlying parallel implementation. In particular, under an appropriate scaling regime with $ρτ=O(1)$, the guaranteed per-iteration convergence rate has the same asymptotic order as that of synchronous RGS. These results provide a unified framework for quantifying the effect of asynchronicity on asynchronous Jacobi/RGS methods in both shared and distributed memory environments.
The fast Fourier transform (FFT) and its real counterparts, the fast algorithms of the discrete cosine transform (DCT) and the discrete sine transform (DST), are frequently employed for the numerical evaluation of the Fourier cosine and Fourier sine transform, thereby acting as parametric quadrature rules based on finitely many samples of a given continuous function. But how accurate are the obtained results? In this paper we study the error occurring if the DCT of type I (DCT-I) and the DST of type I (DST-I) are applied to compute $\int_0^{\infty} f(x)\cos(2πxv)\d x$ and $\int_0^{\infty} f(x)\sin(2πxv)\d x$ for $f \in L^1(\Rp) \cap C(\Rp)$ and $v \in \Rp$, as well as the error occurring if the Chebyshev coefficients $\frac{2}π\int_0^π h(\cosθ)\,\cos(nθ)\dθ$, $n \in \Np$, of a function $h \in C(I)$ on $I := [-1,1]$ are computed by the DCT-I. We present practicable, explicit error estimates under polynomial, exponential, and mixed decay conditions. For the Fourier cosine and sine transform, rational functions lead to a geometric decay in the frequency parameter $P$ combined with an algebraic plateau of order $L^{-3/2}$ in the sampling parameter $L$, while for meromorphic functions such as $\operatorname{sech}(π\,\cdot)$ the critical geometric rate is attained in both parameters. For Chebyshev approximation of a rational function with simple poles in $\C \setminus I$, the quadrature error has the critical geometric rate $ρ_0^{-P}$ with an explicit constant, where $ρ_0 > 1$ is the parameter of the largest Bernstein ellipse of analyticity. Several numerical examples illustrate the theory.
We consider partial differential equation (PDE) inverse problems from sparse, noisy, and corrupted observations, with the aim of recovering unknown coefficient fields and associated state variables in a mesh-free setting. Under such sparse and corrupted observations, standard physics-informed and generative approaches typically treat all samples indiscriminately and therefore lack a principled mechanism for reconciling physical laws with contaminated data. We address this difficulty with a two-stage flow-matching framework. In the first stage, we develop physics-guided conditional flow matching (PG-CFM), which incorporates strong-form PDE information through residual regularization along the generative trajectories together with intermittent global collocation constraints. In the second stage, we introduce energy-regularized flow matching (ERFM), which fine-tunes the Stage-1 model by assigning each observation a physics--data energy score from a frozen teacher and reweighting the flow-matching objective to reduce the influence of high-energy, PDE-inconsistent samples. We show that the resulting Stage-2 objective is equivalent to flow matching under a teacher-induced reweighted data distribution, which gives a population-level interpretation of the robustness mechanism. Numerical experiments on several inverse benchmarks, including Poisson and Navier--Stokes problems, show that the proposed framework yields more accurate coefficient recovery than robust PINN variants and competing generative baselines.
We propose a mixed finite element discretisation for the incompressible Navier-Stokes equations that preserves the evolution laws of both energy and enstrophy, in a stronger sense than previous discretisations. In two dimensions, the evolution law for enstrophy only permits dissipation for thermodynamically isolated systems, leading to a Reynolds-number-independent bound on the velocity gradient that naturally stabilises the scheme, even on severely under-resolved meshes. In three dimensions, the scheme preserves both dissipation and the generation of enstrophy through vortex stretching. We enforce these evolution laws by systematically introducing auxiliary variables into the discretisation. While conforming implementations of these schemes require discrete Stokes complexes with enhanced regularity, we introduce both (i) equivalent reparametrisations and (ii) penalty formulations that require only the typical curl- and div-conforming spaces from the standard discrete de Rham complex. The scheme handles different types of boundary conditions and curved domains. The robust stabilisation properties of the proposed scheme are demonstrated through numerical simulations of a shear flow, a spherical vortex, and flow past an obstacle. We observe numerically that preserving the discrete evolution of enstrophy in this way has a strong stabilising effect on the numerical solution, especially in two dimensions.
We propose a dissipative, irreversible damage model for anisotropic media in large deformations to address the damage process of the myocardium following a cardiac infarction. We model damage to the collagen in the cleavage planes resulting from an increased load to the passive tissue. Starting from variational principles, we derive a quasi-static differential model governing the irreversible evolution of the damage. We apply our model to two test cases. First, a rectangular passive slab geometry to which an indenter-like load is applied. Then, we couple our model with a model for the numerical simulation of left ventricular electromechanics. Our results show that a redistribution of the load within a left ventricle occurring following a cardiac infarction can produce the damage to the collagen scaffold in the infarcted area consistent with experimental results. This is indicative of the fact that passive load-dependent dissipative phenomena are implicated in post-infarction remodelling.
In this paper we propose augmented Lagrangian finite element methods for cavitation in the Reynolds and Stokes models of lubrication, for which the discrete complementarity conditions hold exactly and the augmentation parameter drops out of the method. For the Reynolds equation this follows from nodal quadrature in a mixed piecewise linear method, for which we prove the classical first-order error estimate. We also determine the computable stability threshold of a multiplier-free stabilised alternative. For Stokes flow we show that the constrained scalar is the mechanical pressure when the deviatoric stress, with zero bulk viscosity, is used, but not with the customary incompressible stress. We discretise with a jump-stabilised Crouzeix-Raviart element and piecewise constant pressure, and prove well-posedness, stability in two and three dimensions, and a first-order error estimate. We present numerical examples which verify these properties and show that close pressure profiles in the two models do not imply close cavity predictions, the latter also being sensitive to the end conditions of the computational domain.
We propose a hybridized staggered discontinuous Galerkin--mixed finite element method for the strain gradient elasticity model, a fourth-order singularly perturbed problem governed by a material length scale parameter $ι$ and the Lamé constants $λ$ and $μ$. Introducing the total stress together with the scaled Cauchy and hyper stresses as auxiliary unknowns, we recast the model as a first-order system and discretize the displacement and the total stress by staggered discontinuous Galerkin spaces, the two scaled stresses by Raviart--Thomas pairs, with symmetry imposed strongly on all three stresses. The higher-order Dirichlet boundary condition enters the variational formulation naturally, so that the scheme avoids the numerical boundary layer that limits displacement-based methods. By establishing an inf-sup condition on the symmetric subspace and a discrete Korn inequality, we prove algebraic stability and optimal convergence, both uniform in $ι$ and $λ$. We further develop a hybridized scheme in which local static condensation leaves only two multipliers on the mesh skeleton as globally coupled unknowns. Numerical experiments confirm the predicted convergence rates and the parameter robustness.
We present shape calculus techniques and a scalable automatic shape differentiation framework for multi-phase topology optimisation on unfitted discretisations defined by several level-set functions. First we establish general, exact shape calculus expressions in the discrete case for multi-phase systems by leveraging concepts from convex geometry. To complement this theoretical foundation, we introduce an open-source multi-phase automatic shape differentiation framework based on polytopal cutting. This computational framework is validated against both finite differences and our established exact expressions, matching the latter to near machine precision. Furthermore, the proposed automatic shape differentiation is scalable across distributed computing environments, demonstrating near ideal weak scaling up to 1.65 billion finite elements across 13,824 computer cores. We demonstrate our implementation by solving unfitted multi-phase topology optimisation problems for anisotropic diffusion, linear elasticity, and fluid--structure interaction. Together, these theoretical and computational contributions provide a robust and accessible foundation for advancing multi-phase topology optimisation using unfitted finite element methods. In particular, the methods enable the solution of topology optimisation problems involving multi-phase and multi-physics systems with non-trivial boundary conditions. The open-source software is available at https://github.com/zjwegert/GridapTopOpt.jl.
First-arrival traveltime tomography provides a computationally efficient means of imaging subsurface velocity structure, while boundary-conforming deformed grids allow rugged topography to be represented without abandoning a logically structured mesh. In matrix-free Gauss-Newton inversion, however, the linearized full-metric Eikonal transport and its transpose must be applied repeatedly at a fixed background model. Directional sweeping consequently revisits the same state-dependent dependency structure for every Krylov right-hand side. We recast the frozen tangent transport as a directed graph and distinguish loss of the background-traveltime ordering from genuine algebraic cyclicity. Strongly connected components identify the irreducible part of the transport, whereas the remaining dependencies admit exact scalar substitution after reordering. The condensation graph and local block factorizations are constructed once for each frozen source state and then reused for both primal and transpose applications. Numerical experiments on deformed grids and a three-dimensional tomography problem show that the resulting block-triangular formulation preserves the frozen sensitivity action while substantially reducing the repeated transport cost. By eliminating this dominant inner-iteration cost, the proposed method markedly accelerates Gauss-Newton computation and makes full-metric matrix-free inversion on deformed grids practical at much lower cost.
Simulating elastic wave propagation in heterogeneous media presents significant mathematical and computational challenges, since resolving fine-scale material variations requires extremely fine spatial discretizations, while mixed stress--displacement formulations typically lead to large saddle-point systems that are not well suited for efficient time integration. Moreover, conventional multiscale reductions generally produce nontrivial coarse mass matrices, thereby requiring additional linear solves at every time step and diminishing the advantages of explicit schemes. To address these difficulties, we develop a time explicit multiscale method based on a multipoint stress control volume discretization. The central innovation is a unified spatial--temporal reduction strategy that simultaneously removes two dominant algebraic bottlenecks: the complex global fine-scale mixed solve and the repeated coarse-scale mass solve arising in conventional multiscale time stepping. More specifically, the stress and auxiliary rotation variables are eliminated locally, leading to a symmetric positive definite reduced stiffness operator, upon which a multiscale space is constructed using a density-weighted local spectral problem. The associated projected mass bilinear form preserves the physical mass inner product and gives an identity mass matrix under a density-orthonormal auxiliary basis, leading naturally to an explicit central-difference scheme. We establish discrete energy stability under a coarse-scale CFL condition and derive convergence estimates for both the density-weighted displacement and the locally recovered stress. Numerical experiments in heterogeneous media validate the theoretical results and demonstrate the effectiveness of the proposed method for elastic wave propagation.
Given a real square matrix $A$ and a nonempty closed target set of complex matrices $\mathcal E$, invariant by complex conjugation and containing real matrices, does the subset of $\mathcal{E}$ consisting of the matrices nearest to $A$ in the Frobenius distance always contain a real matrix? We give negative answers when $\mathcal{E}$ is either the set of normal matrices (solving a question posed by N. Higham) or the set of matrices whose eigenvalues lie in the closed left half-plane (solving a question posed by the first author and F. Poloni). For the latter problem, we also argue that the ratio between the real and complex distances is unbounded in every dimension $n\geq 3$, and this holds for every distance induced by a unitarily invariant norm. For both problems and $n \geq 3$, we show that a real input $A$ uniformly drawn from the unit sphere $\|A\|_F=1$ has no real minimizer over $\mathcal{E}$ with probability strictly between $0$ and $1$. The paper is complemented by some further results that are valid for more general target sets $\mathcal{E}$.
This paper introduces an online evolve--filter--relax (EFR) strategy for long-term stability of Operator Inference (OpInf) reduced-order models (ROMs) of complex fluid flow simulations. The main novelty of the new EFR-OpInf strategy is the use of online (i.e., at the learned ROM online evaluation level) spatial filtering inspired from large eddy simulation to significantly improve long-term stability and predictive performance of standard OpInf. Furthermore, the EFR-OpInf strategy reduces, and in some cases even eliminates, the need for standard $L^2$ regularization, while providing physical interpretability for the OpInf hyperparameters. The EFR-OpInf framework is modular, readily integrated into existing OpInf workflows, and accommodates a user-selected ROM filtering strategy. We demonstrate the EFR-OpInf's effectiveness using a fully non-intrusive projection-based ROM filter, a ROM differential filter, and a hybrid projection-differential ROM filter. The new EFR-OpInf models are evaluated on a high-Péclet-number convection--diffusion--reaction problem that embeds the variation in one parameter and two unsteady Navier--Stokes problems that focus on predictions beyond a training horizon: a transitional two-dimensional flow past a cylinder and a three-dimensional turbulent minimal channel flow. Across these three scenarios, EFR-OpInf can reduce prediction errors by up to an order of magnitude relative to standard OpInf. Moreover, EFR-OpInf remains stable over long prediction horizons in cases where standard OpInf diverges. Depending on the ROM filter used, EFR-OpInf's computational cost is comparable to that of standard OpInf.
We study torsional rigidity as a function of the labeled vertices of a convex polygon. Starting from the distributed second shape derivative, we derive the Hessian with respect to vertex coordinates. At a regular polygon, dihedral symmetry makes this matrix block circulant in radial-tangential coordinates, reducing its spectrum to the eigenvalues of Hermitian matrices of order two. We also derive an exact second-variation Galerkin identity and guaranteed functional residual majorants. Finite elements approximate the PDE solutions entering the Hessian, and FLINT/Arb provides the interval arithmetic needed for certification. In the scale-invariant setting, we certify exactly four zero eigenvalues generated by similarities and $2n-4$ strictly negative eigenvalues for $5\leq n\leq25$. The regular polygons in this range are therefore strict local maximizers, modulo similarities, of torsional rigidity divided by area squared. The observed Hessian error decreases nearly quadratically with the mesh size; the transmission regularity needed to prove this rate is stated separately as a conjecture.
A structure-preserving symmetric Galerkin boundary element method (SGBEM) is developed for the anisotropic multi-phase Mullins-Sekerka problem in $\mathbb R^2$. The bulk harmonic problems are eliminated through a discrete Dirichlet-to-Neumann operator assembled from region-wise SGBEM discretization, leading to a boundary-only formulation coupled to the evolution of a polygonal curve network with triple junctions. The Galerkin approximation of the resulting boundary problem is shown to be uniquely solvable, and the fully discrete scheme is unconditionally stable with respect to the time step as long as the evolving polygonal network remains in an admissible class. The accuracy of the proposed method is demonstrated by a convergence test under coupled space--time refinement against a similarity-rescaled exact solution consisting of three concentric circles.
We study two-point fractional boundary value problems with uncertain input data, where randomness may enter through the coefficients and boundary conditions. To quantify the resulting uncertainty in the solution, we employ the generalized polynomial chaos (gPC) framework and develop a stochastic Galerkin formulation of the problem. A particular focus of this work is the convergence analysis of the resulting approximation. Rather than imposing assumptions directly on the stochastic coefficients appearing in the gPC representation, we introduce minimal regularity assumptions on the input data and use them to establish the properties required for the convergence analysis. Based on these results, we prove the convergence of the stochastic Galerkin approximation to the corresponding gPC solution. Numerical experiments are presented to illustrate the theoretical findings and to investigate the influence of random coefficients and boundary conditions on the statistical behavior of the solution.
We derive localized maximum-norm bounds for bending moments computed by the Hellan-Herrmann-Johnson method for the clamped Kirchhoff plate problem. The main difficulty is to localize the discrete stress without leaving the HHJ space or violating its kernel constraint. Using symmetric-curl potentials, we construct a kernel-preserving localization and connect local Green stresses with a global discrete Green stress. This argument separates the local interpolation error from a weaker global pollution term.Under explicit regularity assumptions on the auxiliary problems, the resulting estimates give optimal pointwise convergence for bending-moment values and elementwise first derivatives. No logarithmic loss occurs for positive polynomial degrees, whereas the lowest-order value estimate retains a logarithmic factor. Under additional reflection symmetry, symmetric recovery improves bending-moment values for even polynomial degrees and first derivatives for odd degrees. The proved gains are one third and one half of an order, respectively. Numerical experiments confirm this parity dependence and exhibit gains close to one full order, exceeding those established theoretically.
This paper studies the strong convergence order of structure-preserving stochastic theta Milstein methods for a class of index-$1$ stochastic differential algebraic equations (SDAEs) with time-dependent singular matrices and non-globally Lipschitz coefficients. The singular matrix is allowed to vary in time while preserving a fixed differential algebraic splitting, and the drift and diffusion coefficients may exhibit superlinear growth. By exploiting the index-$1$ algebraic-differential decomposition of the exact solution, we identify the Milstein coefficient of the reduced stochastic differential equation directly in the original SDAE variables and establish the well-posedness and constraint preserving property of the proposed method for $θ\in[1/2,1]$. Under a coupled monotonicity condition and suitable polynomial regularity assumptions, the method is proved to preserve the algebraic constraints at all time levels and to converge with strong order one in the root mean square norm. Numerical experiments confirm the structure-preserving property and the theoretical convergence order.
We develop and analyze a class of Hermite-type methods for linear advection-diffusioon equations on one-dimensional periodic domains. The methods evolve both nodal values and cell averages and are therefore referred to as hybrid-variable (HV) methods. Using Hermite interpolation theory, we construct the HV methods to arbitrary order of accuracy, and prove that the method is supraconvergent in the sense that the spatial order of accuracy of the method is larger than the local truncation error. In the second part of the paper, we prove using the theory of positive trigonometric series that all central HV methods for linear diffusion equations and linear advection-diffusion equations are stable at the semi-discretized level. Both the supraconvergence property and the stability of central HV schemes are verified by extensive numerical examples.
A tensor randomized average block coordinate descent method with heavy-ball momentum is proposed for the tensor least squares problem with respect to the t-product. Theoretical analysis for tensor block coordinate descent method is established and average techniques are applied to avoid the computation of tensor Moore--Penrose inverse. To further accelerate convergence, a heavy-ball momentum scheme is incorporated into the block coordinate descent method, where the step size and momentum parameters are adaptively determined by a two-dimensional minimal residual projection. Theoretical analyses give the convergence of the new method and provide an improved bound on the linear convergence rate. Numerical experiments further verify the efficiency of the proposed method in terms of the number of iterations and the better video recovery performance.