Papers

We develop a boundary element method for the Laplace equation with Dirichlet boundary conditions in a bounded three-dimensional domain. The method uses flat panels of maximum diameter \(h\) on an interior surrogate boundary and corrects the imposed boundary values by a first-order Taylor expansion along the panel normals. The correction uses the normal flux already present in the Calderón formulation and introduces no additional unknown field. Under uniform assumptions on the surrogate geometry and boundary operators, the penalised form is coercive when \(\gamma\delta_h\le c_0<4\), where \(\gamma\) is the penalty parameter and \(\delta_h\) is the maximum normal offset. The error estimate separates trace approximation from a geometry term of size \(\gamma^{1/2}\delta_h^2\|u\|_{H^3(\Omega)}\). With \(\gamma\simeq h^{-1}\), the sufficient condition \(\delta_h\lesssim h^{(k+1)/2}\) yields an error of order \(h^{k+1/2}\) for \(P_k/P_{k-1}\) Dirichlet/flux spaces under the stated regularity and uniform trace bounds. Thus an offset of order \(h^2\) supports the estimate through cubic Dirichlet approximation. We analyse two penalty variants and present numerical experiments on a sphere and an ellipsoid to examine convergence and sensitivity to geometry, penalty scaling and quadrature.

Physics-Informed Neural Networks (PINNs) build neural representations of time-dependent PDE solutions, naturally incorporating physics knowledge and observational data, which makes them well suited to both forward and inverse PDE problems. PINNs, however, are known to suffer from spectral bias and lack of causality. Neuro-Spectral Architectures (NeuSA), a recently proposed alternative to PINNs, mitigate both issues, but their numerical integration becomes unstable for stiff differential equations arising in many relevant physical problems. This study proposes Neuro-Spectral Exponential Time Differencing Architectures (NEXT), which combines the spectral representation of the PDE solution in NeuSA with high-order exponential integrators. Within this approach, the linear stiff part of the vector field induced by the PDE is integrated exactly through matrix exponentials, while the possibly nonlinear remainder is modeled by a neural network. The effectiveness of NEXT is verified through benchmark experiments on a set of stiff PDEs, in which NEXT is stable and accurate while NeuSA diverges numerically. It is also shown that NEXT can be applied to inverse problems, where the model has to learn unknown parameters or boundary conditions from sparse data. All code used in this work is publicly available at: https://github.com/marcioh2m/next.git .

Low-order refined (LOR) preconditioning replaces a high-order operator with a low-order discretisation on a refined nodal mesh. For tensor-product elements, the two operators are spectrally equivalent with bounds independent of the polynomial order $P$, but the construction does not extend directly to simplex and mixed-element discretisations. This work makes two contributions: it extends LOR preconditioning to simplex and mixed-element discretisations, including triangular, tetrahedral, and prismatic elements, and establishes a generalised Vandermonde transformation linking the modal and nodal LOR formulations, showing that the resulting preconditioned spectra and Krylov convergence are independent of the high-order basis. Numerical experiments show controlled iteration growth on triangular meshes despite increasing condition number, and controlled iteration counts up to $P=5$ on tetrahedral, prismatic, and mixed-element meshes. A single algebraic multigrid V-cycle per outer iteration gives the best balance of iteration count and cost. The method is applied to a production incompressible Navier-Stokes simulation of a race-car front-wing and wheel configuration, discretised on a mesh of $2.87\times10^6$ mixed prismatic and tetrahedral elements giving $32.2\times10^6$ pressure degrees of freedom at $P=3$. LOR reduces the mean pressure conjugate gradient (CG) iteration count from $235.3$ to $5.5$, and the pressure-solve time over 1000 timesteps by 16.1%, relative to the default production static-condensation diagonal preconditioner in Nektar++. The one-time cost of constructing the LOR preconditioner is amortised over the production simulation.

This proof-of-concept paper demonstrates how the Kalman filter, as a brain imaging method, can be enhanced by incorporating raw diffusion tensor imaging (DTI) data without tractography, in addition to EEG recordings. An efficient algorithm is designed to streamline the otherwise tedious process of inferring connectivity between brain regions from DTI data. The DTI-based evolution model is applied to the Kalman filter and the Standardized Kalman filter, and these are compared with their conventional counterparts using the random walk evolution model. The numerical experiments use synthetic somatosensory and auditory evoked potentials. The results show that the DTI-modelled Kalman filter can obtain deep activity in the brainstem and thalamus. Moreover, the connection model limits the estimated spreads to the correct or nearest anatomical brain region.

We study a fully discrete approximation theory for the two-dimensional incompressible Euler equations on $\mathbb{T}^{2}$, in vorticity form, driven by Kraichnan-type transport noise and with $\mathbb{H}^{-1}(\mathbb{T}^2)$-valued initial vorticity. The main new ingredient is a finite-dimensional coercivity estimate for the torus Kraichnan model, which captures at the discrete level the regularizing mechanism induced by the rough transport noise. Together with a new compactness strategy based on continuous equation-based interpolants and localized discrete estimates, this coercivity estimate gives tightness of the approximations and shows that subsequential limit is a weak martingale solution of the Euler equations.

The diffusion synthetic accelerated (DSA) source iteration is a standard solver for the radiative transfer equation. Using Fourier analysis, a convergence rate \(\rho_\infty(c)\le0.2247\,c\) with \(c\) the maximum ratio of scattering to total cross section has been established for an infinite homogeneous medium. For slab geometry with inflow boundary conditions and arbitrary bounded cross sections we prove that the spectral radius of the DSA iteration is at most \(\rho_\infty(c)\). We show that this convergence rate carries over to a variational discretization of the DSA iteration on every conforming tensor-product Galerkin space whose angular factor contains the constants. The analysis rests on an exact min–max characterization of the spectral radius. From it we derive a checkable sufficient condition for a given rate. We verify this condition using a suitable energy-stable projection, a weighted angular average whose weight is chosen so that the condition holds with the rate \(\rho_\infty(c)\). Since the projection maps discrete spaces into discrete spaces, the argument applies verbatim to the discrete iteration. As a consequence, the condition number of the preconditioned system is uniformly bounded by \(1+0.29\,c\), which implies rapid convergence of the preconditioned conjugate gradients method.

The present work proposes a new operator splitting technique. The method leverages ideas used in PDE-constrained domain decomposition, adapting them to the objective of decoupling equations defined on the same geometric domain. The methodology is applied here to the resolution of coupled flow and deformation in the linear quasi-static Biot problem. An overview of its main characteristics is provided, also in comparison with the well established fixed-stress split for poroelasticity.

We propose a novel domain decomposition framework based on energy minimization, applicable to a wide class of linear and nonlinear partial differential equations. This energy-minimizing domain decomposition (EMDD) method addresses problems whose solutions are minimizers of an energy functional, either unconstrained or constrained to a manifold, and thus covers source and eigenvalue problems in a single formulation. Besides introducing the general framework, we derive the algorithmic realization of the method for four classes of model problems, namely linear and nonlinear source and eigenvalue problems, and derive the computational complexity of the method. Numerical experiments demonstrate the robustness and efficiency of the method across all four model problems and show that it compares favourably against existing model-specific domain decomposition methods.

In this work we construct a new semi-implicit (SI) hyperbolic and thermodynamically compatible (HTC) scheme for the compressible Euler equations of gasdynamics. The challenge therein lies in solving the entropy inequality as the primary evolution equation instead of the total energy conservation law, which in our framework is an additional conservation law that needs to be fulfilled by the numerical scheme as the consequence of a compatible discretization. The proposed method has a predictor-corrector character, where we first solve the isentropic problem without considering the entropy production. Therein, the entropy and the nonlinear convective terms in the momentum flux are discretized explicitly, while the mass flux and the pressure terms are taken implicitly. The discrete momentum equation is then inserted into the discrete mass conservation equation, leading to a mildly nonlinear system for the pressure, which can be efficiently solved with a nested Newton-type algorithm. Based on this prediction, a non-negative production term in the entropy is added, ensuring a cell entropy inequality for the numerical scheme. To guarantee global conservation of total energy at the discrete level, a new modified global Abgrall-type flux correction is introduced. The presented numerical method allows time steps that are restricted only by the material wave speed and not by the sound speed. This makes the scheme particularly useful for flows in the low-Mach-number regime. Numerical results are shown in one and two space dimensions to assess the theoretical properties of the new scheme.

We develop a balanced truncation method for linear delay differential systems based on a finite-rank alternating-direction implicit (ADI) iteration for the associated operator Lyapunov equations. The method acts directly on the infinite-dimensional history-state realization and does not discretize the delay interval during Gramian approximation or balancing. The finite-rank reachability and observability factors are represented exactly in Takenaka–Malmquist coordinates, reducing the required computations to finite-dimensional linear algebra. These factors yield a finite-dimensional, delay-free reduced model. For real system data, the implementation remains in real arithmetic. Under suitable conditions on the ADI shifts and a gap in the Hankel singular values, we prove convergence in $\mathcal H_\infty$ to the reduced transfer function obtained by exact balanced truncation of the infinite-dimensional system. The construction extends to finitely many discrete delays. Numerical experiments validate the method against independent reference computations. The complete MATLAB implementation used in the experiments is provided as supplementary material.

We propose a vanishing-viscosity minimizing movement (VVMM) Fourier spectral method for scalar hyperbolic conservation laws under periodic boundary conditions. The method is a fully discrete space–time residual-minimizing scheme, Fourier spectral in space and polynomial in time, stabilized by vanishing viscosity and heat-kernel postprocessing. A basic analytical difficulty in Fourier discretizations of nonlinear conservation laws is the interaction between Fourier projection and the nonlinear flux; the resulting projection residual contains unresolved high-frequency components. VVMM addresses this issue by selecting, on each time slab, the numerical solution that minimizes {a penalized cut-off viscous conservation laws residual} and then applies heat-kernel postprocessing to damp unresolved modes. We prove a fully discrete $L^1$ error estimate, in any fixed space dimension, that separates the vanishing-viscosity error, the temporal interpolation error, and the spatial Fourier truncation error. {Under exact minimization of the penalized cut-off residual and suitable parameter coupling,} the estimate yields, for each fixed $\gamma>0$, an $L^1$ convergence rate $N^{-1/2+\gamma}$, with a constant depending on $\gamma$. Numerical experiments in one and two dimensions illustrate stable shock capturing without visible Gibbs oscillations and convergence behavior consistent with the theory.

We investigate the use of neural network solvers for infinity and $p$-Laplace problems, which are fundamental in nonlinear analysis and have practical applications. Our approach employs Physics-Informed Neural Networks (PINNs) and Deep Operator Networks (DeepONets) to address computational challenges associated with large $p$ values, ranging from $2$ to $1000$, on various 2D and 3D domains. Our method offers advantages over traditional physics-based solvers, especially in three dimensions where mesh-based solvers become very costly for these problems. We also establish conditional convergence results for PINN approximations of both problems and a universal approximation result for DeepONet on the parametric $p$-Poisson problem. We demonstrate the effectiveness of these neural network solvers through numerical experiments and compare their performance with conventional methods.

We propose a Sparsely Connected Rank-Inspired Neural Network (SCRINN) for function approximation and the numerical solution of partial differential equations (PDEs). By introducing a structured sparse connectivity architecture and an orthogonality- and normalization-guided loss function, SCRINN constructs an approximately orthonormal neural basis while promoting progressive rank growth of its evaluation matrix and reducing the number of trainable parameters. A two-stage framework is developed, in which an approximately orthonormal neural basis is first constructed, and the output weights are subsequently determined using least-squares methods. Numerical experiments on function approximation and a variety of steady-state, time-dependent, and nonlinear PDEs demonstrate that SCRINN consistently achieves high approximation accuracy with substantially lower training costs. The code is publicly available on GitHub.

We study an elliptic optimal control problem governed by Darcy's equation in heterogeneous porous media, with pointwise box constraints on the control. The objective functional is formulated via a Radon measure, which allows the desired pressure state to be tracked on observation sets of varying dimension, including points, curves, and subdomains, within a single formulation. The state and adjoint equations are discretized by a symmetric interior penalty discontinuous Galerkin method, yielding a locally mass-conservative approximation that is robust across strong permeability discontinuities, while the control is approximated by piecewise constants. The state, adjoint, and control are retained as primary unknowns in a single coupled optimality system and solved monolithically by a primal-dual active set strategy. We establish stability and well-posedness of the discretization and derive a priori $L^2$ error estimates for both variables. The principal difficulty is the reduced regularity of the adjoint state induced by the measure-valued tracking data. The analysis controls the resulting adjoint-control coupling through an intermediate adjoint driven by the continuous optimal state. Numerical experiments confirm the predicted convergence rates for point, curve, and subdomain observations, exhibit mesh-independent primal-dual active set iteration counts, and demonstrate local mass conservation.

We investigate fully nonlinear water waves and fluid-particle dynamics generated by the horizontal motion of a seabed obstacle with prescribed time-dependent velocity. The governing equations are the full Euler equations, formulated in a time-dependent conformal domain that simultaneously maps the moving free surface and seabed onto fixed boundaries. The main contribution of this work is a Lagrangian formulation for computing particle trajectories in the resulting genuinely unsteady conformal domain. We derive a closed-form trajectory system in the canonical domain in which the time dependence of the conformal map is entirely represented by the real and imaginary parts of an analytic function that can be evaluated spectrally from the surface. We apply the formulation to waves generated by horizontal submarine landslide motion and characterize particle displacements throughout the fluid as functions of the initial particle position and Froude number. Our findings identify distinct regions in which particle motion is predominantly associated with the moving seabed, the generated wave, or the combined action of both. The results reveal a transition in the dominant mechanism driving particle motion: at low Froude numbers, particle displacements are primarily associated with the moving seabed, whereas at high Froude numbers the generated wave becomes increasingly dominant, particularly near the free surface and away from the obstacle path.The numerical predictions are benchmarked against laboratory data available in the literature, showing good agreement and providing a quantitative assessment of the model accuracy. Moreover, laboratory topographies beyond those considered here can be readily incorporated into the numerical framework through a Hermite interpolation procedure.

We design and analyze a nonoverlapping additive Schwarz preconditioner for interior penalty discontinuous Galerkin discretizations of anisotropic elliptic problems. The preconditioned method coupled with a Krylov subspace iteration is shown to be independent of the highly discontinuous (and anisotropic) jump coefficients as well as the subdomain size. To increase efficacy, various auxiliary spaces are considered to reduce the size of the coarse grid operator. We demonstrate how to modify the additive Schwarz preconditioner such that it is applicable to the nonsymmetric IPDG schemes. Several numerical experiments verify the theory and validate the robustness of the preconditioner.

Smoothed Particle Hydrodynamics (SPH) is well suited to a range of problems, particularly those involving large deformations in fluid dynamics. In recent years, in addition to the advancement of SPH formulations, the development of differentiable solvers has also progressed. However, unified frameworks that flexibly accommodate diverse numerical schemes, including advanced and implicit methods, while supporting continuous extension and updating remain limited. In this study, a high-accuracy and composable differentiable SPH solver, CraftSPH, has been developed. In CraftSPH, major computational operations are implemented as independent modules, which can be combined according to the intended purpose of constructing a solver. This design enables different computational schemes to be intuitively constructed from combinations of common components and allows newly developed high-accuracy discretization methods to be flexibly incorporated and extended. Furthermore, each module supports automatic differentiation, enabling applications to physical-parameter estimation and to optimization that combines SPH solvers with deep learning models. Accordingly, CraftSPH provides an SPH computational framework in which high-accuracy discretization, explicit and implicit computations, and automatic differentiation can be handled in a unified and extensible manner. To demonstrate the applicability of CraftSPH to a wide range of problems, forward analyses are conducted for Poiseuille flow, two- and three-dimensional dam-break problems, and a rising-bubble problem. In addition, the performance of automatic differentiation is evaluated through inverse problems involving parameter estimation for the Taylor-Green vortex and lid-driven cavity flow.

Calderón preconditioning restores mesh-independent conditioning for boundary-element discretizations of the Maxwell electric field integral equation, but assembling the dual barycentric operator can dominate the computational cost. We show that this preconditioner can be assembled once at a nominal geometry and reused across connectivity-preserving shape perturbations. Exact combinatorial alignment of the RWG-BC degrees of freedom transports the duality pairing, leaving only the dual-mesh operator frozen. Under uniform discrete-stability and geometric operator-continuity assumptions, a bi-parametric operator-preconditioning analysis yields a spectral-condition-number bound uniform in the mesh size and, after tensorization, in the stochastic-Galerkin dimension. Numerical experiments on smooth and non-smooth scatterers show that the frozen preconditioner closely tracks freshly assembled Calderón preconditioners over substantial shape variations. A Monte Carlo campaign on a NASA almond, whose admissibility is screened before each shape is constructed so that no realization is rejected, separates the two strategies by less than the variation between realizations, and mesh refinement there leaves the frozen iteration count unchanged. Reuse yields roughly 30x measured speedup per new realization in our dense implementation; a cost-model estimate based on published compressed-assembly timings gives roughly 2x.

In this paper, we establish the first asymptotic convergence guarantees for the Muon algorithm through a more accurate proxy for the Newton-Schultz iteration than the typical matrix sign function. We prove that, for appropriate choices of hyperparameters, the iterates satisfy $\lim_{k\to\infty}\|\nabla f(x_k)\|=0$, and, under a global Polyak-Łojasiewicz condition, that the sequence of function values converges linearly. The key insight is that the regularization, implicit in Muon's Newton-Schulz implementation, induces a bounded preconditioner, exposing Muon as a preconditioned Polyak heavy-ball method and enabling a classical Lyapunov analysis. This observation naturally motivates applying the same preconditioning structure to the Nesterov gradient evaluation. We formalize this idea by introducing Muesterov, a Nesterov-based variant of Muon, and prove that it enjoys the same convergence guarantees, extending the theoretical framework beyond the heavy-ball setting. Numerical experiments on a scalar cross-entropy problem corroborate the theory and illuminate the joint role of the learning rate and the Newton-Schulz regularizer in controlling convergence. Preliminary numerical simulations training the nanoGPT dataset provide intuition regarding the relevance of the observations in this paper to practical applications.

We address the inverse problem of reconstructing the initial temperature distribution in a one-dimensional non-homogeneous heat equation with Dirichlet boundary conditions from a finite number of pointwise-in-time measurements at a single spatial location. We propose an explicit approximation scheme that incorporates the source term and utilizes a refined sequence of measurement times. Unlike previous approaches that often rely on exponentially growing observation times, our chosen time sequence allows for efficient recovery within a fixed time horizon. We establish algebraic convergence rates under suitable regularity assumptions on the initial data and the forcing, and we validate the theoretical sharpness with numerical experiments for varying numbers of measurements.

High-dimensional extreme eigenvalue problems often arise from molecular vibrational models, electronic structure calculations, and quantum mechanics. Directly solving these problems suffers from the curse of dimensionality. Instead of tackling the problem in high-dimensional ambient space, we reformulate the problem through low-rank tensor formats, which can significantly reduce the computational cost and storage. Specifically, we consider the Rayleigh–Ritz problem on bounded-rank tensors in the tensor train format, which enables a more flexible choice of rank parameters. Moreover, we consider a smooth parametrization for bounded-rank tensors by introducing slack variables, leading to a smooth manifold structure. The original Rayleigh–Ritz problem is therefore transferred to an optimization problem on the manifold. We develop the Riemannian geometry and propose optimization methods for solving extreme eigenvalue problems. In practice, the Kronecker-product structure in Hamiltonian and PDE operators is employed to simplify the computation. Numerical experiments on the harmonic oscillator, the Laplace operator, the Schrödinger equation, the layered cluster problem, and the Bose–Einstein model demonstrate that the proposed method achieves accuracy comparable to full-space eigensolvers with reduced computation time and avoids storing full vectors. The results also show numerical rank reduction during iteration and robustness when the rank parameter is over-estimated. In addition, the fixed-point iterations illustrate that the proposed methods are able to serve as an inner eigensolver for solving nonlinear eigenvalue problems.

In science and engineering, many physical laws and scientific principles naturally arise as variational problems of minimizing an energy functional over an appropriate functional space. In many cases, it is more advantageous to directly discover the minimizer of the energy functional than to solve the associated Euler-Lagrange equations. Conventional numerical solvers, such as the finite-element method (FEM), often lack the flexibility to incorporate prior knowledge or noisy observational data. In recent years, machine learning approaches, particularly kernel-based methods, have attracted growing attention in scientific computing. Unlike FEM, which requires domain discretization and mesh generation, kernel-based methods construct solutions from scattered nodes, thereby avoiding the complexity of meshing, especially in high-dimensional or geometrically complex domains. This works presents a two-step procedure for variational energy minimization based on the optimal recovery formulation in a reproducing kernel Hilbert space (RKHS), providing a unified framework for seamlessly incorporating both physical constraints and noisy data. On the computational side, we employ the sparse Cholesky decomposition for Matérn kernels to alleviate the well-known cubic complexity bottleneck of kernel methods. On the theoretical side, we rigorously establish the existence and convergence theory of the proposed method using $\Gamma$-convergence theory. Numerical experiments on benchmark problems demonstrate the efficiency, robustness and accuracy of the approach.

This paper develops a divergence-free, inf-sup stable, optimally convergent and pressure-robust finite element method for the three-dimensional Stokes problem on curved domains. The geometry is approximated by an isoparametric tetrahedral mesh, while the velocity space is obtained from Scott–Vogelius spaces on an Alfeld split by the Piola transform. The discrete inf-sup condition is proved using suitable face bubble functions. The insufficient accuracy of quadrature rules on curved triangular interfaces leads to a consistency error that may cause suboptimal convergence. The remedy is to introduce a suitable consistency correction without stabilization terms. Moreover, commuting operators on curved domains are constructed from local commuting interpolants and shown to be globally conforming. The operators are used to approximate the load and obtain a pressure-robust discretization. Numerical examples are provided to validate the theoretical results.

We give a constructive numerical treatment of Monge's optimal transport problem in 1D, 2D, 3D, following Evans-Gangbo: the map is built from a Kantorovich potential $u$ and density $a$ solving $-\mathrm{div}(a\nabla u) = f^+-f^-$, via $\dot T = -a\nabla u(T)/[(1-t)f^+(T)+tf^-(T)]$, $T(0,x)=x$. In 1D this gives a closed form. In 2D–left unfinished in an earlier report–we obtain $(u,a)$ via a convex Beckmann flow problem (second-order-cone programming), cross-checked against a dual LP, then build the map by direct integration or a ray-based 1D "clock" ODE with explicit density per ray. Validation uses strong duality, mass conservation, and vanishing-regularization limits, on a translation-invariant case and two 2D examples (a radial cone, a checkerboard) checked by symmetry; both show boundary-concentrated mass and non-unique flux where supports touch. We extend to 3D (theory unchanged), matching the 2D answer to six digits and reproducing the same phenomena on a radial ball, exposing a limit of the 2D visualization style. Movies illustrating the transport maps are available as ancillary files on the arXiv page of this submission.

We consider a class of SPDEs driven by a standard real-valued Brownian motion, with drift and diffusion coefficients such that there exists a unique mild solution taking values in the interval $[-1,1]$ almost surely. To preserve this qualitative property of the exact solution, we propose a domain preserving Lie–Trotter splitting scheme: for any choice of the time-step size, the numerical solution takes values in the interval $[-1,1]$ almost surely. Furthermore, we prove mean-square convergence with rate $1/2-$ for the domain preserving Lie–Trotter scheme. These theoretical results are illustrated with numerical experiments.

The Brinkman penalisation method embeds boundary-value problems on complex domains into a simple computational box by modeling the solid region as a strongly dissipative medium, avoiding body-fitted mesh generation. We show that multi-symplectic Hamiltonian PDEs regularised by Brinkman-type penalisation retain a multi-conformal symplectic structure under a compatibility condition linking the symplectic matrix and the penalisation projection. This yields an exact local conservation law, under which the multi-symplectic two-form is conserved in the fluid region and decays exponentially inside the solid. The linear wave equation with Brinkman friction and Maxwell's equations with artificial Ohmic conductivity satisfy this condition, with explicit modified Hamiltonian densities. Building on this, we propose (i) structure-preserving numerical integrators via Strang splitting that satisfy a discrete conformal conservation law, and (ii) conformal symplectic neural operators that interleave exact dissipative flows with learnable multi-symplectic evolution operators, allowing geometry-dependent operator learning. Numerical experiments on wave and electromagnetic scattering demonstrate that our methods reproduce correct local energy budgets and avoid unphysical energy drift, providing a principled framework for physics-consistent scientific machine learning on complex domains.

In the present work, we study the properties of the sparse entropic quadrature method, a recently proposed method for the closure of systems of moment equations in rarefied gas dynamics. We prove the properties of the method regarding its well-posedness, hyperbolicity, and stability for first-order spatial convection schemes. We apply the method to two model rarefied gas dynamics problems and compare the results with those given by a discrete velocity method and show that it is capable of reproducing rarefied gas effects whilst having significantly fewer degrees of freedom than the discrete velocity method.

We set up and investigate an approximation scheme for a variational theory of cohesive fracture in antiplane shear deformations. We consider a family of functionals acting on scalar valued functions in general dimensions depending on a small parameter $0 < \varepsilon \ll 1$ and on two fields: the elastic part of the displacement field and an eigendeformation field that describes the inelastic response of the material beyond the elastic regime. We measure the inelastic contributions of the latter in terms of a non-local energy functional. Our main result shows that, as $\varepsilon \to 0$, the approximate functionals $\Gamma$-converge to a cohesive zone model.

Geometry-dependent PDEs define operators whose inputs and outputs may each consist of an oriented interface and a function on that interface, whereas most existing neural operators are formulated on fixed domains. In this paper, we first establish an approximation theorem for continuous operators between general state sets. We then instantiate the theorem with the Interface Autoencoder (IAE), which represents interface–function states in a fixed box using finite projection codes and decodes them by zero-set extraction and restriction. For tensor cosine codes, we prove a uniform $C^0$ reconstruction rate and tube-local gradient convergence. These estimates yield Hausdorff and graph-Hausdorff error bounds, together with a reference-based topology certificate. The IAE accommodates both cross-space geometric operators, mapping domain geometries to the associated PDE solution fields, and within-space geometric operators, mapping between interface–function states. Its effectiveness is demonstrated through numerical experiments on the Poisson equation, interface merging in Hele–Shaw flow, and two-phase Stokes flow with surfactant.

Convergence of the bilevel policy iteration algorithm for the value function, satisfying the Hamilton-Jacobi-Isaacs (HJI) equation, is analyzed in the presence of control constraints. Both first and second-order HJI equations corresponding to the deterministic and stochastic systems are considered. For numerical tests, a semi-implicit upwind scheme for backward HJI PDEs is applied and a thorough comparison between the solutions to first and second-order Hamilton-Jacobi-Bellman (HJB) equations is presented for both the unconstrained and the constrained control cases.