Papers

We reconstruct a piecewise-smooth source in a Poisson equation on a semi-infinite strip from noisy solution values measured along an interior line. Applying the one-dimensional Dirichlet Laplacian to the observation reduces the inverse problem to regularized second-order differentiation followed by a boundedly invertible correction. The differentiated trace and the source differ by an analytic smoothing term and therefore have the same interior singular support and jump data. This structure motivates a localized Fourier extension method that combines two staggered detection partitions, GTSVD-regularized local derivative coefficients, mollified conjugate Fourier sums, and structure-aligned numerical differentiation. The source is then recovered by an exponentially decaying spectral correction. For an exact partition, the reconstruction inherits the piecewise differentiation rate $\mathcal O(\delta^{(\bar s-2)/\bar s})$; local peak and partition perturbation estimates describe the additional effect of breakpoint errors. Comparisons with full-grid total-variation regularization and truncated Fourier inversion show competitive high-noise performance and a pronounced low-noise advantage of the localized method. Repeated Gaussian-noise tests confirm robustness in the moderate- and low-noise regimes.

Gaussian elimination is one of the oldest algorithms in mathematics, and the most popular method for solving an unstructured linear system. Its stability in finite precision is controlled by its growth factor, which measures how large the entries produced during elimination can become. Understanding the worst-case behavior of this quantity has been a central problem in numerical analysis since the 1940s. Here we make a significant leap in that understanding, settling several open problems. In particular, we determine the asymptotic behavior of the maximum growth factor under complete and rook pivoting, proving that both are quasi-polynomial in dimension. We also show that the exponential growth under partial pivoting persists for sparse matrices and that randomized partial pivoting suffers the same instability. In contrast, we show that every non-singular matrix has a row permutation with polynomial growth, though finding the optimal row permutation is NP-hard.

We propose a local Legendre frame method for the accurate computation of Caputo fractional derivatives of order \(0<\alpha<1\). On each local subinterval, the function is represented by a restricted Legendre frame obtained from scaled Legendre polynomials on an extended interval. The local coefficients are computed from equispaced samples by an exponentially weighted GTSVD regularization. The Caputo derivative is then evaluated by applying the weakly singular fractional integral to the derivatives of the local frame basis functions. Since these derivatives are polynomials, the corresponding Caputo weights can be written in terms of finite weighted moments, so that the singular kernel is treated analytically rather than by a low-order quadrature rule. For uniform partitions, the history weights have a block-dependent structure and can be reused efficiently. The error analysis separates the local frame reconstruction from the Caputo integration. In particular, the Caputo error is bounded by the derivative reconstruction error, while the latter is obtained from the \(L^2\) reconstruction error and a weighted smoothness bound of the GTSVD approximation through an interpolation argument. For analytic local functions with exponential coefficient decay, this leads to exponential-type convergence of the derivative and hence of the Caputo approximation. Numerical experiments confirm the accuracy of the exact moment weights, the effectiveness of the local weighted reconstruction, and the efficiency of the block implementation.

This paper concerns the null control of quasi-linear parabolic systems where the diffusion coefficient depends on the gradient of the state variable. In our main theoretical result, with some assumptions on the regularity and growth of the diffusion coefficient and regular initial data, we prove that local null controllability holds. To this purpose, we consider the null controllability problem for the linearized system, we deduce new estimates on the control and the state and, then, we apply a Local Inversion Theorem. We also formulate an iterative algorithm of the quasi-Newton kind for the computation of a null control and an associated state. We apply this method to some numerical approximations of the problem and illustrate the results with several experiments.

We introduce adaptive resolution frames (AR-frames), a class of weighted frames equipped with a partitioned index structure that models multiple resolution levels in Hilbert spaces. We establish a block-triangular representation of the associated operators and prove the convergence of an iterative reconstruction method with explicit error estimates. We further derive optimal block-diagonal preconditioners that minimize the iteration error and obtain stability results showing that AR-frames and their canonical reconstruction operators remain stable under sufficiently small weighted perturbations. These results provide a unified framework for efficient reconstruction and perturbation analysis of AR-frames.

Convolution integrals widely exist in applications, and to enable fast and accurate computations, this paper introduces two general multi-stage neural operator learning frameworks. The first, Deep Collocation Neural Operator (DCNO), is a supervised approach that iteratively refines the operator approximation by learning residuals from input-output data pairs. The second, Deep Galerkin Neural Operator (DGNO), is an unsupervised framework applicable when the target operator can be represented by a PDE, leveraging the weak form of the PDE residual for training. Both methods progressively construct basis operators through multiple training stages to enrich the approximation space, leading to significantly improved accuracy over standard one-shot operator learning. We provide theoretical analysis for their approximation capabilities and implement them for learning convolutions. Extensive numerical experiments demonstrate that both DCNO and DGNO achieve high accuracy, approaching machine precision under single float for convolution problems, and offer substantial efficiency gains for numerous queries or parametric variations compared to traditional solvers. We also extend these frameworks to handle multi-input operator learning scenarios involving variations in both the density and kernel of a convolution.

This work establishes uniform error estimates for classical numerical schemes applied to the generalized Korteweg-de Vries (gKdV) equation in the dispersionless limit regime, specifically before the development of dispersive shock wave (DSW). We analyze the Crank-Nicolson (CN) and a Lawson-type Runge-Kutta (Lawson-RK) methods, when discretized in space via the Fourier pseudo-spectral method. We prove that both fully discrete schemes achieve optimal second-order temporal accuracy and spectral spatial accuracy, with error constants being uniform in the vanishing dispersion parameter $\varepsilon$. The analysis also addresses the unique solvability of the CN scheme in the dispersionless limit regime. Theoretical findings are supported by numerical experiments, demonstrating the $\varepsilon$-uniform accuracy and the necessity of step size restriction for CN. The study in general validates the classical methods for solving the dispersionless gKdV equation before DSW.

The classical AT1 phase-field model contains an intrinsic energy barrier for crack nucle ation, which makes the predicted strength depend on the fracture toughness and the regularization length. For tensile-dominated brittle fracture, this barrier is shifted by mapping the Rankine criterion, evaluated on the effective stress, onto a state-dependent active-energy threshold. The prescribed tensile strength then controls crack nucleation, while the AT1 crack-density functional, stiffness degradation, and degraded stress response remain unchanged. Since the threshold depends on the current stress state, the field equations are derived from a restricted variational principle. A microforce formulation identifies the barrier shift as a dissipative resistance and provides the corresponding lower bound on the regularization length. In one-dimensional tension, closed-form solutions recover the prescribed peak strength and give a cosine-type localization profile that ap proaches the classical AT1 profile as the shift vanishes. Numerical examples show that, within the admissible range, the nucleation load is nearly insensitive to the regularization length and the predicted multiaxial nucleation states follow the Rankine envelope. Under overall compression, crack nucleation remains associated with local tensile stress concen trations. The formulation also captures the transition from strength-controlled failure for small flaws to the LEFM limit for large cracks.

Combining classical group theory and partial differential equation theory, this paper investigates the symmetry group $\operatorname{Sym}(u)$ of the unique solution $u$ to the second-order linear elliptic boundary value problem on an $n$-dimensional bounded domain $\Omega$ $-\sum_{i,j=1}^{n} a_{ij}(x)u_{x_ix_j} + \sum_{i=1}^{n} b_i(x)u_{x_i} + c(x)u = f(x), x\in \Omega$, $u(x) = h(x), x\in \partial \Omega$, The following symmetry groups are defined and characterized respectively: the symmetry group $\operatorname{Sym}(A)$ of the second-order coefficient matrix function $A(x)=(a_{ij}(x))_{n\times n}$; the symmetry group $\operatorname{Sym}(b)$ of the first-order coefficient column vector function $b(x)=(b_{1}(x),b_{2}(x),\cdots,b_{n}(x))^{T}$; the symmetry group $\operatorname{Sym}(c)$ of the zero-order coefficient function $c(x)$; the symmetry group $\operatorname{Sym}(f)$ of the internal source function $f(x)$; and the symmetry group $\operatorname{Sym}(h)$ of the boundary source function $h(x)$. This paper rigorously proves that the common symmetry group $\operatorname{Sym}(A)\cap\operatorname{Sym}(b) \cap\operatorname{Sym}(c) \cap\operatorname{Sym}(f) \cap\operatorname{Sym}(h)$ is a subgroup of $\operatorname{Sym}(u)$. In addition, if the common symmetry group contains several mirror symmetry elements, the original second-order linear elliptic boundary value problem on the entire domain $\Omega$ can be reduced to the corresponding boundary value problem on a certain subdomain. It is strictly proven in this paper that the new boundary condition imposed on the boundary of the subdomain is the homogeneous generalized Neumann boundary condition. The linear finite element method is used to numerically solve the second-order linear elliptic boundary value problem on the subdomain, thereby achieving domain reduction and significantly reducing the computational cost.

We study Sobolev approximation on bounded domains by linearized shallow neural networks whose inner parameters are prescribed independently of the target function. Our main step is a one-dimensional construction for analytic activations. We prove that quasi-Chebyshev parameter sets with univariate resolution $m$ generate fixed feature spaces attaining the sharp $H^r$-to-$H^s$ approximation order $m^{-(r-s)}$ for a class of analytic activations satisfying a quantitative non-cancellation condition on their Taylor coefficients. Combining this result with the ridge-function lifting theorem in [SIAM J. Math. Anal. 30 (1998), pp. 155-189] and its extension to arbitrary quasi-uniform direction sets established in this work, we construct tensor-product-type parameter sets that attain the sharp rate $$\|f-f_n\|_{L^2(\Omega)}\lesssim n^{-\frac rd}\|f\|_{H^r(\Omega)},\quad f\in H^r(\Omega)$$ for all $r>0$. In contrast to the finite-difference construction in [Neural Comput. 8 (1996), pp. 164-177], whose explicit admissibility condition may require an extremely small parameter scale, the proposed parameter sets remain distributed over fixed intervals and are therefore more amenable to practical computation.

We propose a novel, robust, and second-order-accurate parameterized penalty-projection ensemble algorithm for incompressible Navier–Stokes flow problems. The resulting linearized algorithm, based on the second-order Backward Differentiation Formula (BDF-2), is computationally efficient because it shares the same coefficient matrix across all realizations for each subproblem at every time step. To enhance robustness in convection-dominated flows, the scheme incorporates Ensemble Eddy Viscosity (EEV) regularization. In addition, it is equipped with grad-div stabilization parameter $\gamma$, which controls the splitting error; under the assumptions of the analysis, the splitting error decreases and vanishes asymptotically as $\gamma\to\infty$. We establish the stability of the proposed scheme and rigorously prove its optimal convergence by demonstrating that, as $\gamma\to\infty$, the scheme converges to an equivalent coupled formulation. We further validate the method through a series of numerical experiments designed to verify the theoretically predicted convergence rates and evaluate its performance on benchmark convection-dominated problems. The numerical results are in excellent agreement with the theoretical analysis and confirm the effectiveness of the proposed scheme.

This paper develops a computational framework to optimize the process parameters such that the bond quality between extruded polymer filaments is maximized in fused filament fabrication (FFF). A transient heat transfer analysis providing an estimate of the temperature profile of the filaments is coupled with a sintering neck growth model to assess the bond quality that occurs at the interfaces between adjacent filaments. Predicting the variability in the FFF process is essential for achieving proactive quality control of the manufactured part; however, the models used to predict the variability are affected by assumptions and approximations. This paper systematically quantifies the uncertainty in the bond quality model prediction due to various sources of uncertainty, both aleatory and epistemic, and includes the uncertainty and the model discrepancy in the process parameter optimization. Variance-based sensitivity analysis based on Sobol indices is used to quantify the relative contributions of the different uncertainty sources to the uncertainty in the bond quality. A Gaussian process (GP) surrogate model is constructed to compute and include the model discrepancy within the optimization. Physical experiments are conducted for calibration and validation of the physics model and also for validation of the optimum solution. The results show that the proposed formulation for process parameter optimization under uncertainty results in high bond quality between adjoining filaments of the FFF product.

This work presents a data-driven methodology for multi-objective optimization under uncertainty of process parameters in the fused filament fabrication (FFF) process. The proposed approach optimizes the process parameters with the objectives of minimizing the geometric inaccuracy and maximizing the filament bond quality of the manufactured part. First, experiments are conducted to collect data pertaining to the part quality. Then, Bayesian neural network (BNN) models are constructed to predict the geometric inaccuracy and bond quality as functions of the process parameters. The BNN model captures the model uncertainty caused by the lack of knowledge about model parameters (neuron weights) and the input variability due to the intrinsic randomness in the input parameters. Using the stochastic predictions from these models, different robustness-based design optimization formulations are investigated, wherein process parameters such as nozzle temperature, nozzle speed, and layer thickness are optimized under uncertainty for different multi-objective scenarios. Epistemic uncertainty in the prediction model and the aleatory uncertainty in the input is considered in the optimization. Finally, Pareto surfaces are constructed to estimate the tradeoffs between the objectives. Both the BNN models and the effectiveness of the proposed optimization methodology are validated using the actual manufacturing of the parts.

We prove an impossibility result for finite-range locally computable a posteriori error estimators for the conforming finite element discretization of an elliptic interface problem. For a class of interface problems with a checkerboard cross-point and $\{1,M\}$-valued coefficients, we construct two problem instances on the same interface-fitted mesh. The two instances share a piecewise-constant load for which the finite element solution and the data oscillation both vanish. The ratio of their exact energy errors grows at least proportionally to $M^{1/4}$. Locality together with efficiency forces identical estimator values for these instances. The product of reliability and efficiency constants is therefore bounded below by a constant multiple of $M^{1/4}$. Consequently, any locally computable error estimator for interface problems, whether of residual, equilibrated, or recovery type, cannot be simultaneously reliable and efficient with contrast-independent constants.

We introduce a bilingual domain-specific language (DSL) for modelling compute kernels within a generic solver for hyperbolic partial differential equations (PDEs). Users express PDE terms, i.e.~the underlying physics, in a familiar native language such as C or C++, while the numerical scheme is specified in a Python-embedded DSL, DSLHyPE. DSLHyPE's compiler lowers the Python description to MLIR and introduces a translation pass that integrates it with native code likewise mapped to MLIR. Our approach keeps the numerical representation and the physics implementation separate for as long as possible, while delegating optimization to the compiler through existing MLIR optimization passes. This separation of concerns benefits researchers developing numerical schemes on top of existing PDE implementations or with applications involving nonlinear systems whose PDE terms must solve PDEs themselves. We demonstrate the feasibility of the approach using a gravitational-wave solver and a matter-evolution solver on x86 processors and H200 GPUs.

We consider plane-wave discontinuous Galerkin (PWDG) approximations of the Helmholtz equation with adaptive local propagation directions. The directions are chosen by minimizing a weighted residual on the mesh skeleton. We study two formulations: in Part A the PWDG system is solved for the coefficients at fixed directions, while in Part B the same residual is minimized jointly over coefficients and directions, with the coefficients eliminated by variable projection. Complex angles unify propagating and evanescent Trefftz waves. The discrete system is normalized in the Trefftz-DG norm, and a local Cauchy-trace Gramian is used to remove numerically dependent directions. On straight edges the residual integrals are exact. For all-Dirichlet problems, the residual equals the squared DG error and provides local adaptive indicators. We prove local quadratic growth of the reduced direction functional near an identifiable zero-residual solution of fixed rank. Numerically, an exact circular DtN test shows that the dominant phase direction of a Hankel wave can be recovered accurately even when a small plane-wave fan gives a larger field error. We also separate the effects of the coefficient basis, trace cutoff, and trace-spectrum arithmetic on the high-$p$ DtN floor. Finally, for finite sums of plane waves, residual-based ENRICH–MOVE continuation recovers all directions to roundoff for $M=1,\ldots,19$. At $M=20$ the automatic birth step enters a false basin, whereas a nearby birth again reaches roundoff. The results indicate that direction adaptation is most effective for solutions of low directional complexity.

The artificial compressibility method (ACM) relaxes the incompressibility constraint in the Navier-Stokes equations by introducing a perturbation term proportional to the time derivative of the pressure, scaled by a small positive parameter $\varepsilon$. Consequently, the ACM system inherently involves two small parameters: the fluid viscosity $\nu$ and the artificial compressibility parameter $\varepsilon$. While existing temporal analyses of ACM provide insights into its behavior, rigorous fully discrete error estimates that are robust with respect to both parameters remain a significant gap in the literature. In this paper, we propose a second-order Crank-Nicolson fully discrete ACM scheme and establish its parameter-robust optimal error estimates. To enhance stability and ensure robustness, we incorporate two distinct grad-div stabilization strategies: one facilitates the decoupling of velocity and pressure computations, while the other guarantees robustness with respect to both $\nu$ and $\varepsilon$. For spatial discretization, we employ the Scott-Vogelius finite element pair, which is crucial for the decoupling and error analysis. The resulting parameter-uniform bounds are crucial for ensuring the long-time accuracy of the scheme, circumventing the exponential dependence on the Reynolds number typically introduced by Grönwall's lemma. Numerical experiments are provided to validate the theoretical findings and demonstrate the efficiency of the proposed methods.

Bayesian targets may converge under model refinement even when the exact sensitivities used by gradient-based samplers do not. We study this probability–sensitivity mismatch and its consequences for Metropolized Hamiltonian proposals. A vanishing-amplitude wiggly-energy model first gives the basic analytic obstruction: the potential perturbation tends to zero while its classical derivative is of order $r_\varepsilon/\varepsilon$. We then show that the same scaling arises naturally in a periodic elliptic inverse problem, where homogenization makes the forward map and Gaussian likelihood converge while differentiation with respect to a microscopic scale parameter retains an $O(1)$ oscillatory contribution. This provides a PDE origin for the single-scale wiggly mechanism. The main construction concerns a more demanding nested Weierstrass hierarchy, interpreted as an analytically tractable prototype for repeated corrector contributions across geometrically separated scales. There all previously resolved scales persist, adjacent classical-force increments grow geometrically like $(ab)^N$, and the limiting rough component may fail to possess a classical derivative. In this self-similar setting the matched Jackson quotient is structurally adapted to the refinement through dilation covariance and exact finite closure. Uniform negative-log-likelihood approximation yields explicit total-variation, Hellinger, and bounded quantity-of-interest bounds. Measurable kick–drift–kick maps remain exact after Metropolis correction, local field convergence propagates to fixed-length proposals and kernels, and the first classical HMC half-kick can have no fixed-step refinement limit. Numerical experiments on scale-structured inverse problems test the resulting resolution-stability mechanism across one- and two-dimensional inverse problems.

The inverse Born series provides a general perturbative framework for representing nonlinear inverse maps between infinite-dimensional function spaces and has found numerous applications in inverse problems governed by partial differential equations and integral equations. Motivated by its operator-theoretic structure, we develop a systematic framework for constructing neural operators that approximate the operator expansions arising in the inverse Born series. Our approach combines the analytical representation of the inverse Born series with the expressive power of neural operators, yielding data-driven approximations of the nonlinear inverse map while preserving the underlying operator structure. The proposed framework is applicable to a broad class of inverse problems and is presented in a general functional-analytic setting. To demonstrate its practical performance, we consider two representative examples: inverse scattering and the Calderón (electrical impedance tomography) problem. Numerical experiments show that the constructed neural operators accurately approximate the inverse Born expansions and produce high-quality reconstructions across a range of test cases. These results indicate that the proposed methodology provides an effective and computationally efficient approach for learning nonlinear inverse operators and suggests a promising direction for integrating classical operator expansions with modern neural operator architectures.

We introduce Ordered Diffusion Kernels (ODKs), a novel class of local kernels that can approximate the infinitesimal generator of an arbitrary Itô Stochastic Differential Equation (SDE). ODKs are designed to be applied to data sampled from dynamical systems where little dynamical information is available a priori. The Laplacian of classical diffusion kernels approximates the Laplace-Beltrami operator on the underlying manifold; adjusting the normalisation introduces an advection term that depends on the sampling density; recently, TMDmap generalised this normalisation to target an arbitrary measure, but at the cost of coupling advection to diffusion. More general local kernels can learn arbitrary second-order elliptic operators but are formulated in terms of known velocity fields — making the first step of any analysis a potentially ill-posed inference problem. To formulate ODK, we first relax the problem of potential estimation to the more tractable task of inferring an ordering of the data, which we represent through an ordering function. We prove ODK's Laplacian converges to the infinitesimal generator of a gradient-flow SDE with state-dependent isotropic diffusion, without coupling advection and diffusion. We provide various extensions of ODK to: arbitrary drifts via local ordering functions; anisotropic diffusions via a Strang splitting scheme; multiple ordering functions; and self-tuning bandwidths. In addition, we introduce two loss functions which exploit the structure of ODKs to solve a non-parametric inference problem. We validate this framework on synthetic data from deterministic and stochastic systems, demonstrating accurate recovery of operators, velocity fields, extrinsic curvature, and spatially dependent drift and diffusion coefficients.

Transport-dominated problems remain challenging for data-driven methods, which often exhibit severe numerical oscillations near shocks or steep gradients due to globally supported basis functions or overly smooth hypothesis spaces. Gegenbauer reconstruction has shown promise in mitigating such oscillations, but its effectiveness critically depends on the reconstruction parameters, particularly the weight parameter $\lambda$ and truncation order $m$. For data-driven models, variations in governing problems, training data, and model architectures make systematic parameter selection particularly challenging. To address this issue, we propose a physics-informed machine-learning framework that predicts probability distributions over candidate Gegenbauer parameter pairs, enabling probabilistically weighted reconstruction while accounting for parameter uncertainty. A two-stage strategy is adopted, in which a general predictor is first pre-trained and then fine-tuned for target problems to balance accuracy and computational cost. The framework is evaluated for reduced-order and neural operator models, represented by POD-Galerkin and DeepONet, respectively. Numerical experiments on one- and two-dimensional transport-dominated problems show that the framework learns effective spatially adaptive parameter distributions. Compared with conventional reconstruction strategies, it reduces numerical errors by up to one to two orders of magnitude and achieves a more favorable accuracy–cost trade-off than problem-specific model retraining.

The complex-step method is traditionally derived from the Taylor expansion of an analytic function and is widely used as a numerical technique for derivative approximation. We present an alternative formulation based on the Cauchy–Riemann equations and show that the classical complex-step relation arises naturally from the harmonic structure of holomorphic functions. In particular, the complex-step method admits two complementary harmonic interpretations: as a Cauchy problem, in which the derivative is identified with the normal datum of the imaginary component on the real axis, and as a reconstruction problem in a strip, in which the finite imaginary perturbation provides the upper-boundary data. The latter formulation leads explicitly to the strip Poisson and conjugate Poisson kernels and their derivatives. A related harmonic reconstruction framework in the upper half-plane leads to the Poisson, conjugate Poisson, and Cauchy kernels as elementary reconstruction operators for harmonic and holomorphic functions. Extending this reconstruction from ordinary boundary functions to finite measures yields the classical Stieltjes transform and its inversion formula. The same measure-theoretic structure appears in spectral theory, where scalar matrix elements of the resolvent are Stieltjes transforms of the associated spectral measures. These results establish a common complex-analytic structure connecting complex-step differentiation, harmonic reconstruction, Stieltjes inversion, and spectral reconstruction, while distinguishing the boundary-value problems through which the corresponding information is recovered.

The numerical simulation of time-harmonic wave propagation in elastic media plays an important role in applications such as geophysics and non-destructive testing. Accurate discretization of the elastic Helmholtz equation at high frequencies is challenging due to numerical dispersion and pollution effects. In this work, we develop an asymptotic dispersion correction for a Marker-And-Cell (MAC) discretization of the isotropic elastic Helmholtz equation. We characterize the discrete dispersion relation of the scheme and determine the leading-order term in the dispersion error in both two and three spatial dimensions. Based on this analysis, we derive a correction that is asymptotically optimal in the limit of vanishing mesh size. The proposed approach improves the agreement between the discrete and continuous wave propagation properties while preserving the structure of the underlying discretization. We also establish a connection between the factorization of the dispersion relation and the structure of the grad-div operator symbol, providing additional insight into the algebraic structure of the elastic problem. Numerical experiments finally demonstrate a substantial reduction of relative errors and confirm the effectiveness of the proposed correction. We further provide numerical evidence that the corrected discretization improves the convergence behavior of multigrid solvers.

We establish a discrete-to-continuum variational theory for area-preserving parameterization of surfaces based on the stretch energy. For orientation-preserving diffeomorphisms between compact, connected, oriented Riemannian surfaces of equal total area, we prove that area-preserving maps are precisely the critical points of the stretch energy under boundary-fixing variations. We then establish discrete-to-continuum convergence of the optimal energy values: within a uniformly geometrically controlled admissible class, the discrete infimum converges to the smooth minimum with second-order accuracy, while discrete almost minimizers have first-order decay of their $L^2$ area distortion. Thus, the discrete problem not only approximates the smooth stretch energy but also recovers area preservation in the refinement limit. We further derive the $L^2$-gradient flow of the stretch energy and its simplicial counterpart, whose projected quasi-implicit iteration is shown to be globally convergent. The resulting method applies to open and closed surfaces of several topological types, and numerical experiments on benchmark surfaces illustrate its improvement over existing methods.

Estimation of eddy-current losses in conducting non-laminated components of electrical devices requires expensive three-dimensional simulations. Various approximations are therefore used in practice to reduce the computational cost in the early design phase. We review some approaches and discuss their modelling assumptions and resulting approximations. In particular, we identify eddy-current reaction fields as a significant contribution that should be accounted for globally. These reaction fields can be approximately reconstructed from two-dimensional magnetostatic simulations by solving a single linearized time-periodic problem. We further discuss different strategies for solving this post-processing problem. Numerical results demonstrate improved loss prediction compared to standard post-processing at moderate additional cost.

We derive duality-based $\textit{a posteriori}$ error identities for a broad class of subgradient flows induced by time-dependent convex integral functionals. Starting from the Brézis-Ekeland-Nayroles principle, we identify an unsteady primal energy functional and derive its Fenchel dual formulation, including strong duality and the corresponding optimality system under general normal-integrand assumptions. This Fenchel duality framework is used to derive $\textit{a posteriori}$ error identities for subgradient flows. In doing so, we depart from the usual duality-based $\textit{a posteriori}$ error control framework in the unsteady setting, since the Brézis-Ekeland-Nayroles formulation reveals the following unsteady feature: the minimal primal value and the maximal dual value are both prescribed by the initial datum. This allows us to pass from a combined primal-dual gap identity to separate primal and dual gap identities. These identities quantify the primal and dual errors independently and admit representations in terms of generalized Bregman divergences and, under a spatial convex conjugation formula, as non-negative time-space integral quantities suitable for localization. The abstract framework is applied to a number of variational problems of physical interest, including the unsteady heat equation, the unsteady Stokes equations, the unsteady Navier-Lamé equations, the unsteady Bingham flow through a pipe, the unsteady obstacle problem, and the unsteady elasto-plastic torsion problem.

This paper presents an improved full-wave solver for finite array antennas, and uses this solver to examine asymptotic properties of active impedance for large regular arrays. The improved solver is based on a preconditioning scheme that has been adapted for use on more general finite geometries, and an improved data structure. These two improvements result in a fast and stable solver with a lower memory footprint. By investigating the active impedance of finite arrays it is found that, even for arrays with 1000 elements, asymptotic behavior may differ from infinite arrays. Tied to these properties, different predictors on how the active impedance of the center element in the array behaves are presented. The two best predictors work very well for the wideband arrays investigated, and may be used to decide when array approximations, such as unit cell methods, are appropriate, rather than general statements on array size.

In this paper, we investigate the convergence of an implicit defect-correction solver for a viscous discretization based on the alpha-damping scheme for unstructured grids. We show that significantly more robust iterative convergence is achieved by evaluating the damping term at the midpoint between two adjacent cell centers (edge midpoint) rather than at the face centroid. A one-dimensional Fourier analysis reveals that the implicit solver tends to be stable when the damping term in the residual is smaller than that used to construct the Jacobian. This observation suggests that the solver can be stabilized by effectively reducing the magnitude of the damping term in the residual - an effect achieved by the edge-midpoint evaluation. Robust convergence is demonstrated numerically for two-dimensional viscous-flow problems on highly irregular mixed-element and triangular grids.

We develop a quantum algorithm for slow analytic Hamiltonians $\widetilde H(t)=H(t/T)$ with $\|H(s)\|\leq\alpha$ that achieves nearly additive query complexity and low gate overhead. Our main technical contribution is a periodic Gevrey extension of $H(s)$, together with Fourier component decay and truncation bounds that enable an efficient finite-dimensional simulation. Combined with Floquet embedding and optimal time-independent Hamiltonian simulation technique, this gives query complexity $\widetilde{\mathcal O}\!\left(\alpha T+\log(1/\varepsilon)\right)$ and additional gate complexity $\widetilde{\mathcal O}\!\left((\alpha T+\log(1/\varepsilon))^2\log(1/\varepsilon)\right)$, assuming coherent access to $H'(s)$ and endpoint derivatives. For slow analytic control Hamiltonians, only block encodings of the time-independent control operators are required, with the same query complexity and lower gate overhead. Our method also extends to Gevrey Hamiltonians and improves the precision dependence for simulating slow analytic semi-dissipative linear differential equations.

This paper considers European multi-asset option pricing under Lévy and affine characteristic-function models. The main obstruction is the curse of dimensionality: direct multidimensional COS pricing forms tensor-product coefficient arrays whose size grows exponentially with the number of assets. We study and extend COS-TT-CHF, a low-rank construction that uses TT-cross to compress sampled characteristic-function tensors into tensor-train COS coefficients for arithmetic basket and min/max option pricing. Once built, the compressed representation gives fast post-setup strike-grid and selected component Delta/Vega calculations. The numerical study compares with adaptive-quadrature Fourier benchmarks, direct COS, a tensor-Fourier min-option benchmark, and quasi-Monte Carlo (QMC) references based on randomized Sobol points. The reported timings show a low-dimensional crossover against direct COS as the benchmark moves from $d=2$ to $d=4$, favorable timings against the tensor-Fourier min-option benchmark from $d=3$ onward, and favorable timings against the QMC common-Heston reference already at $d=2$. The reported tests reach $d=30$ for GBM and $d=20$ for VG, NIG, and common-Heston benchmark families, with accuracy, rank, runtime, control-sensitivity, and component Delta/Vega diagnostics reported throughout.