Papers

Alefeld (1981) recast Halley's cubic convergence as Newton's method on $g=f/\sqrt{f'}$, linked to the Schwarzian derivative. We generalize this to matrix gradient fields $F=\nabla\phi$ ($X\in\R^{m\times n}$) via a matrix Schwarzian derivative interpreted through information geometry ($\alpha$-connections). Simplifying Palmore (1994), we construct this operator square-root-free derivative from the third Fréchet derivative of the Newton map. Our main theorem proves local cubic convergence with an explicit error constant, without self-adjointness, commutativity, or gradient-field assumptions. A matrix-free algorithm (Hessian-vector products only) validates the theory. We contrast this with a power-Newton family (whose naive matrix extension fails) and the matrix Laguerre family (which requires an operator square root). A case study on the Oja-type system $\dot X=AXB-XBX^TAX$ reveals that convergence depends on target eigenvalue gaps, cubic gains grow with ill-conditioning, and 70-digit tests confirm exact theoretical orders alongside a working-precision accuracy budget. Finally, we note coupled formulations excel primarily when sequential deflation is inapplicable.

Although machine learning can be used to predict the evolution of physical systems from data, a formulation that learns only the system state at each time leaves the temporal and dynamical structure of the solution to be resolved within a broad hypothesis space. We introduce Learning$^2$, a representation-level framework that structures this space by coupling a primary representation to a second representation through a known physical transformation. The resulting cross-representation constraint restricts the effective hypothesis space and provides an ante-hoc, physically interpretable criterion for excluding solutions that satisfy the primary representation alone. We instantiate Learning$^2$ through EuLaNet, an Eulerian–Lagrangian representation for fluid dynamics. Given the velocity state $u(\mathbf{x},t)$, EuLaNet constructs its induced Lagrangian flow map $X(\mathbf{a},t)$ through $\dot{X}(\mathbf{a},t)=u(X(\mathbf{a},t),t)$, from which material transport and finite-time deformation are derived. The resulting representation couples the predicted state to the dynamical consequences it induces, providing a second consistency criterion beyond state-level agreement. We formalize this construction through an effective hypothesis space $\mathcal{H}_{L^2}\subseteq\mathcal{H}$ and define the conditions under which a consequence representation provides discriminative constraints on candidate solutions. EuLaNet is implemented as a model-independent representation module, separating the physical constraint from the downstream learning architecture. This construction provides an ante-hoc mechanism for physically interpretable constraint in scientific learning and offers a basis for developing and evaluating broader classes of Learning$^2$ architectures. The implementation is open-sourced to support the development and extension of the architecture across scientific domains.

Recent developments have shown that eigenvalue localization regions can often be significantly sharpened by identifying and excluding subsets that are guaranteed to contain no eigenvalues, leading to the concept of inclusion-exclusion localization sets. In this paper, we introduce a new family of inclusion-exclusion eigenvalue localization sets based on the S-strict diagonal dominance framework and the corresponding CKV localization regions. For an arbitrary nonempty subset S of the index set, we construct novel exclusion regions and derive the associated S-SDD inclusion-exclusion localization sets. We prove that every eigenvalue of a complex matrix belongs to the proposed regions and define a refined localization set obtained by intersecting the corresponding inclusion-exclusion regions over all nonempty subsets S. The proposed framework unifies and extends several existing results. In particular, we show that the classical Gershgorin inclusion-exclusion localization set and the Dashnic-Zusmanovich inclusion-exclusion localization set arise naturally as special cases of our construction. Consequently, the new localization set provides, in general, a sharper eigenvalue inclusion region than both of these previously known localization regions. The theory is further extended to matrix polynomials, yielding corresponding inclusion-exclusion localization results for polynomial eigenvalue problems. Several numerical examples illustrate the effectiveness of the proposed approach and demonstrate the improvement achieved over existing localization sets.

The stochastic Landau–Lifshitz–Gilbert (sLLG) equation is a strongly nonlinear stochastic PDE with a non-convex pointwise constraint arising in the theory of micromagnetics. We analyse a fully discrete, structure-preserving finite element approximation of the sLLG equation with coloured multiplicative Stratonovich noise on a bounded interval. The method utilises continuous piecewise affine finite elements, mass lumping, and midpoint time discretisation to preserve the unit-length constraint exactly at the finite element nodes. Under suitable regularity assumptions on the initial data and the noise, we establish uniform higher-moment stability and develop an error analysis for the scheme. The analysis exploits the geometric structure of the equation and the stochastic midpoint discretisation. For every $\gamma\in(0,\frac12)$, we prove first-order spatial convergence and temporal convergence of order $\gamma$ in the natural discrete energy norm, locally in mean square on events of arbitrarily large probability and, consequently, in probability. To the best of our knowledge, this is the first convergence-rate result for a fully discrete structure-preserving finite element scheme solving the stochastic Landau–Lifshitz–Gilbert equation.

In this paper, we present a high-order, stable summation-by-parts (SBP) finite difference method for solving the Vlasov-Maxwell system in a 2D2V phase space. Central SBP operators for the advection terms are not stable when the solution becomes non-smooth and fine-scale filamentary structures develop, as is typical in high-dimensional Vlasov-Maxwell simulations. To address this issue, the method is stabilized using high-order upwind SBP operators. High-order explicit Runge-Kutta methods are employed for time integration. We prove that the fully discrete scheme exactly conserves mass and preserves momentum up to truncation error. Furthermore, we present a matrix-free implementation of the method on modern GPU architectures. A range of challenging benchmark problems is solved to demonstrate the accuracy, robustness, and performance of the proposed scheme.

Let $\mathbb D$ be the normalized ridge dictionary generated by $\operatorname{ReLU}^k$ on a bounded Lipschitz domain $\Omega\subset\mathbb R^d$. We determine the sharp algebraic rate of finite $n$-term approximation from $\mathbb D$ when the outer $\ell^1$ coefficient budget is independent of $n$. More precisely, let $d\ge3$, $k\in\mathbb N_+$, $0\le m\le k$, $0<p<1$, and $0<q\le1$. For the unit ball of the critical Besov space $B_{p,q}^{k+d/p}(\Omega)$, measured in $H^m(\Omega)$, the optimal algebraic exponent is \[ \min\left\{\frac{k-m+d/2}{d-1}, \frac{k-m+d/2+1/p-1/2}{d}\right\}. \] We prove two-sided estimates in the three regimes determined by the two branches of this minimum and give the corresponding logarithmic factors at and beyond the transition. The estimates coincide without logarithmic loss in the strict angular regime and at the transition when $q\le[1/2+(k-m+d/2)/(d-1)]^{-1}$. The upper estimate combines critical wavelet sparsity, localized Fourier–Radon representations, and a stable allocation of directions and biases. The lower estimates arise from two distinct obstructions, namely radial ridge approximation on the direction sphere and Gevrey localization in the joint direction–bias space. The critical smoothness is exactly the endpoint at which Besov regularity provides a uniform shallow-network variation bound without additional scale decay. The resulting representation exponent is strictly larger than the degree-limited exponent for fixed-degree isotropic finite elements under a parameter-count comparison. This last statement concerns best approximation, not training or computational complexity.

We construct a nondegenerate conforming straight-line triangulation of a polygonal disk on which some vertex data admit no continuously differentiable piecewise polynomial interpolant of total degree at most three. The triangulation has 41 vertices and 56 triangles, and the graph induced by its interior vertices is a tree with eight arms of length two. Every cubic $C^1$ spline on this triangulation satisfies an explicit linear relation with integer coefficients among its vertex values. We derive the relation by a weighted sum of Bernstein–Bézier smoothness equations and verify it using rational coordinates and weights. This disproves Alfeld's conjecture on vertex interpolation. The nonexistence proof does not require a matrix rank computation. A separate computation in exact arithmetic gives rank 40 for the vertex evaluation map and dimension 107 for the spline space, attaining the classical dimension lower bound for this triangulation.

We consider a Type-II variational principle driven by geometric rough path with split boundary conditions naturally suited to adjoint systems. From this rough variational principle we derive rough Hamilton's equations, establish their pathwise conservation laws and associated Hamilton–Jacobi equation. We then specialise the framework to rough adjoint systems, obtaining pathwise conservation and quasi-conservation laws that underpin adjoint sensitivity analysis with respect to initial conditions and parameters. On the discrete side, we construct a rough Galerkin discretisation of the rough Type-II variational principle and show that it generates a symplectic flow with discrete analogues of the continuous conservation laws. We establish it's equivalence to a class of Rough Symplectic Partitioned Runge–Kutta (RSPRK) methods and analyse its convergence and naturality properties. Lastly, we perform numerical experiments to validate the predicted convergence rates and demonstrate that RSPRK methods preserve the adjoint conservation laws to machine precision, yielding more accurate and stable gradients in optimisation problems than non-symplectic alternatives.

We analyze an inverse problem for determining the piston speed and the associated flow field from a prescribed leading shock and the initial data in a shock tube. The gas flow is described by the isentropic Euler equations (i.e., the $p$-system), while the trajectory of the leading shock is prescribed as a given Lipschitz curve. Under an Oleĭnik-type entropy condition on the leading shock, we develop a modified wavefront tracking scheme to construct the flow field behind the shock. This construction enables us to determine the corresponding piston speed and the associated flow field.

Strongly enforced Dirichlet boundary conditions require highly refined near-wall meshes to resolve steep velocity and thermal gradients. This introduces high computational costs, especially for practical flow simulations. Weakly enforced boundary conditions alleviate this burden by acting as a variationally consistent near-wall model. By allowing a controlled slip at the solid wall, weak enforcement recovers accurate flow quantities on coarse boundary-layer meshes across both incompressible and compressible regimes. Furthermore, weak boundary conditions serve as the fundamental enabling technology for immersogeometric analysis. Because the weak operator evaluates boundary integrals independently of the background mesh, high-fidelity flow analysis can be performed directly on complex geometries without fitting a mesh to the surface. This capability has facilitated direct geometry-to-analysis workflows for boundary-representation CAD models, raw point clouds, photogrammetric reconstructions, and segmented medical images. This review examines the unified mathematical development of the weak boundary condition framework from scalar advection-diffusion equations to the full Navier-Stokes equations, illustrating its versatility and robustness across incompressible and compressible flows, whether in traditional boundary-fitted, sliding-interface, or advanced immersogeometric applications.

Randomly pivoted Cholesky, QR, and LU are iterative algorithms that form structured low-rank approximations of a matrix by sampling columns, or rows and columns, from the residual. We give convergence rates for these methods that depend on the decay of the singular values of the original matrix. The rates hold, up to a constant, when pivots are drawn from an approximation to the exact distribution. This extra degree of freedom can be used to trade exact sampling for cheaper approximations. We use it to establish rates for variants of randomly pivoted LU, including one that picks pivots using a sketch of the input matrix.

We numerically implement full quantum circuits of optimal tapered Quantum Phase Estimation (in the form of an implementable approximate bandwidth-limited DPSS taper) for the purpose of minimum eigenvalue computation, where the Hamiltonian time-evolution unitary is implemented using standard Trotterization, and the taper state is constructed using an optimized bandwidth-limited taper. We illustrate, on representative quantum Hamiltonians, with tractable exact classical numerical quantum circuit simulations, how tapered QPE performs with respect to optimal phase sampling rate, absolute error, Shannon entropy of the optimal phase distribution, when the tQPE algorithm parameters are changed. Those algorithm parameters are Trotter error, finite sampling, imperfect initial states, total evolution time, and digitization error (number of phase qubits). The reported numerical experiments include up to $\ell=10$ phase qubits of precision with $m=3$ additional phase-register error suppression qubits, on a $4$-qubit Heisenberg quantum magnet model Hamiltonian, using up to 12th-order Trotterization. We numerically show i) the steady-state optimal phase sampling rate is determined by the initial state overlap with the ground-state, ii) with respect to phase sampling probability, or overall eigenvalue sampling error rate, there is no strong evolution time dependence for tQPE unless there is high Trotter error, iii) the approximated time-evolution unitaries in tQPE, like in standard QPE, can result in substantial ringing, which leads to non-physical eigenvalue estimates.

We introduce the weighted Laplace space $H_w$, an RKHS of Laplace transforms on $(0,\infty)$, and study spectral measures in its dual space $H'_w$. For conforming FEM discretizations of the Dirichlet Laplacian on bounded Lipschitz domains, we prove the dual-norm inequality $\|\mu_h\|_{H'_w} \leq \|\mu\|_{H'_w}$, where $\mu = \sum_k\delta_{\lambda_k}$ and $\mu_h = \sum_k \delta_{\lambda_{k,h}}$. The proof combines min-max monotonicity of FEM eigenvalues with a heat-trace representation of the dual norm. We then analyze $H_w$-adapted rational approximation of shifted symbols $\phi(x)=(x+\kappa^2)^{-\beta}$ and give a conditional transfer principle for estimates proved in the corresponding weighted Laplace pre-image norm. Via dual pairing, the norm inequality yields uniform bounds for finite spectral sums and related transformed observables.

The Fourier-cosine (COS) method of Fang and Oosterlee (2008) recovers densities and computes expectations semi-analytically from characteristic functions (ch.f.s). In higher dimensions, both expansion synthesis and COS coefficient-tensor construction scale exponentially with dimension. We address this curse of dimensionality with two tensor-train (TT) methods. A straightforward TT decomposition of the COS coefficient tensor gives COS-TT, which removes the exponential cost of online synthesis but not of offline decomposition. Our main contribution, COS-TT-CHF, instead compresses a ch.f. sample tensor and maps it to the original Fourier–cosine coefficient tensor through an alternative COS representation. Apart from the model-dependent cost of one ch.f. evaluation, both offline and online costs scale linearly with dimension, while error grows only algebraically. COS-TT-CHF introduces two additional errors: frequency-domain truncation and discretization errors. Both are controlled by parameter-selection rules from theoretical error analysis. We further propose an approach for expectations of nonseparable functions of additive scalar aggregates by computing the ch.f.s of the aggregates using tensorized COS methods, reducing the original multidimensional problem to a one-dimensional COS calculation. Closed-form formulas are derived for the TT-contraction integrals arising when the aggregate is a weighted sum, as in European basket option pricing. Experiments under geometric Brownian motion (GBM) and variance gamma (VG) show that, with the parameter-selection rules, COS-TT-CHF remains accurate through 150 dimensions for GBM and 100 for VG within a 30-minute offline computation budget on a laptop. As a by-product, our Fourier-cosine TT constructions yield cosine-basis functional TTs (FTTs) under weaker assumptions than in spectral FTT literature.

There is a close relationship between divided differences and hypergeometric series, and recent studies have shown that divided differences can be used effectively to derive terminating hypergeometric identities. In this paper, we apply this approach to terminating zero-balanced hypergeometric series. Using explicit product evaluations arising from divided differences, we give alternative proofs of zero-balanced ${}_3F_2(1)$ and ${}_4F_3(1)$ summation formulas and then extend the argument to the general terminating ${}_{r+1}F_r(1)$ case. We also combine the Lagrange representation with the Leibniz rule for divided differences to derive a finite convolution transformation for a terminating ${}_{r+2}F_{r+1}(1)$ series. Thus, the divided-difference approach not only reproduces the known zero-balanced summation formulas but also yields a finite convolution transformation for terminating hypergeometric series, together with its zero-balanced specialization.

The Levin method transforms the evaluation of a highly oscillatory integral into the solution of a first-order linear ODE for a slowly varying auxiliary function. This ODE is typically approximated by collocation, after which the integral value is recovered from the auxiliary function at the endpoints. The present work develops a new extension of the Levin method for the summation of one-dimensional and multidimensional infinite oscillatory series. The summation problem is transformed into the solution of a functional equation involving transformed arguments of the unknown function. The resulting approach is particularly attractive in the multidimensional setting, where the range of existing numerical methods is relatively limited.

Sloshing is an inevitable phenomenon in an ocean-going vessel that can have adverse effects. In this work, the mitigation of sloshing is investigated using multiple thin porous baffles of various configurations in a partially filled swaying tank. The boundary value problem is solved within the framework of a linearized potential flow theory using the scaled boundary finite element method (SBFEM). The flow through the thin porous baffles is assumed to follow Darcy's law. The computational domain is divided into a minimum number of subdomains due to the presence of porous baffles and to ensure star convexity. Higher-order polynomials are used along each subdomain edge to represent the unknown field, i.e., velocity potential. The developed numerical model is validated with the known results in the literature. Subsequently, various results, such as the amplification factor and the forces of the tank wall, are presented and discussed for the effect of configuration, porosity, slosh tank width, depth of baffle submergence and the space between adjacent baffles. From the parametric study, it is observed that top-mounted baffles enhance sloshing suppression by $50\%$ compared to bottom-mounted vertical baffles, considering all sloshing modes. Assessing the overall effectiveness, top-mounted convex baffle configuration emerges as the most efficient configuration for sloshing suppression, achieving a well-balanced reduction across all modes.

A second order finite difference scheme is analyzed for the nonlocal Cahn-Hilliard system with variable mobility and the singular Flory-Huggins logarithmic potential. The temporal discretization employs a modified Crank-Nicolson formulation, while the mobility is treated explicitly to guarantee ellipticity and reduce computational cost. To strictly enforce the positivity-preserving property at the discrete level, a nonlinear regularization term is added into the numerical scheme. Rigorous analysis shows that the scheme is uniquely solvable and stable under a modified numerical energy. The well-known analytical challenge associated with the variable mobility has to be handled carefully. To overcome this, a convergence framework is constructed based on two non-standard techniques: (1) a higher-order asymptotic expansion (extending up to the third order in time and fourth order in space) to retain a sufficiently order of accuracy; (2) a two-stage error strategy, wherein a rough error estimate first guarantees the discrete boundedness of the mobility and nonlinear terms, followed by a refined error analysis that yields the optimal convergence rate. Numerical experiments validate the theoretical results and demonstrate the robustness of the proposed scheme.

We analyse symmetric interior-penalty discontinuous Galerkin approximations of elliptic stochastic partial differential equations driven by spatial Gaussian white noise on quasi-uniform polytopic meshes. The low regularity of the exact random field precludes the usual SIP consistency argument. We instead restrict the same continuum white-noise functional canonically to the discontinuous polynomial space and separate the error into stochastic forcing-restriction and deterministic discretisation components. For every fixed polynomial degree $k\geq1$, we prove the sharp two-sided estimate \[ \|Y-Y_h\|_{L^2(\Xi;L^2(\Omega))} \simeq h^{2-d/2}, \] giving rates $h$ in two dimensions and $h^{1/2}$ in three dimensions. A finite-dimensional lower bound establishes optimality. The analysis allows polytopic elements with arbitrarily many and arbitrarily small faces, including admissible non-convex elements. In two dimensions we also treat polygons with a single reentrant corner. Although the deterministic approximation exhibits an angle-dependent loss of regularity, the singular resolvent contribution is finite rank, and the sharp first-order stochastic rate is retained. The discontinuous structure yields an exact element-local sampler for the restricted white noise. Numerical experiments on irregular polygonal, non-convex polytopic and reentrant domains support the predicted convergence behaviour.

We consider conforming $\Pone$ and lowest-order Crouzeix–Raviart (CR) approximations of the Dirichlet Poisson problem in two and three dimensions. We obtain computable a priori $L^2$-error bounds on anisotropic simplicial meshes without a global $H^2$-regularity assumption. For piecewise constant data, the Marini relation gives an equilibrated lowest-order Raviart–Thomas (RT) flux from the scalar CR solution. The square of the resulting Marini constant is the largest generalised eigenvalue of a problem involving only scalar conforming and CR matrices. Combined with the elementwise Poincaré constant, it controls both the conforming and CR errors. The CR $L^2$ estimate uses positivity of the discrete CR–conforming gap operator rather than a nonconforming Aubin–Nitsche argument. We also derive an exact decomposition of the Marini defect into the CR–conforming energy gap and an explicit geometric term with no aspect-ratio factor. On the algebraically graded L-shaped mesh family considered here, $q>3/2$ yields $\kappa_{M,h}=O(h)$ and computable $O(h^2)$ $L^2$-operator bounds.

We study randomized constructions of quasi-uniform designs on a compact metric-measure space satisfying two-sided polynomial ball-growth conditions. Given integers $N\le M$, we first draw $M$ independent candidate points and then retain the first $N$ points of a farthest-point traversal of the candidate set. We prove probabilistic non-asymptotic bounds on the mesh ratio and show that $M=\Theta(N\log N)$ is the sharp order of the candidate-pool size required for bounded mesh ratio: a sufficiently large multiple of $N\log(N/\delta)$ suffices with probability at least $1-\delta$, whereas, if $M=o(N\log N)$, the mesh ratio diverges in probability for every procedure that selects $N$ points from the same independent candidate pool. If $M/(N\log M)\to\infty$, the upper bound on the mesh ratio for the exact farthest-point sampling (FPS) tends to $2$. We also quantify the effect of approximate FPS and analyze direct and fast implementations in spaces of bounded doubling dimension. Finally, using a single infinite stream of candidate points and increasing candidate budgets, we construct an almost surely quasi-uniform nested sequence; under supercritical oversampling and exact FPS, or more generally when the approximation factors tend to $1$, its mesh-ratio limit superior equals $2$ for compact positive-volume subsets of Euclidean space.

We study the numerical approximation of high-dimensional integrals using randomized quasi-Monte Carlo (RQMC) methods, with a focus on scrambled Sobol' sequences. While asymptotically faster than Monte Carlo, classical RQMC suffers from error bounds that grow exponentially in the dimension $s$, and its root mean squared error (RMSE) convergence rate is generally no better than $O(N^{-3/2})$. To overcome these limitations, we combine two recent developments: Suzuki's coarse scrambling and the median trick. Coarse scrambling randomizes Sobol' sequences according to their generating base polynomials and significantly reduces the maximal gain coefficient. By appropriately choosing the base polynomials, we show that the coefficient can be made uniformly bounded in $s$, yielding an $O(N^{-1/2})$ RMSE for $L^2$ integrands with a dimension-independent constant. The median trick then enables near-optimal convergence for function classes beyond $L^2$: for $L^p$ integrands with $p\in(1,2)$, the median of independent coarsely scrambled estimates achieves an error of $O(N^{-1+1/p})$ with high probability; for integrands in the Haar wavelet space $\mathcal H_{\mathrm{wav},\alpha,s,p,q}$ with $\alpha_p:=\alpha-(1/p-1/2)_+>0$, we prove a high-probability error bound of $O(N^{-\alpha_p-1/2+\varepsilon})$ for any $\varepsilon>0$; the bound improves to $O(N^{-r-\alpha_p-1/2+\varepsilon})$ when the integrand has dominating mixed derivatives of order $r\ge1$ that belong to $\mathcal H_{\mathrm{wav},\alpha,s,p,q}$. The latter two bounds are uniform in $s$ under suitable conditions on the ANOVA components of the integrands. Numerical experiments confirm the predicted convergence rates and demonstrate the robustness of the proposed approach.

We introduce a novel structure-preserving numerical method for a first-order hyperbolic approximation system, which approximates the solutions of general fourth-order partial differential equations. By employing an implicit-explicit (IMEX) splitting between the stiff and non-stiff terms, we rigorously prove that the proposed scheme is energy-consistent and positivity-preserving under a CFL-type condition. Furthermore, we show that the method remains robustly stable in asymptotic regimes, with a time-step restriction that is entirely independent of the relaxation parameters. Finally, we present a series of numerical examples for thin film equations to validate the theoretical properties and efficacy of the scheme.

We consider a Volterra integral equation with multiscale and nonlinear exponent kernel. We propose the high-order smoothing conditions, under which the initial singularity of the solutions can be eliminated up to any prescribed order. This indicates the application of the multiscale nature of the kernel on local modification of the solutions. Then a discontinuous high-order collocation method with arbitrary polynomial degree is developed and analyzed on uniform or graded meshes based on the solution regularity. This work serves as a comprehensive extension and complement to [Zheng, Qiu and Stynes, SIAM J. Numer. Anal., to appear] in both mathematical and numerical aspects.

We develop an efficient $C^1$ Bernstein quasi-Trefftz discretization for heterogeneous high-frequency Helmholtz problems on unstructured triangular meshes. The method first imposes strong $C^1$ continuity through Bernstein–Bézier smoothness relations on local macro-patches and then reduces the resulting conforming polynomial space by enforcing projected Helmholtz residual constraints. On two-triangle patches, this compresses the local dimension from quadratic growth in the polynomial degree to a trace-sized space of dimension $2p+1$ under the natural rank condition. Variable matrix-valued coefficients are handled through local polynomial projection, allowing coefficient-approximation effects to be separated from discretization error. A graph-residual formulation couples the reduced patch spaces and produces a sparse global system in compressed coordinates. The implementation combines explicit $C^1$ continuation, batched coefficient projection, batched residual construction, and stable local kernel extraction by QR factorization. Numerical experiments on heterogeneous unstructured meshes demonstrate high-order accuracy, roundoff-level $C^1$ conformity, robust high-frequency resolution, and substantial reductions in local and global computational cost. Turning-point and penetrable-scattering examples further illustrate the flexibility of the approach.

Elliptic homogenization is used to determine coarse-grained properties of materials with features on small scales. When these small scale features have rapid, periodic fluctuations, the solution field corresponding to a homogenized constitutive relation closely resembles the true solution based on the heterogeneous material. This homogenized behavior of the material is computed from a cell problem, where a cell is defined to be one period of the fluctuating material. In the context of linear elliptic partial differential equations, the homogenized constitutive relation is defined simply by a constant coefficient tensor, but for nonlinear problems, the homogenized response depends on the macroscopic state and/or its gradient, thus requiring solutions to parametric cell problems. When computing a numerical solution with the homogenized constitutive relation, it is useful to have a differentiable representation of the solution to the cell problem, as derivatives of the homogenized constitutive relation are required in Newton iterations for the macroscopic state field. In this work, we use the Deep Ritz method to solve the parametric cell problems that arise from nonlinear homogenization. First, we exploit the variational structure of the cell problem, then we discretize the dependence of the cell response on both space and the macroscopic state with a neural network. Enforcing boundary conditions on the cell response strongly, we next use the parametric Deep Ritz method to simultaneously solve the cell problem over a range of macroscopic states. We show that this method is accurate, efficient, and offers a continuous and differentiable representation of the cell response over the macroscopic state and gradient. We then show that our parametric representation of the cell response significantly expedites macroscale solutions when compared to a traditional $\text{FE}^2$ scheme.

Thermoacoustic instability in premixed flames results from the coupling between the flame's heat release, as determined by combustion parameters, and surrounding acoustics, as determined by combustor geometry. A primary instability results in flame flattening as intrinsic flame instability modes are stabilised. Secondary thermoacoustic instability results in a parametric flame instability and drastic growth of acoustic amplitudes. Due to their relative expense, numerical simulations of these phenomena remain scarce. In this work, Direct Numerical Simulations (DNS) of thermoacoustically unstable idealised premixed flames in a tube with acoustically closed upstream and open downstream ends are presented. Results herein demonstrate nonlinear saturation of the primary instability as the flame flattens as well as an oscillating flame fingering characteristic of the unsteady Rayleigh-Taylor effect. To reduce computational cost, we perform DNS only on the region surrounding the flame. Acoustics at in- and outflows are described using the Navier-Stokes Characteristic Boundary Condition (NSCBC) method to model their delayed reentry into the domain in a formulation referred to as the Acoustic Delay Characteristic Boundary Condition (ADCBC) method. A new Averaged Proportional and Integral Linear Relaxation (APILR) method is also introduced, which modifies the Classic Linear Relaxation (CLR) method to maintain time-averaged values of inflow velocity and outflow pressure. Here, an integral control term is used to remove non-zero equilibrium time-averaged inflow velocities which impinge control over flame position. Both new methods demonstrate their capability in inert and counterflow flames test cases. These methods enable the numerical simulation of combustion instabilities at significantly reduced computational expense.

We consider the problem of minimizing error in quantized matrix multiplication $C=AB$. Scalar quantization of the factors introduces rounding errors whose scale depends on the maximum absolute entries – the ranges – of their rows and columns. These ranges determine the quantization grid steps. To reduce the error, we optimize over product-preserving transformations that alter the factor ranges and grid steps without changing $C$. Specifically, we seek the smallest leading expected squared error over invertible inner changes of basis and orthogonal outer rotations. Under independent, zero-mean subtractive dither noise on an unbounded lattice, we prove the output-only bound $E_{\rm lead} \ge (c_A+c_B)/K \Vert AB\Vert_*^2$, where $K$ is the inner dimension, $c_A$ and $c_B$ are normalized noise variances, and $\Vert AB\Vert_*$ is the nuclear norm. The bound is tight: an SVD-aligned Hadamard construction attains the infimum whenever a Hadamard matrix of order $K$ exists, including every power of two, while an SVD-aligned DCT construction is within a factor of two for every $K$. Without outer rotations, Gram-matrix balancing minimizes factorization energy, and finite-set flattening achieves the bound within $C\log(K(m+n))$. For power-of-two $K$, conditional expectations deterministically select the Hadamard signs in $O((m+n)K^2)$ exact-real operations. Synthetic experiments verify both constructions and illustrate the tradeoff between regularization and conditioning. These results characterize the full-gauge optimum and quantify the cost of preserving row and column indices.

We investigate the best $L_2$ approximation of mixed Sobolev spaces by shallow neural networks with $n$ neurons and general activation functions. We first establish an activation-independent Fourier-block principle: if an activation has univariate approximation order $\rho$ in the sense of the Fourier-block property, then the global approximation rate has algebraic order $\min\{\alpha,\rho\}$ for target functions of mixed smoothness $\alpha$, up to explicit logarithmic factors. To verify this property for concrete activations, we introduce a structured univariate approximation condition that implies the Fourier-block property with explicit parameters. For $\mathrm{ReLU}^k$, a matching algebraic lower bound identifies $\min\{\alpha,k+1\}$ as the optimal algebraic approximation exponent in any dimension, up to logarithmic factors in the upper bound. The framework also yields the exponent $\min\{\alpha,k+1\}$ for cardinal B-splines and soft-$\mathrm{ReLU}^k$, and the full mixed-smoothness exponent $\alpha$ for ELU and cosine activations, again up to logarithmic~factors.

We develop an energy-optimal generalized scalar auxiliary variable (EOP-GSAV) framework for nonlinear index-one port-Hamiltonian differential-algebraic equations (pH-DAEs). Exploiting the port-Hamiltonian structure, we separate the nonlinear effort, interconnection, and dissipation terms from a constant implicit core. The resulting BDF-1 and BDF-2 schemes require one linear solve per time step with a reusable factorization, while retaining discrete passivity and accurate tracking of the Hamiltonian. The schemes are compared with the implicit midpoint method equipped with full, modified, and frozen-Jacobian Newton iterations. Numerical experiments ranging from a strongly state-dependent nonlinear stress test to large-scale benchmarks demonstrate robust and competitive performance, with substantial efficiency gains in matched-accuracy regimes. A comparison with SUNDIALS IDA shows comparable work-precision behavior at equal order despite a non-specialized Python/SciPy implementation, while unrestricted variable-order adaptive IDA is faster in the high-accuracy regime.