Papers

We present a high-order entropy-stable framework for the compressible Euler and Navier-Stokes equations on moving domains. The space-time mapping describing the domain motion is recast in an arbitrary Lagrangian Eulerian (ALE) formulation, in which the physical inviscid fluxes and the contributions induced by mesh motion are treated in a unified manner. At the continuous level, we prove that the proposed moving-wall boundary conditions are entropy conservative for the Euler equations and entropy stable for the Navier-Stokes equations. The no-slip condition is formulated in terms of the velocity relative to the moving wall, yielding a bounded inviscid contribution to the entropy balance, while the viscous terms contribute only entropy dissipation. Using diagonal norm summation-by-parts (SBP) operators together with appropriate numerical fluxes, these properties are extended to the semi-discrete formulation, resulting in nonlinear stability in the $L^2$ sense. The accuracy, robustness, and scalability of the proposed method are demonstrated in practice through an extensive set of numerical experiments, ranging from canonical two-dimensional verification cases to large-scale turbulent and supersonic simulations involving moving boundaries and fluid-structure interaction. The results confirm the suitability of high-order entropy-stable schemes for complex moving-domain problems across a broad range of flow regimes and multiphysics applications. Because the analysis relies on the SBP property rather than on a particular discretization, the framework naturally extends to a broad class of methods based on diagonal-norm SBP operators, including finite volume, finite element, and flux reconstruction schemes.

We present a natively multiscale operator learning method for the surrogate modeling of (numerical solvers for) multiscale partial differential equations (PDEs). The primary novelty of our method lies in a novel multiscale kernel frame function approximation technique. Leveraging this new kernel frame technique, we cast the operator learning problem as one of learning frame coefficients of output functions as a function of frame coefficients of input functions. The generalization step then automatically allows for a multiscale decomposition of the output functions. Our method is applicable to both tensor-product grids and point clouds. We present interpolation proofs, error estimates, and numerical convergence rates for our frame approximation. We the demonstrate the applicability of our method for the surrogate modeling of inherently multiscale PDEs. The new multiscale frame kernel method is significantly more accurate than popular neural operators on challenging problems from the literature, while simultaneously admitting an a posteriori multiscale decomposition upon generalization.

We study the exact discretization of the $L_2$-norm in the space spanned by the first $N$ Rademacher functions. It is shown that the sufficient number of nodes for discretization with the minimal number of nodes is equal to $N$ or $N+1$ and depends on the dimension $N$. The connection with Hadamard matrices and the Hadamard conjecture is demonstrated.

A novel dual-primal finite element tearing and interconnecting (FETI-DP) domain decomposition method (DDM) is introduced for solving Maxwell's equations. The proposed method is built upon a two-channel transmission condition (TC) enforcing both tangential-field and normal-flux continuity across subdomain interfaces. Two interface channels are separately constructed from the Faraday and Ampere-Maxwell equations. Simultaneously enforcing tangential-field and normal-flux continuity reduces interface jumps by several orders of magnitude relative to the widely used Robin TC. The resulting global interface equation exhibits markedly improved iterative convergence. Numerical experiments across different partitioning strategies, subdomain counts and volumes, mesh resolutions, and geometries demonstrate drastic reductions in iteration counts, while accuracy, massive parallelism, and scalability are preserved.

Multiway data such as image collections and video is ubiquitous, but the usual approach of flattening them into matrices discards the cross-mode structure that often carries the signal. The t-product and its generalization, the $\star_{\mathbf{M}}$-product, give a matrix-mimetic tensor algebra with a tensor SVD whose truncation is optimal in the Frobenius norm, just as in the matrix case. Much real data also has internal reflective symmetry: frontal faces, manufactured parts, and leaves are all bilaterally symmetric. We define a symmetry-preserving $\star_{\mathbf{M}}$-SVD that extends the matrix symmetry-preserving SVD of Shah and Sorensen to the $\star_{\mathbf{M}}$-algebra. When a tensor's transform-domain frontal slices are reflectively symmetric, the left basis of its $\star_{\mathbf{M}}$-SVD is symmetric too, so only its top half must be stored. Bilaterally symmetric images, turned on their side and stored as lateral slices, give such a tensor. This construction keeps the symmetric part of each image and stores only half the basis. We then recognize new images by projecting them onto this basis and matching to the nearest training image in coefficient space. Across ten datasets of faces, leaves, butterfly wings, and other objects, the symmetric basis matches the recognition rate of ordinary tensor SVD at 2-11x less basis storage, with the exception of MUCT. On faces under varying illumination it exceeds the best rate the ordinary tensor SVD attains at any storage level.

In this essay I reflect on the status of three of the Millennium Prize problems: the Riemann Hypothesis, P vs. NP, and the solvability of the Navier-Stokes equations.

Rapid data assimilation is required for real-time prediction and control in digital twins. For many physical systems, the data assimilation task requires the solution of a physics-constrained inverse problem, which is often computationally intractable in real time using traditional physics solvers. This work presents a reduced-basis neural operator approach to enable real-time inverse problem solutions in the digital twin setting. Our approach specifically targets the large class of problems with spatiotemporal dynamics governed by partial differential equations (PDEs) that are parameterized nonlinearly with respect to model parameters $m$, and linearly with respect to inversion parameters $q$. Based on this physical structure, our neural operator approximates the nonlinear map from the model parameters $m$ to the parameter-to-observable operator $\mathcal{F}(m)$ in a reduced subspace. Since the output of the neural operator is the parameter-to-observable operator itself (manifested as a matrix), we refer to this approach as NEural Matrix Operator (NEMO). With NEMO, for new given $m$, we enable a closed-form inverse problem solution for $q$ in a reduced subspace. We apply NEMO in two real-world applications: contaminant transport initial condition identification, and hypersonic vehicle load identification. We show that NEMO delivers high quality inverse problem solutions for data assimilation in real time, with over three orders of magnitude speedup compared to constructing the reduced operator with the PDE solver. Further, NEMO demonstrates comparable inverse performance to a state-of-the-art multiple-input neural operator, while reducing online computational complexity by over an order of magnitude and providing real-time uncertainty quantification.

For the Landau-Lifshitz equation, the discrete linear systems obtained by our semi-implicit method possess the following properties: they are large-sparse systems with non-Hermitian yet positive-definite coefficient matrices. To solve these systems efficiently, we apply the Hermitian/skew-Hermitian splitting (HSS) method and its inexact variant (IHSS). Numerical experiments in one and three dimensions show that the spectral radius of the HSS iteration remains below its theoretical upper bound and strictly below one for the tested grid resolutions and damping parameters. Moreover, the theoretical bound closely follows the actual spectral radius, providing an accurate estimate of the convergence behavior. The IHSS results demonstrate effective convergence for the tested cases and show that its efficiency is sensitive to the splitting parameter. Overall, the two semi-implicit schemes exhibit comparable convergence behavior.

Iterative diagonalization is the dominant cost of plane-wave density-functional theory (DFT), with search-space orthogonalization scaling particularly quickly with problem size and the number of target states. We present a randomized block Davidson-type eigensolver that replaces Euclidean orthogonalization with randomized Gram-Schmidt in a sketched inner product, requiring only a single pass over the basis while keeping its conditioning bounded independently of the input vectors. This modification changes only the Rayleigh-Ritz step, which becomes a definite generalized Hermitian eigenproblem. Ritz extraction remains exact, preserving true Ritz pairs and the interlacing property that makes each band energy an upper bound on the true one. The method is implemented in mixed precision for CPUs and GPUs from a single Julia code, interfaces matrix-free with DFTK, and is released in the open-source RandESC library. On sparse test problems with a fixed number of eigenpairs, the sketched solver overtakes its deterministic counterpart beyond matrix dimensions of about $2\times 10^4$ and is $25\%$ faster at $5\times 10^5$. In full self-consistent field DFT calculations, however, both Davidson variants outperform the locally optimal block preconditioned conjugate gradient (LOBPCG) reference only by $5$ to $11\%$ in total time, while the additional benefit of sketching is limited. As the number of requested states grows with system size, orthogonalization savings are offset by the generalized eigenproblem. Therefore, the regime in which sketching pays off is set by how the number of wanted states scales with the problem dimension, not by the eigensolver as such.

Profitability metrics, such as the single-product gross profit margin, constitute the standard toolbox to prioritize products in sales and operations planning (S&OP). These metrics leverage the premise that production lines can be committed independently, which endows them with a simple per-ton ranking that a planner can compute and use easily in practice. However, coupled plants, especially polymer plants whose parallel lines draw simultaneously on a single bulk feed, can fundamentally undermine this premise and mislead to value-destroying production plans. We propose a framework built on what we call the Augmented Gross Profit for Product Clusters (AGPPC), which converts the per-ton ranking metric into a per-coupled-hour ranking that carries a per-instance optimality certificate. We present the theory of AGPPC, which combines the co-production column concept with a fluid relaxation of the planning problem, alongside an integrated bilinear mixed integer program of the plant and a McCormick-linearized baseline that provides certified bounds for it, and we demonstrate its effectiveness on real operational data from an Indonesian polymer producer, where the coupling-aware metric alone raises operating profit by 7.6% and the full optimization by 28% over current practice.

In this work, we develop an accurate numerical homogenization method for computing effective Hamiltonians of quasiperiodic Hamilton–Jacobi equations (QHJEs) with convex Hamiltonians of the form $H(x,p) = |p|^k/{k}-f(x), ~k>1$, where $f$ is quasiperiodic. Computing effective Hamiltonians in the quasiperiodic setting requires solving QHJEs posed on the whole space. Their solutions generally possess neither translational symmetry nor decay and may exhibit low regularity. These features pose substantial challenges for numerical computation. To address these difficulties, we introduce a quasiperiodic boundary condition, which allows the original whole-space problem to be treated on a bounded domain while preserving quasiperiodicity at the boundary. We then propose an SL–FPR scheme that combines a semi-Lagrangian (SL) approximation with the finite points recovery (FPR) method and establish stability and error estimates for the resulting scheme. We also extend the quasiperiodic homogenization result from the quadratic case to general $k>1$ and apply the proposed method to accurately approximate the corresponding effective Hamiltonians. Numerical experiments illustrate the convergence and applicability of the method and validate the extended homogenization results.

We study three-dimensional wave scattering by high-index dielectric resonators with Kerr-type nonlinearity under plane-wave incidence, together with the associated nonlinear dielectric scattering resonances, through a nonlinear Lippmann-Schwinger equation. Using a Lyapunov-Schmidt reduction near a simple eigenmode of the Newtonian potential, we obtain a local decomposition of the scattered wave into a resonant contribution and a controlled remainder and derive an explicit nonlinear equation for the resonant modal coefficient for sufficiently small contrast parameter and locally small resonant amplitude and incident field. A second reduction covers the regime of incident fields of order one and weak nonlinearity. Our results extend for the first time the linear modal decomposition to nonlinear wave-scattering problems. We complement them by developing a numerical framework based on Nyström discretization, real Newton iteration, and pseudo-arclength continuation. For single resonators of several geometries, our computations confirm the predicted high-contrast resonance scaling and the nonlinear modal approximation, and exhibit the multivalued incident-wave response. For a mirror-symmetric resonator dimer, we track symmetric, antisymmetric, and symmetry-broken resonance branches over a broad range of separations. Our computations show that the symmetry-breaking threshold increases as the gap between the resonators decreases, in agreement with the leading-order bifurcation theory.

We present a samplet-based framework for the efficient numerical solution of saddle-point systems arising from conditionally positive definite (CPD) kernel approximation in general and universal Kriging in particular. The vanishing moment property of samplets as well as the particular structure of the associated scaling distributions, which correspond to discrete orthogonal polynomials, allow for a numerically favorable representation of these saddle-point systems. Concretely, they enable a natural null-space reduction of the indefinite saddle-point system to a (large) linear system for the detail coefficients and a small triangular system for the polynomial coefficients. We derive error bounds for the approximation by polyharmonic splines in Beppo-Levi spaces and show that the detail coefficients span precisely the subspace on which the CPD kernel is positive definite, rendering the reduced block symmetric positive definite. In view of the quasi-sparsity of the samplet-transformed kernel matrix for asymptotically smooth kernels, the resulting method achieves O(N log N) cost for the assembly and the storage of the saddle-point system. The reduced system can efficiently be solved by a sparse Cholesky factorization. We illustrate the framework with three applications, namely Gaussian process regression with generalized covariances, landmark-based image registration via samplet-compressed thin plate splines, and three-dimensional mesh deformation.

With the rise of scientific machine learning, physics-informed neural networks (PINNs) have been extensively applied to a wide range of problems. Nevertheless, most theoretical analyses of PINNs remain confined to linear equations, and a substantial gap persists between PINNs and classical numerical analysis for nonlinear problems. To address this issue, we propose a robust training framework for PINNs solving nonlinear partial differential equations, together with a rigorous error analysis. Specifically, for the $p$-Laplace equation, we introduce a novel loss formulation that combines a dual residual loss measured in $W^{-1,p'}$ with a boundary loss measured in a fractional Sobolev norm $W^{1-\frac{1}{p},p}$. This formulation is designed to accommodate the limited regularity of weak solutions and enables us to establish rigorous a priori and a posteriori error estimates. Moreover, the proposed framework and its analysis are extended to a parametric setting in which the exponent, the source term, and the boundary condition may all vary with the parameters. Finally, we present numerical experiments that substantiate our theoretical findings.

We develop eigenvalue bounds for symmetric, block-tridiagonal multiple saddle-point linear systems, preconditioned with block-triangular matrices, based on approximate Schur complements. Irrespective on the number of blocks, we prove that all complex eigenvalues, with nontrivial imaginary part, are strictly contained in a circle within the complex plane, with center 1. The real and positive eigenvalues are bounded in terms of the extremal roots of a sequences of parametric polynomials. Numerical results reveal that the bounds describe very well the eigenvalue distribution of the preconditioned matrix.

This paper develops a fully discrete finite element method for three-dimensional incompressible magnetohydrodynamic (MHD) flows with nonlinear $p$-Laplace viscosity. The scheme combines spatial finite elements with a semi-implicit Euler time discretisation. Convergence to a weak solution is proved using a time-translation compactness argument together with Minty's monotonicity method. Under additional regularity assumptions, we derive unconditional error estimates for both the velocity and magnetic field, with no coupling restriction between the time step and mesh size. The framework extends finite element analysis of incompressible MHD systems to non-Newtonian shear-thickening fluids ($p>2$), while recovering the classical Newtonian case ($p=2$). Finally, numerical simulations are provided to validate the theoretical convergence rates and demonstrate the robustness of the proposed method.

In previous work by several authors, the behavior of the condition numbers of banded Toeplitz matrices was studied as the matrix size tends to infinity. In the present contribution, two main directions are pursued. As a first step, we extend this study to block Toeplitz matrices with blocks of fixed size $N$. As in the scalar case, we show that even when the symbol generates a Fredholm infinite Toeplitz operator, the condition numbers of the finite matrices may grow at least exponentially. Upper and lower bounds for the condition numbers are obtained, and examples showing that they may grow arbitrarily fast are presented. Then, as a second step, we apply the developed theory to the stability analysis of space-time Galerkin methods, where in time an Isogeometric approach is used with regularity $r$, $1\le r\le p-1$, $p$ being the employed polynomial degree. These stability issues are related exactly to the conditioning of block Toeplitz-like matrices with blocks of size $N=p-r$. Specific examples are treated in detail and related numerical experiments are presented and critically discussed. We finally present a short list of relevant open problems.

Phase-field modeling of brittle fracture removes the need to track cracks explicitly by recasting their evolution as the minimization of an energy functional. In return it requires a discretization dense enough to resolve a localization band whose width is set by a regularization length and whose path is not known in advance. We propose a mesh-free discretization in which a single neural network represents the displacement and phase fields and is trained by minimizing the incremental energy directly. The coordinates enter the network through a multiresolution feature encoding built from $C^1$ quadratic B-spline grids, so the finest scale the representation can express is set by choice rather than reached through slow training, and the energy is estimated by stratified Monte Carlo integration on points redrawn at every optimizer iteration. This pairing proves critical, since the crack fails to advance both when the integration points are held fixed and when the encoding is too coarse to represent the band, while each ingredient tolerates a wide range of settings once the other is in place. Because the representation is globally $C^1$, the second- and the fourth-order fracture energy densities run on the identical discretization. Across six problems, from single-edge-notched tension and shear to a thick-walled ring on a single spline patch, the computed load-displacement curves follow staggered finite element references at matched regularization length, with peak loads within about 1% on the single-edge-notched tests and within 8% where the crack pattern changes topology. On a public benchmark dataset of random multi-crack configurations the method classifies the active or dormant state of 90% of the seeded cracks in twenty zero-shot runs, where the deep Ritz baseline of the dataset authors fails.

Splines, which are piecewise polynomial functions with given smoothness, are used in numerous applications such as computer aided geometric design, algebraic geometry and for implementation of the finite element method. In many applications, it is important to know the dimension of the space of splines over a domain, but this continues to be an open problem. Strang's conjecture predicts the dimension of the space of bivariate $C^1$ splines of degree at most $k$ over a polygonal domain with no holes. For the cases, $k=4$ and $k \geq 5$ the conjecture has been proven true. For the case $k=2$, John Morgan and L. Ridgway Scott provided a counterexample. In this project, we begin by relating the space of bivariate $C^1$ splines of degree at most $k$ to various other function spaces. With the aid of some intrinsic operators on these spaces, we obtain a complex that yields an equivalence statement to Strang's conjecture. The equivalent statement admits a computational approach to the conjecture and thus we present an NGSolve code to check Strang's conjecture for a given mesh. In addition, we prove that if the simplicial complex associated to a triangulation is strongly collapsible (with a non-collinear condition), then Strang's conjecture holds for $k=2$.

We investigate an inverse coefficient problem for a multidimensional free-boundary model of ductal carcinoma in situ (DCIS), in which the tumor interface is governed by the nonlinear coupling of nutrient concentration, tissue pressure and curvature, and the unknown nutrient consumption function is recovered from a temporal trace of the nutrient concentration obtained by needle aspiration biopsy. For the forward problem, we establish uniform local well-posedness over an admissible class of consumption functions. The inverse problem is recast as a fixed-point problem: approximating the admissible set by finite-dimensional spaces yields discrete iteration operators, for which we prove the existence of fixed points, and the strong convergence of a subsequence of discrete fixed points to a fixed point of the continuous operator, which solves the inverse problem under a consistency condition. To approximate these fixed points, we develop a homotopy-continuation method combining a linearly convergent Picard iteration with a cubical Sperner search, without differentiating an objective functional or computing an adjoint state. Several numerical experiments on radially symmetric and non-symmetric DCIS models corroborate the theoretical findings.

Immersed finite element (IFE) methods provide an effective framework for solving interface problems on unfitted meshes. The basic idea underlying IFE methods is to modify standard finite element spaces by enforcing interface conditions at certain points. It is known that, for nodal degrees of freedom, unisolvence of linear IFE basis functions on triangular elements is subject to a non-obtuse-angle condition for scalar diffusion coefficients. This limitation is more severe for tensor diffusion coefficients, for which the non-obtuse-angle condition is no longer sufficient for unisolvence. Even on rectangular meshes, conventional bilinear IFE basis functions may not be unisolvent for tensor diffusion coefficients. In this paper, we develop and analyze a nodal isoparametric IFE method on quadrilateral and hexahedral meshes and show that the unisolvence issue can be overcome by appropriately selecting the enforcement point of the discrete flux condition. The key observation is that, unlike triangular elements, the discrete flux in quadrilateral and hexahedral elements is not constant, which provides the flexibility to select such an enforcement point to ensure unisolvence. We provide a systematic procedure for this selection that not only ensures unisolvence of the IFE basis functions on general quadrilateral and hexahedral elements with either scalar or tensor-valued diffusion coefficients but also preserves the optimal approximation properties of the resulting IFE space. The proposed IFE method offers several advantages: flexibility for complex geometries, the absence of angle restrictions, and applicability to tensor diffusion coefficients, thereby overcoming the limitations of existing nodal IFE methods on rectangular and triangular meshes. Optimal error estimates are established and confirmed by numerical experiments.

We propose a hybrid sharp-diffuse interface representation for modeling intermittent and isolated topological transitions during Cahn-Hilliard phase coarsening. Away from topological events, the evolution is approximated by the Mullins-Sekerka sharp-interface limit system and computed using a boundary integral formulation. When diffuse transition layers overlap, a novel interface surgery algorithm resolves topology changes through a localized Cahn-Hilliard pseudo-time evolution, after which the sharp-interface calculation is resumed. We develop the mathematical formulation underlying this decomposition, present a simple two-dimensional numerical implementation, and simulate a mass-exchange problem with interface coalescence. By localizing the diffuse evolution to topological events, the method achieves a speedup of two to three orders of magnitude over conventional diffuse-interface simulations.

This paper is concerned with a structure-preserving neural network-based framework for the Cahn-Hilliard equation in mixed form. We employ generalized transferable neural networks (GTransNet) for spatial approximation and stabilized backward differentiation formulas (BDF) for temporal discretization. The resulting first- and second-order in time GTransNet-BDF schemes are shown to conserve mass and satisfy energy stability at the time-discrete level. The schemes are implemented by a collocation-based method, in which a least-squares system with constant coefficient matrix needs to be solved at each time step. The solution of this system, which determines the output-layer weights of the network, violates mass conservation due to the expected nonzero least-squares residual. To overcome this issue, we introduce a novel post-processing mass-conserving projection that enforces the mass constraint through a minimization problem, whose solution can be computed at negligible computational cost. A key advantage of the proposed method lies in its predetermined hidden layers and mesh-free nature, making the method applicable to complex domains, variable mobility, and long-time simulations. Extensive numerical experiments in two and three dimensions verify convergence, mass conservation, and energy dissipation as well as demonstrate the accuracy and robustness of the proposed GTransNet-BDF schemes.

In this work, we investigate implicit-explicit (IMEX) arbitrary high-order Deferred Correction (DeC) methods for the approximation of ordinary differential equations (ODEs). Such schemes are characterized by an iterative procedure that increases the order of accuracy by one at each iteration. More precisely, we study an efficient modification based on the introduction of interpolation processes between consecutive iterations, with the aim of systematically matching the accuracy achieved at each iteration with the order of the discretization employed. On the one hand, this modification leads to computational advantages, since the low-order iterations are performed on cheaper lower-order discretization structures; on the other hand, it endows the methods with a natural $p$-adaptive character, which is particularly appealing in the context of practical applications. We investigate this modification for two families of DeC schemes, providing numerical validation, efficiency assessments, and stability region plots. The numerical validation includes several examples involving stiff ODEs and partial differential equations (PDEs) with high-order spatial derivatives. The ability of the modified schemes to provide high-fidelity results at reduced computational cost, as well as the effectiveness of the adaptive strategy, is demonstrated through the numerical experiments.

In this report, we recall Christiansen's [2011, Num. Math.] discrete elasticity complex that uses spaces of polynomials and distributions. Using the degrees of freedom of those spaces, we write a complex using vector-valued chains and co-chains that is isomorphic to Christiansen's complex. The operators of the new complex are algebraic, but use the geometry of the simplicial complex. Finally, we define a de Rham-like map from the smooth elasticity complex to this discrete complex.

The ultra-weak variational formulation (UWVF) is known to be the Trefftz discontinuous Galerkin scheme obtained from the standard PWDG flux family by setting \(\alpha=\beta=\delta=1/2\). We give a short structural explanation of this fact. On every interior mesh edge, the two PWDG jump terms at these parameter values are exactly one half of the sum of the squared norms of the two oriented impedance-transmission defects used by the UWVF. Thus the Trefftz-DG norm is the natural Riesz norm of the global transmission defect. In PWDG coordinates its Gramian is the anti-Hermitian part of the discrete matrix, and Riesz normalization puts the operator in the normal form \(K+\ii I\), with \(K\) Hermitian. The same argument applied directly to the discrete UWVF yields the complementary normal form \(\tfrac12 I+\ii K_{\rm U}\). This identifies the classical equivalence as a relation between broken-field and impedance-trace coordinates for the same edge-coupling problem, rather than only as a particular choice of penalties.

We develop a Monte Carlo lattice Boltzmann method to compute statistical solutions to the three-dimensional incompressible Navier–Stokes and Euler equations. Entropic space-time adaptive relaxation of the higher order kinetic moments yields stable numerical solutions with decreasing viscosity. We provide a convergence analysis that is conditional on four explicitly stated assumptions regarding the discrete dynamics. Under diffusive scaling, the laws of the discrete ensemble converge along a subsequence to a limit satisfying the Foias–Temam Liouville formulation of the Navier–Stokes equations. Consequently, provided the structure function scaling holds uniformly, the vanishing viscosity limit of these measures satisfies the multi-point statistical Euler hierarchy of Fjordholm, Mishra, and Weber. The limit measures inherit the known weak-strong uniqueness principle on the interval of existence of a strong Euler solution. Under explicit scaling assumptions, a Kuznetsov-type argument yields a fractional 1-Wasserstein convergence rate. We present three-dimensional computations of time-dependent statistical solutions along the inviscid limit of the incompressible Navier–Stokes equations together with convergence measurements in the Wasserstein metric. Numerical experiments on a randomized Taylor–Green vortex with 24-dimensional initial uncertainty recover Kolmogorov's K41 scaling for energy spectra and structure functions, exhibit the failure of pathwise strong convergence, and yield Wasserstein convergence rates of about $0.5$ at the onset of turbulence. Finally, error measurements with respect to spectral hyperviscosity computations indicate that the computed limit measure is independent of the numerical regularization.

Iteratively reweighted least squares (IRLS) methods constitute a natural approach to nuclear norm minimization, but their convergence rates and the role of the weight operator have remained poorly understood. This paper establishes sharp convergence rates for IRLS methods for constrained nuclear norm minimization in low-rank recovery. A central ingredient is a new majorization analysis for the smoothed nuclear norm: we prove that the harmonic-mean weight operator defines a valid global quadratic majorizer. Furthermore, we show that this weight operator is optimal within the family of power-mean weights, clarifying why it improves over classical one-sided reweighting schemes that use only row- or column-space information. Under a Schatten-1 null space property, we prove global linear convergence of IRLS algorithms using a variety of weight operators, including the harmonic-mean weights. For IRLS with harmonic-mean weights, we prove a dimension-independent, locally linear convergence rate. We provide a counterexample showing that this dimension-independent local rate cannot in general be obtained for IRLS algorithms using one-sided weight operators, which predominate in the literature. Numerical experiments corroborate the theoretical results and illustrate the practical advantage of harmonic-mean reweighting across square, rectangular, and adversarially initialized recovery problems.

The human diaphragm is vital for respiration and its simulation can aid in further understanding complications arising in individuals that have been subjected to mechanical ventilation. The diaphragm geometry is very thin and hence its numerical simulation poses a challenge for most numerical methods. In this paper, we compute the deformation of real diaphragm geometries which are modelled as linearly elastic and are subject to realistic boundary conditions. The method used is the unfitted least-squares radial basis function partition of unity method, which is well suited to solve problems on thin geometries given its unfitted nature and the shape and refinement of the cylindrical patches that are used. We derive a theoretical proof showing that the method converges with high order to the solution of this problem.

In this paper, we introduce the Inertial Manifold Neural Operator (IMNO) for solving dissipative time-dependent partial differential equations (PDEs). The long-time dynamics of such systems often exhibit an effective low-dimensional structure due to dissipation. Unlike standard neural operator architectures such as the Fourier Neural Operator (FNO), IMNO explicitly leverages the low-dimensional structure to achieve better physical interpretability, accuracy, and stability in long-horizon autoregressive training and prediction for nonlinear dissipative PDEs. For shift-equivariant PDEs, we further introduce a shift-equivariant variant (IMNO-SE) of the proposed neural operator, ensuring that a spatial shift in the input induces the same spatial shift in the output. This symmetry-preserving inductive bias substantially improves its performance in shift-equivariant PDEs. Extensive benchmark experiments are presented to evaluate IMNO's performance numerically.