Papers

We study the relation between equal-weight circle quadrature, Gauss–Chebyshev quadrature, and planar Minkowski reconstruction. The regular $N$-point grid is an $\mathbb S^1$-design of strength $N-1$; a half-step rotation of the $2M$-point grid projects to the classical $M$-point Gauss–Chebyshev rule. A common family of quadrature-error functionals gives exact identities for integration error, mixed-area error, and support-function error. For $N\ge3$, the centered Minkowski polygon $P_N$ with perimeter $2\pi$ has the exact Hausdorff distance from the unit disk $B^2$ \[ d_H(P_N,B^2)=1-\frac{\pi}{N}\cot\frac{\pi}{N} =\frac{\pi^2}{3N^2}+O(N^{-4}). \] With Steiner centering, we obtain $d_H(K_\mu,K_\nu)\le\pi W_1(\mu,\nu)$ for balanced, nondegenerate probability measures, where $W_1$ uses geodesic distance. Its exponent is sharp even for smooth positive densities near the uniform measure. The identity $W_1(\nu_N,\sigma)=\pi/(2N)$ yields a general $O(N^{-1})$ bound for the regular grid. The sharper $O(N^{-2})$ rate follows from the sparse discrepancy spectrum and the inverse multiplier $(1-k^2)^{-1}$. Finally, a perturbation estimate gives sufficient conditions for preserving the quadratic rate and its leading constant.

For $0<s\leq 1/2$ and $d\geq 1$, the $\widetilde H^s$-projection onto piecewise linear finite element spaces on quasi-uniform meshes of a bounded Lipschitz domain in $\mathbb{R}^d$ with a uniform exterior ball condition is almost stable in the maximum norm. The mesh domain need not coincide with the original domain.

We consider two-level hybrid Schwarz domain-decomposition preconditioners for finite element discretizations of the high-frequency Helmholtz equation (with wavenumber $k$). To date, all existing $k$-explicit theory for these preconditioners considers coarse spaces in the $\textit{asymptotic regime}$, i.e., where the Galerkin solution is $k$-uniformly quasi-optimal. We consider the situation where the fine and coarse spaces consist, respectively, of piecewise polynomials of degree $p$ on meshes of width $h$ and $H_c$. We prove results for the coarse space in the $\textit{preasymptotic regime}$, under the familiar mesh threshold that $(H_c k)^{2p} \rho(k)$ is sufficiently small, where $\rho(k)$ is the norm of the solution operator (such that $\rho(k)\sim k$ for the Helmholtz equation in free space). We prove that, in this regime, the fixed-point iteration with this preconditioner converges in a number of iterations that grows at most like $\log k$. We give numerical experiments demonstrating this theory. Furthermore, these experiments indicate that our theoretical results are sharp, in that the number of iterations required for these preconditioners to converge with $(H_c k)^{2p} \rho(k) \gg 1$ grows rapidly with $k$.

We propose and analyze two second-order basis-update & Galerkin (BUG) time integration methods for dynamical low-rank approximation on Tucker tensors and tree tensor networks. Both are built from a sequential pre-computation sweep, after which the differential equations for all basis matrices and connecting tensors are evolved by a Galerkin method. Further, the proposed methods are rank-adaptive by construction. The first, the second-order parallel BUG integrator, solves all these differential equations fully in parallel, which is favorable on parallel architectures, followed by a sequential augmentation and truncation step. The second, the second-order augmented BUG integrator, gives up full parallelism but conserves norm and energy for Schrödinger equations and dissipates energy for gradient flows up to the truncation tolerance. For both integrators, we prove a second-order error bound that is robust with respect to small singular values of the matricization of connecting tensors. Numerical experiments for radiative transfer and quantum spin systems validate the theoretical findings.

We study uniform approximation for operator learning on compact subsets of Banach spaces. We bound the recovery dimensions and the number of network parameters in terms of compact-set complexity, the target gauge, and recovery stability. For bounded linear operators with a separable domain, we characterize gauge approximability by coordinate decompositions with a common vanishing gauge and arbitrarily small additive errors. Their partial sums give a fixed family of nonlinear encoder–decoder pairs. We prove that nonlinear recovery with a common vanishing gauge and a rank bound depending only on the model dimension is equivalent to the uniform approximation property (UAP). We also obtain linear recovery estimates with the same rank bound and an operator version for the linear factors of nonlinear targets. For Barron operators taking values in a Banach space $Y$ of type $p\in(1,2]$ and admitting a Bochner-integrable density, sampling gives the rate $N^{-(1-1/p)}$ for the expected $L^p(\lambda;Y)$ error, where $\lambda$ is a probability measure on the input compact set. For fixed finite-dimensional factors and a class of Hölder maps, we obtain matching upper and lower parameter bounds at fixed sufficiently large depth. Classical sampling and truncation examples provide explicit recovery maps, and a nonlinear integral operator illustrates the parameter bounds.

We develop a localized hybrid finite element–random feature method (FEM–RFM) for elliptic interface problems on unfitted background meshes. The central idea is to use the random feature method (RFM) only in the narrow band formed by interface-cut cells of a uniform Cartesian background mesh,while retaining standard $Q_p$ finite elements on the remaining uncut cells. Within the band, local RFM representations are used on the two physical sides of the interface. The prescribed jumps in the solution and normal flux are incorporated through residual equations involving the two RFM representations on the physical interface $\Gamma$, while the finite element method (FEM) and RFM components are coupled across a separate grid-aligned artificial interface $\Gamma_c$ through trace and variational normal-flux residuals. A Schur reduction eliminates the interior finite element unknowns through sparse solves while retaining their Galerkin equations as equality constraints. The remaining band and interface residuals are minimized in a reduced least-squares problem involving only the coupling nodal values and random-feature coefficients. The method is assessed numerically on a range of two- and three-dimensional interface problems involving high coefficient contrasts, nonsmooth and multiple interfaces, and irregular outer boundaries, together with an application-oriented thermal-homogenization problem for a periodic bicontinuous composite.

We propose a symmetric Nitsche method for the biharmonic Dirichlet problem on a smooth domain whose boundary may cut arbitrarily through the computational mesh. Bulk least-squares stabilization in a boundary strip and tangential boundary penalties control the consistency terms through an explicit lifting of the two Dirichlet traces. This gives coercivity on $H^4(\Omega)$ without inverse inequalities and permits a fixed penalty independent of the cut configuration for sufficiently thin strips. Under suitable regularity assumptions, splines of degree $p\geq4$ give energy and $L^2$ error bounds of orders $h^{p-1}$ and $h^{p+1}$. We also analyze a face-stabilized extension to $C^1$ splines of degree $p=2$ and $C^2$ splines of degree $p=3$. Exact-ellipse experiments examine convergence, cut stability, and parameter sensitivity.

In this paper, we establish a quantitative finite-time approximation of scattering fields for nonlinear Alfvén waves in three-dimensional ideal incompressible magnetohydrodynamics near a strong constant magnetic background. Using the global weighted estimates, we prove that truncating the scattering profiles at time $T$ introduces an error in both the $L^\infty$ norm and the $H^{N_*+1}$ norm of order $$ O\bigl(\varepsilon^2(R+T)^{-\delta}\bigr), $$ where $\varepsilon\in (0,1)$ measures the size of the initial perturbation, $R\ge100$ is the scale parameter in the weighted energy, $\delta\in(0,2/3)$ is the exponent related to weights, and $N_*+1$ is the Sobolev regularity order assumed for the initial data. The pointwise estimate follows from the decay of the nonlocal pressure along nonlinear characteristics, while the Sobolev estimate combines higher-order weighted pressure bounds with estimates for derivatives of the characteristic flows. Numerical experiments measure pointwise and Sobolev pressure-integral tails to verify the quadratic dependence of the initial perturbation amplitude.

This paper develops a unified framework for establishing original energy dissipation of equilibrium-preserving exponential time differencing Runge–Kutta (ETDRK) methods for $L^2$ gradient flows. The analysis is based on an auxiliary continuous evolution representation of a general ETDRK method, which allows the energy variation to be controlled without imposing positive-definiteness conditions on the Butcher coefficients or introducing stabilization terms. For a broad class of equilibrium-preserving ETDRK schemes, we derive an explicitly computable time-step restriction. Explicit restrictions are obtained for representative schemes of orders two through five, and the framework also applies to arbitrarily high-order interpolation-based ETDRK methods. For the latter class, the resulting admissible step sizes are substantially larger than those obtained from the previous interpolation-specific analysis. Numerical experiments for the Allen–Cahn equation confirm the theoretical dissipation restrictions.

Operator learning aims to learn mappings between function spaces and has been widely used for PDE-related problems. This raises two fundamental questions. First, is there a connection between standard numerical methods and neural operators, given that both are designed to approximate PDE solution operators? Second, when a neural operator is trained on numerically generated data, does it approximate the underlying infinite dimension operator, the training data themselves, or the discretized numerical scheme used to generate them? A key hypothesis in the paper is that the trained neural operator may mimic the discretized numerical schemes that generate the training data. The hypothesis has been partially and preliminary verified by the linear approximation principle proposed in Chen Chen 1995 paper and the following scaling law papers. However, neural operator approximation from the perspective of numerical discretization remains unexplored. In this work, we consider a class of PDEs and derive neural operator architectures from the numerical schemes. Firstly, the numerical basis functions, such as finite element basis functions, can be approximated or directly constructed by learnable neural networks. Then the finite difference schemes for temporal discretization can be unrolled into iterative neural architectures that generate the coefficients of the basis of the resulting neural operator. By leveraging numerical analysis, we establish convergence error estimates for these numerically derived neural operators. Notably, the network size depends polynomial on the approximation error in $L_2$, improving the nested exponential dependence derived in the general neural operator approximation and aligning with the numerical observations. Finally, numerical experiments demonstrate reduced errors and improved parameter efficiency.

The development of accurate and efficient computational methods for fluid-structure interaction (FSI) problems remains an active research area due to its broad engineering and biomedical applications. Partitioned loosely coupled (LC) schemes are attractive for large-scale simulations because of their modularity and low computational cost. However, their applicability is limited by stability issues in high added-mass regimes, where fluid and structure densities are comparable, as in hemodynamics. In our recent work, we proposed an LC algorithm based on Dirichlet-Neumann interface coupling, the LC-DN-{\alpha} scheme, which improves the stability of the coupling through a stabilization parameter {\alpha}. Although theoretical stability has been demonstrated under suitable constraints on {\alpha}, its optimal value depends on the physical and numerical characteristics of the FSI problem and is difficult to determine in practice. In this work, we propose a data-driven strategy to automatically estimate the optimal value of {\alpha}. We develop a feed-forward neural network to predict {\alpha} from the descriptors of the FSI problem. The model is trained on a dataset of baseline FSI problems in large added-mass regimes and evaluated on more realistic configurations to assess robustness and generalization. Numerical results show that the proposed model generalizes across the considered FSI scenarios and provides robust estimates of {\alpha} under moderate variations in the problem setup, including changes in geometry and loading conditions, that do not completely alter the interaction mechanics represented in the training set. The methodology is intentionally developed and validated on simplified FSI problems, providing a proof of concept for data-driven strategies to improve the efficiency of partitioned FSI simulations and supporting their extension to more complex and realistic scenarios.

We introduce multilevel Picard (MLP) approximations of McKean-Vlasov stochastic differential equations (SDEs) with jumps of finite or infinite activity. Under Lipschitz and integrability assumptions on the coefficients, we show that the MLP algorithm does not suffer from the curse of dimensionality when approximating the solution of the McKean-Vlasov SDE. The latter means that its computational cost grows at most polynomially in both the state-space dimension of the SDE and the reciprocal of the prescribed error tolerance. In two numerical experiments, we demonstrate the practical applicability of the MLP algorithm for two different McKean-Vlasov SDEs in dimensions up to 200.

Stagnation in the summation of floating-point numbers occurs when the running sum, or the accumulator, does not change its value despite addends being greater than zero. Here we define different types of stagnation - commonly referred to as swamping or absorption, where one operand in an addition is nearly or completely rounded off—in floating-point summation, and propose algorithms for detecting it in numerical code that uses floating-point arithmetic. Whilst cancellation in floating-point arithmetic has been studied, for example by Lam et al. in `Dynamic floating point cancellation detection', stagnation has not been exploited as a property for detecting and reporting accuracy loss in floating-point operations. We have built prototype dynamic stagnation detection tools in two languages, MATLAB and Julia, and applied them to applications ranging from basic summation of a vector, to matrix multiplication, to a fluid flow simulation model. With this tool, users are able to identify areas of their floating-point code where stagnation occurs - in those areas of the code one may need to switch to a higher precision number format, apply stochastic rounding, or use different summation algorithms such as compensated summation. To the best of our knowledge we provide a first rigorous definition of stagnation in floating-point summation.

We develop a regularized finite element method for the fractional conductivity equation, recovering the conductivity on an exterior observation region from noisy Dirichlet-to-Neumann measurements and using this reconstruction to determine the interior conductivity. Mesh-localized inputs and mass-lumped Tikhonov regularization discretize the exterior determination principle of Covi–Railo–Zimmermann (2026, Calc. Var. PDE). We prove exterior \(L^2\)- and \(L^\infty\)-convergence with error estimates, including the logarithmic correction at the critical Hölder exponent. The exterior reconstruction enables the Liouville reduction to the Calderón problem for the fractional Schrödinger equation. Separated measurements recover the interior potential, after which overlapping measurements determine the exterior potential tail. Together with the recovered exterior conductivity, the tail supplies the exterior Cauchy data for a Tikhonov-regularized finite element least-squares continuation. We prove \(H^s\)-convergence of the transformed conductivity and \(L^2\)-convergence of the interior conductivity, together with logarithmic state-bias and conductivity error estimates for full observation spaces. The admissible reconstructions remain uniformly positive and bounded. We present one- and two-dimensional numerical experiments illustrating the accuracy of the exterior and interior reconstructions and their sensitivity to measurement noise and discretization.

Gaussian elimination with partial pivoting (GEPP) remains the most widely used solver for dense linear systems $A \mathbf x = \mathbf b$ for $A \in \mathbb C^{n\times n}$. We study the permutation $\pi = \pi(A)$ that arises in the GEPP factorization $PA = LU$, encoded by the permutation matrix factor $P = P_\pi$. When the input matrix is random, so is $\pi$. For random scalar butterfly matrices of size $2^n$ (a recursively defined family originally introduced to eliminate the need for pivoting altogether), we give the exact GEPP factorization and fully classify the induced permutation as an element of a $2$-Sylow subgroup of $S_{2^n}$ contained in the separable permutations. Moreover, the uniform-angle model induces the uniform distribution on this subgroup. For the GOE, GUE, and iid Bernoulli models, the induced permutation is never exactly uniform for $n \ge 2$. We give the precise rate of departure from uniformity at the leading pivot for the GOE and GUE, and give evidence that this non-uniformity vanishes asymptotically in the permuton sense. In contrast, for banded random matrices of sublinear bandwidth, including the tridiagonal $\beta$-Hermite ensembles, the induced permutation converges to the diagonal permuton. We further show that the resulting pivot probabilities are sensitive to implementation choices: standard LAPACK routines compare complex pivot candidates using the $\ell^1$ rather than $\ell^2$ norm, changing the GUE(2) pivot probability from $1/\sqrt3$ to $2/3$. Together these results establish a new connection between random matrix theory and permutation combinatorics through numerical linear algebra.

The suitable moments of the Maxwell velocity distribution function recover the conserved variable and flux vectors of the Euler equations of gas dynamics. However, the Maxwellian has infinite support in velocity space, which poses challenges for the numerical analysis of kinetic schemes based on it. Although Maxwellian is a natural choice, it is not the only distribution function that preserves these key moments. In his seminal paper, Perthame (1990) introduced a compactly supported hat function as an alternative. This paper presents the formulation and analysis of a contact discontinuity capturing Boltzmann scheme based on peculiar velocity, constructed using the hat function. Owing to the compact support of the hat function in velocity space, the numerical analysis of the proposed scheme is considerably simplified, allowing for a direct proof of its positivity. The scheme, extended to second-order accuracy without compromising the positivity-preserving property of its first-order counterpart, is tested on several benchmark inviscid flow problems in one- and two-dimensions.

We extend the BDF2-type finite-element integrator proposed in [M. Aldé, M. Feischl, D. Praetorius; BDF2-type integrator for Landau-Lifshitz-Gilbert equation in micromagnetics: unconditional weak convergence to weak solutions. (Math. Comp., 2026, doi:10.1090/mcom/4210)] to the Landau-Lifshitz-Gilbert (LLG) equation with lower-order effective-field contributions and current-induced non-conservative torques arising in micromagnetic applications. To incorporate these additional terms while retaining a linear solve at each time step, we employ an implicit-explicit (IMEX) time splitting. Our main theorem proves unconditional weak convergence of the fully discrete approximations towards weak solutions of the LLG equation. The analysis removes the mild but artificial CFL-type restrictions required at the first and final time steps in [Aldé M., Feischl M., Praetorius D.; Math. Comp., 2026]. Furthermore, in the presence of the Dzyaloshinskii-Moriya interaction (DMI), the convergence is likewise unconditional and does not require the CFL-type restriction imposed in [G. Hrkac, C.-M. Pfeiler, D. Praetorius, M. Ruggeri, A. Segatti, and B. Stiftner. Convergent tangent plane integrators for the simulation of chiral magnetic skyrmion dynamics. (Adv. Comput. Math., 2019. doi: 10.1007/s10444- 019-09667-z)]. Numerical experiments support the theoretical findings.

We propose a neural algorithm for sampling from distributions specified by unnormalized Boltzmann densities. Our approach is based on the Jordan–Kinderlehrer–Otto scheme for the Kullback–Leibler divergence in the Wasserstein–Fisher–Rao geometry (WFR JKO scheme). Our contributions are twofold. First, we prove that, for any fixed step size, the exact WFR JKO iterates converge exponentially fast to the target as the number of iterations tends to infinity. Notably, this result requires no structural assumptions on the target, such as log-concavity or a logarithmic Sobolev inequality. Second, we develop a neural implementation of the WFR JKO scheme that parametrizes its transport and reaction components using reweighted normalizing flows. Numerical experiments on challenging multimodal targets demonstrate the promising performance of the proposed method.

The Fourier Neural Operator (FNO) learns solution operators of partial differential equations (PDEs) through Fourier-space kernel parameterization, but frequency truncation can limit the learning of high-frequency variations. AM-FNO and SirenFNO generate kernels for all grid modes from spectral coordinates using shared networks, making coordinate encoding and generator design important. Recent work on implicit neural representations (INRs) has proposed constructing frequency interactions through explicit feature composition rather than relying on subsequent MLPs to form them implicitly. Building on this approach, we propose CAFE+FNO, which incorporates Content-Aware Frequency Encoding+ (CAFE+) into Fourier kernel generation. CAFE+ combines Fourier–Chebyshev features through parallel affine branches and a Hadamard product, forming interactions within and across the two feature families. A kernel MLP maps the resulting representation of each normalized spectral coordinate to a complex channel-mixing matrix. Each layer shares its generator across all stored modes, making the number of trainable parameters independent of the number of modes for a fixed architecture. We compare CAFE+FNO with existing FNO variants on five PDE benchmarks and conduct ablation studies on basis configuration, multiplicative composition, and bandwidth learnability. Code and experimental configurations are available at https://github.com/fabsk101/CAFEPlusFNO.git.

We propose a novel unfitted finite element method based on polynomial splines over hierarchical T-meshes (PHT-splines) for elliptic partial differential equations (PDEs) posed on domains whose boundaries are described implicitly by a level-set function. Exploiting the enhanced smoothness of the spline space, we formulate the method by minimizing the $L^2$ residual of the PDE over an extended neighborhood of the physical domain. This formulation allows the method to be implemented entirely on the background mesh using standard element-wise quadrature. It yields a symmetric positive-definite linear system and avoids the geometric construction of curved cut subregions in boundary elements, as well as the associated nontrivial numerical integration. We establish well-posedness of the numerical scheme and derive optimal a priori error estimates in the $H^1$ norm. In addition, we employ a residual-based a posteriori error indicator to guide local adaptive refinement of the hierarchical T-meshes. Numerical experiments corroborate the theoretical results and demonstrate the accuracy and efficiency of the proposed approach.

In the framework of multiple orthogonal polynomials (MOPs) and their extension to the multivariate setting, a unified mechanism for generating structured families of multiple orthogonal polynomials across diverse multidimensional domains is still lacking. In this work, we bridge this gap by developing an extended Koornwinder-type methodology for constructing multivariate multiple orthogonal polynomials. We introduce two distinct coupling schemes: combining a univariate MOPs family with a standard orthogonal family, which enables the definition of both Type I and Type II bivariate systems along with their dual biorthogonality relations, and coupling two univariate MOPs families to form Type II multivariate systems. Using this general framework, we provide the first explicit formulations of multiple orthogonal polynomials on a variety of bivariate regions, including bounded domains such as a parabolic domain and the square $[0,1]^2$, non-standard unbounded geometries like the positive quadrant $(\mathbb{R}_0^+)^2$ and the wedge $\mathbb{V}^2$, as well as new constructions on the triangle $T$. Furthermore, explicit multiple orthogonal systems on the $d$-dimensional simplex $T^d$ and on the $d$-dimensional cone $\mathbb{V}^d$ are given.

We present a new class of structure-preserving, arbitrarily high-order accurate space-time schemes using unstructured simplicial meshes in space for a class of nonlinear acoustic wave propagation systems, in which the nonlinearity is encoded in the constitutive relations between the state and dual variables, while the differential operators retain the standard first-order acoustic wave structure. At the fully discrete level, the schemes satisfy the following properties: I) local conservation of the state variables; II) a local energy/entropy balance on each space-time element, with a consistent, single-valued numerical energy/entropy flux and nonnegative local energy/entropy dissipation, and exact energy conservation when the stabilization terms are omitted; III) weak preservation of the curl involution of the velocity field, even in the presence of discontinuities and with the shock-capturing terms activated. The key ingredients of the new schemes are: i) a discontinuous Galerkin (DG) discretization in space, in which the dual variables are computed as the L2 projections onto the DG space of the energy derivatives with respect to the conserved variables; ii) numerical fluxes and a bubble-weighted, element-local viscous regularization expressed in terms of the dual variables; iii) a continuous-discontinuous Petrov-Galerkin method in time. Properties I)-III) are first proved mathematically and are then verified by suitable numerical tests.

This paper derives an equivalent formulation of the Costabel FEM-BEM coupling for the time-harmonic Maxwell transmission problem with a Lipschitz interface, and constructs a block-diagonal preconditioner for the resulting discrete system. The reformulated coupling is posed on a product of trace spaces and is obtained by introducing a trace lifting operator and the Steklov–Poincaré operator. At the algebraic level, it corresponds to eliminating the interior FEM unknowns, resulting in a discrete system that involves the Schur complement of the interior FEM block. Assuming that the frequency lies outside the resonance sets, we prove that the reformulated coupling is equivalent to the Costabel coupling at both the continuous and discrete levels. Using the operator preconditioning framework, we build a preconditioner from the Galerkin matrix induced by the Maxwell single layer boundary integral operator. The spectral condition number of the preconditioned matrix is shown to be uniformly bounded. Numerical experiments illustrate the practical effectiveness of the proposed approach.

We study linear Lagrange interpolation on exact geodesic spherical triangles obtained by radial projection of chordal triangles. Factorising the element map into an affine map and a radial correction separates chordal anisotropy from radial distortion. We derive exact identities for the surface measure, Dirichlet energy, and directional derivatives, and use them to transfer planar anisotropic interpolation estimates to the sphere. The resulting local error estimates retain the two chordal directional length scales and keep the geometric factors explicit, with constants independent of the element size, shape, and sphere radius. In the gradient estimate, the maximum-angle factor multiplies only the longer directional scale. Four explicit triangle families quantify the loss caused by replacing these scales with the diameter and establish the sharp maximum-angle dependence. They also show that, within the same two-term directional bound with a constant independent of the element, the maximum-angle factor cannot be moved to the shorter scale or have its exponent reduced below one. The fourth family proves the optimality of the square-root radial-transfer factor while retaining the same transported directional scale. The local estimates extend to compatible spherical triangulations. As an application, they yield an anisotropic energy-error estimate for the conforming approximation of the mean-zero Laplace–Beltrami problem with exact integration.

We develop and analyze a locking-free virtual element method for the three-dimensional quasi-static electroporoelasticity equations. The pressure, electric, magnetic, and displacement fields are discretized by the lowest-order nodal, edge, face, and Stokes-like virtual element spaces, respectively, so that the scheme works on general polyhedral meshes, while the temporal discretization is carried out by a backward Euler scheme. The locking phenomenon is overcome by defining the Stokes-like space on a refined mesh $\widetilde{\mathcal{T}}_h$, on which we can construct an interpolation operator $\mathrm{I}_s$ satisfying $\div(\mathrm{I}_s\bm v)|_K=\Pi^K_0\div\bm v$, which renders the error estimates independent of the Lamé constant $\lambda$. Moreover, since the exact sequence $\curl\bm V^e_h=\bm V^f_{h,0}$ holds, the magnetic field is updated explicitly, reducing the four-field system to a three-field one and lowering the computational cost. Numerical experiments are provided to demonstrate the performance of the proposed method.

This paper proposes and analyzes a numerical method for coupled wave-heat systems that arise, for example, in visco-elasticity with memory, thermo-elastic wave propagation with finite thermal speed, and Maxwell's equations coupled to dispersive material laws. Motivated by the widespread use of explicit leapfrog schemes for wave problems and the parabolic time-step restriction for heat equations, we design a second-order scheme that combines leapfrog time integration for the wave part with a Crank–Nicolson step for the heat part. The coupling is arranged in a Strang-splitting–type fashion so that the overall scheme remains explicit in the coupling and is subject only to a CFL condition of hyperbolic type, rather than the more restrictive parabolic constraint. We introduce an abstract space discretization that covers a broad class of wave–heat systems and accommodates both conforming and stabilized non-conforming discretizations. Using a general error decomposition, we derive error estimates for the fully discrete scheme under the desired CFL condition. Numerical experiments for our model applications confirm the theoretical results.

Physics-informed neural PDE solvers adapt their parameters to satisfy governing equations, yet their representational structure typically remains fixed throughout training. This rigidity is poorly matched to PDE solutions with strongly heterogeneous complexity across space and space–time, leaving capacity insufficient where the physics is difficult and redundant where it is simple. We introduce physics-informed neural plasticity, a paradigm in which the representation itself reshapes during optimization in response to unresolved physics. We instantiate this principle with Representation Capacity Adaptation for PDEs (ReCAP), a Gaussian-localized solver that dynamically redistributes capacity through local enrichment, residual-directed splitting, gate-based pruning, and function-aware merging. ReCAP uses responsibility-weighted error indicators and the geometry of residual energy to determine where and how to refine. To limit the disturbance introduced by splitting, we introduce quiet-child refinement, which initializes new components by transporting the parent representation while controlling instantaneous functional perturbation. We further establish conditional a posteriori reliability and structural-stability guarantees linking localized physics residuals to solution error and stable refinement. Across five challenging 3D and 4D PDE benchmarks against 11 physics-informed solvers, ReCAP achieves the lowest relative $L^2$ error on every problem, reducing error by $10.7\%$–$27.5\%$ relative to the strongest competing result. These results suggest that physics-informed solvers need not merely learn their parameters—they can learn how their representational capacity should be organized.

Thermo-mechanical fracture couples transient heat conduction on a cracked domain with a crack that grows as the temperature and the displacement evolve. Neural energy solvers have been proposed for phase-field fracture and later extended to represent a sharp crack through the network input, but heat conduction on the cracked domain and crack propagation under the resulting thermal stresses have not yet been treated together in these solvers. We present an extended deep energy method for thermo-mechanical crack propagation in which the crack remains a sharp polyline. Two networks represent the temperature and the displacement and receive the crack through a scalar embedding function, discontinuous across the crack and smooth elsewhere, so that both fields can jump across it without a regularization length, and the displacement is enriched near the tip by the Williams expansion with trainable amplitudes. The two fields are obtained by minimizing an incremental conduction functional and the thermoelastic potential energy in a staggered sequence, with Monte Carlo integration on points stratified over background elements, densified near the tip and redrawn during training. The stress intensity factors are extracted by the interaction integral with the area term of Wilson and Yu and checked by a sweep of the contour radius, and the crack advances at the maximum hoop stress angle when the energy release rate of the kink reaches the critical value at the crack-tip temperature. On a stationary thermal edge crack the extracted stress intensity factor agrees with the published value to 0.11%, in a functionally graded shear test initiation agrees with an independent sharp-crack finite element solution to within one load step, and on a notched cruciform specimen the crack paths follow the published solutions under mechanical, thermal and combined loading.

Optimal control of open quantum systems governed by the Lindblad master equation requires the repeated solution of forward and adjoint evolution equations. For large Hilbert spaces, these computations demand numerical methods that are structure-preserving, accurate, and computationally efficient, while remaining compatible with nonsmooth optimization techniques. In this work, a unified framework for optimal control of Lindblad equations is presented. The continuous optimality system is derived from the Pontryagin maximum principle and solved by a sequential quadratic Hamiltonian (SQH) method. Its numerical realization is based on second-order exponential midpoint propagators for the forward and adjoint Lindblad equations with time-dependent Hamiltonians. Full-rank schemes preserving the Hermitian and positive-semidefinite structure of the density matrix are developed together with low-rank formulations that substantially reduce storage requirements and computational cost. Rigorous error estimates are established for the full- and low-rank forward and adjoint propagators. Numerical experiments confirm the predicted convergence rates, demonstrate the effectiveness of the low-rank approximations, and illustrate the performance of the proposed FREM-SQH and LREM-SQH algorithms for optimal control problems with smooth and nonsmooth control costs.

We develop an end-to-end framework for estimating canonical thermal properties of particle-number-conserving second-quantized Hamiltonians using quantum simulations. The framework integrates initial state preparation with mixed Dicke states, particle-number-preserving Trotterized simulation, and generalized quantum signal processing to estimate canonical partition functions and observables. Building on this canonical scheme, we introduce a grand-canonical reconstruction scheme in which the same quantum measurement data can be reused for arbitrary chemical potentials through classical reweighting, thereby avoiding repeated quantum circuit executions as the chemical potential is varied. We prove the existence of system-size-independent temperature regimes in which both the canonical and reconstructed grand-canonical schemes have lower asymptotic time complexity than classical Lanczos methods, up to polynomial factors, for physically relevant Hamiltonians. We numerically validate our theoretical results for the two-dimensional Hubbard model in both weak- and intermediate-coupling regimes.