Papers

Near-field holotomography combines coherent diffraction imaging with tomographic acquisition to recover the three-dimensional complex refractive index of a specimen. Since the measured diffraction intensities are generated by a nonlinear object transmission and wave propagation process, the resulting inverse problem is intrinsically nonlinear and ill-posed. We study a direct variational reconstruction framework based on a fully nonlinear wave propagation model that avoids both intermediate phase retrieval and linearization under the weak-object approximation. We analyze the forward operator in appropriate Banach spaces, establish its Fréchet differentiability, and derive explicit gradient expressions for variational reconstruction. To improve quantitative reconstruction of weak absorption features, we further develop a phase-guided Bregman TV regularization framework that exploits structural correlations between phase and absorption components. This enables multi-material reconstruction without imposing a globally fixed ratio between the two components. We perform numerical studies on synthetic phantoms and experimental data. The results demonstrate stable three-dimensional reconstructions and improved recovery of the absorption contrast compared to state-of-the-art methods.

This article presents a fast, high-order Nyström solver for the two-dimensional Lippmann–Schwinger equation arising from time-harmonic scattering by penetrable, piecewise-smooth heterogeneous media. Relying on high-order evaluation of the Newtonian potential on unstructured grids adapted to interfaces of discontinuity, the methodology achieves high-order accuracy using existing fast algorithms such as the fast multipole method. As an iterative method the solver exhibits rapid convergence when coupled to a preconditioning strategy that exploits a class of structured-grid solvers—fast solution methods offering quasi-linear time and memory complexity and nearly-constant iteration counts, but long limited in accuracy. The preconditioning strategy couples the Nyström discretization—given by high-order quadrature nodes over a (curved) unstructured mesh conforming to the support of the spatially varying contrast—to a uniform Cartesian grid underlying the fast preconditioner via a pair of transfer operators. The resulting preconditioner inherits the frequency-robust behavior of its Cartesian counterpart without sacrificing the geometric flexibility and high-order accuracy of the unstructured discretization. We prove that invertibility of the proposed preconditioner holds under explicit conditions on the mesh sizes and on the Cartesian preconditioner. Numerical experiments demonstrate that the preconditioned system requires significantly fewer GMRES iterations than its unpreconditioned counterpart, with iteration counts almost independent of mesh size and wavenumber, and illustrate the method's robustness for inhomogeneities with piecewise-smooth refractive indices and jump discontinuities across interfaces.

For a mixed finite element method, a Dirac source is first a failure of duality, not of regularity: the conservation equation is tested against a Lebesgue space, and a Dirac measure lies in the dual of none. We therefore remove the measure from the conservation law by a divergence-form splitting. An explicit field whose divergence is the Dirac measure is subtracted from the physical flux, and the modified flux is taken as the mixed unknown, so that only the regular part of the load remains in the conservation equation. Equivalently, and independently of any discretization, the Dirac problem is rewritten as an elliptic equation whose data are in divergence form, generated by a field of L^p. The subtracted field depends on the location of the source alone and not on the coefficient. No coefficient-dependent singular solution and no discrete delta is needed, only the load vector of the RT_0-P_0 system changes, and the source may sit anywhere relative to the mesh: at a vertex, on a coefficient interface, or inside an element. Unless the splitting is matched to the operator at the source, the modified flux lies in L^p for every p<2 but not in L^2, so the flux error analysis has to leave the Hilbert scale. We prove a quasi-best approximation bound for the flux, and with it that on a quasi-uniform family the flux error is exactly of order h^(2/p-1): the matching lower bound comes already from the single element carrying the source. Grading the mesh there restores first-order complexity, N^(-1/2) in the number of elements, and the adaptive computations attain it. The scalar variable is limited only by piecewise constant approximation of the solution, which it attains. We also prove a residual norm equivalence in the Lebesgue scale, yielding a computable L^p estimator, reliable and locally efficient for the mixed flux together with a recovered potential.

An implicit Euler finite-volume scheme for degenerate pseudo-parabolic cross-diffusion equations is proposed and analyzed. The system describes the dynamics of an unsaturated two-phase flow mixture with dynamic capillary pressure in a porous medium. The numerical scheme is based on a two-point flux approximation that preserves the energy structure, ensures the conservation of total mass, and guarantees strict positivity and boundedness of the water saturation. These properties rely on carefully selected mean functions for the nonlinear components. The existence and uniqueness of a discrete solution and additional mesh-uniform bounds are proved. Numerical experiments in two space dimensions illustrate the effects of the dynamic capillary pressure.

Inverse learning within a statistical framework has a wide range of applications. It has garnered significant attention in machine learning, artificial intelligence, and related fields, where the goal is to infer unknown parameters from indirect and noisy observations. This work investigates the stable approximation of $u^{\dagger}$ which solves the equation $Au=g$, with $A$ being a linear operator between appropriate vector spaces. We will consider the domain to be a non-reflexive Banach Space and the co-domain to be a space of real-valued functions on a metric space $X$. The function $g$ is characterized by a finite number of independently and identically distributed data points, which are assumed to follow some unknown probability measure $\rho$. We employ Tikhonov regularization with an arbitrary convex functional to obtain the regularized solution corresponding to the given data point. The convergence analysis is carried out with respect to the Bregman distance, and an upper bound for the error is derived in probability terms. The theoretical findings are then supported by numerical experiments.

Moving least-squares approximation is a popular method for approximating multivariate functions from given discrete data. For higher accuracy higher degree polynomials have to be used, resulting also in higher computational cost and numerical instabilities. Recently, the combination of low-order moving least squares with a multilevel scheme showed superior numerical behavior. In this paper we will prove, amongst other things, that such a combination of moving least-squares with a multilevel scheme indeed leads to improved convergence results, at least if the data sites form a regular grid.

This work develops tensor-train (TT) formulations for solving large-scale three-dimensional linear elasticity problems discretized by isogeometric analysis. By exploiting the tensor-product structure of the basis functions and the low-rank structure of geometry-dependent coefficient fields, the stiffness operator, mass operator, force vector, and displacement solution are represented in TT format. Two solution strategies are investigated: a block-operator formulation, in which the coupled elasticity operator is stored as separated TT blocks, and a single-operator formulation, in which the full coupled system is stored as one monolithic TT operator. A matrix-free three-field TT conjugate-gradient solver is introduced for the block formulation, while AMEn is used for the single-operator formulation. Numerical examples demonstrate substantial compression of both operators and solutions compared with conventional sparse full-grid representations, showing that TT-based formulations provide an efficient and scalable approach for large-scale three-dimensional elasticity simulations.

Contact Hamiltonian dynamics gives dissipative mechanics an intrinsic action variable, but explicit contact splittings reach only kinetic energies whose terms are exactly integrable: frozen-coordinate diagonal metrics (the spherical pendulum, a torus particle) are included, while dense metrics with momentum cross terms, with the double pendulum as flagship, are not. We introduce a projected Pihajoki-contact integrator for this non-separable setting, combining phase-space duplication, symmetric projection onto the physical diagonal, and constant-friction damping half-steps, with the action factor carried by an exact Herglotz update. As in the projected extended-phase-space framework it builds on, the construction needs no binding parameter, returns the copies to the diagonal at every step, and confines the nonlinear solve to the $2n$ projection variables. For constant friction the step rescales $\omega=d\eta$ by the exact factor $e^{-\gamma\tau}$ when the projection is solved exactly (a classical conformally symplectic identity, realized here for this class), while time-symmetry, consistency, and smoothness yield an $O(\tau^3)$ one-step contact-form residual, a bound not specific to the contact form. On the damped double pendulum, spherical pendulum, and torus particle the method is second-order accurate, reproduces the contact decay law, and controls long-time energy and contact drift in coarse or stiff regimes where the Tao baseline and the unprojected average lose the solution. A head-to-head with exact-contactomorphism splittings delimits the niche: where a frozen-coordinate splitting exists it preserves the contact form exactly and wins at matched cost; for the dense double-pendulum metric the realizable alternative is first-order with a prohibitive constant and the projected method prevails. The contact-form estimate is local, one-step, and constant-friction.

We develop a general framework for constructing robust high-order projector-splitting integrators for dynamical low-rank approximation. For a prescribed matrix increment, we show that the standard K-S-L projector-splitting step is equivalent to a reduced K-L step, thereby eliminating the explicit backward S-step. We then establish a central relaxed exactness property of the standard projector-splitting integrator: the rank-$r$ approximation inherits the accuracy of the numerical matrix increment, with an error bound independent of small singular values. The resulting schemes evolve fixed rank-$r$ factors and require neither basis augmentation nor rank truncation. As concrete examples, we combine the framework with selected second- and third-order Runge–Kutta methods to obtain robust high-order projector-splitting integrators. Numerical experiments confirm the predicted uniform convergence rates with respect to small singular values.

Quasi-Monte Carlo sampling is a numerical integration method that uses points with a space-filling property in $[0,1]^s$ designed to give better estimates than plain Monte Carlo methods do. For integrands of bounded variation in the sense of Hardy and Krause, errors of $O(n^{-1+\epsilon})$ for any $\epsilon>0$ are obtained from $n$ sample points. Randomized quasi-Monte Carlo (RQMC) points are individually uniformly distributed but collectively space-filling and then independent replications provide variance estimates. For smooth enough integrands the randomization can give a root mean squared error of $O(n^{-3/2+\epsilon})$. This article explains RQMC for a statistical readership recounting some history and presenting some current directions.

Solving large-scale $\ell_1$-regularized image deblurring problems efficiently while preserving sharp edges remains a significant computational challenge. We propose a wavelet-based multilevel framework that embeds three iterative solvers, Iteratively Reweighted Least Squares (IRLS), Split Bregman (SB), and Majorization-Minimization (MM), within a multilevel V-cycle. Discrete wavelet transforms define the interlevel transfer operators, and regularization parameters are selected automatically by Generalized Cross Validation. Two information transfer strategies are introduced and compared: one transfers only the coarse solution to the fine level, while the other transfers solver-specific auxiliary quantities. Numerical experiments demonstrate substantial computational savings for IRLS, with speedups exceeding an order of magnitude, while MM and SB exhibit more modest computational differences. The experiments generally show that transferring auxiliary iterates performs best with Haar wavelets, whereas transferring only the solution performs best with Daubechies wavelets.

In this work, we propose a new Arrow-Hurwicz iterative scheme designed to solve the steady inductionless magnetohydrodynamics system. The main feature of the proposed scheme is the introduction of new penalty terms in the current density equation. These terms play a central role in effectively controlling the unfavorable mixed terms that commonly arise in such formulations. The proposed method is shown to achieve geometric convergence. Numerical tests affirm the efficiency of the new scheme without compromising accuracy.

We analyze lattice-based algorithms for recovering integer-valued signals from partial discrete Fourier transform (DFT) measurements. These algorithms formulate signal recovery as the problem of finding short vectors in an appropriately constructed lattice. We derive parameter estimates that predict the regimes required for successful recovery and quantify how these estimates depend on the signal length, the error of an initial guess, and the number of sampled DFT coefficients. The analysis also reveals a close connection between the lattice parameters and the precision of the measured DFT coefficients, providing insight into the numerical stability of the inversion. Numerical experiments demonstrate close agreement between the theoretical predictions and observed recovery thresholds over a broad range of problem parameters.

Snow is a complex geomaterial whose macroscopic response is governed by density, temperature, and the topology of its evolving microstructure. Its mechanical behavior spans elastic, plastic, viscous, and failure dominated regimes, imposing significant challenges for numerical methods, which intends to simulate large deformations, evolving free surfaces, and complex boundary interactions. This work presents the first integration of a Modified Cam-Clay constitutive formulation for snow into a purely meshfree strong-form collocation framework based on the Generalized Finite Difference Method. The main methodological contribution is a numerical coupling that combines a global implicit mixed formulation for pressure and velocity with a constitutive return-mapping algorithm. The hydrostatic pressure contribution is obtained from a Poisson equation and subsequently corrected through the Modified Cam-Clay return-mapping procedure, while the deviatoric response is treated semi-implicitly using a numerical viscosity formulation. This partitioned treatment of the volumetric and deviatoric stress contributions enables stable simulations with comparatively large time steps while producing smooth spatial pressure fields. As a result, forces on complex boundary geometries can be evaluated accurately. Numerical results of this coupling illustrate the algorithmic stability of the framework, the effective imposition of boundary conditions, and the suitability of local spatial refinement. The feasibility of applying the framework to vehicle-snow interaction through rigid-body coupling is also illustrated. The presented formulation provides a robust basis for future simulations of dynamic snow loading on vehicle structures.

We introduce Automatic Symbolic Regression (AutoSR), a fully automated system that instantiates Research-Space Symbolic Regression by searching persistent scientific investigations rather than isolated equations. Finite, noisy data often yield numerically competitive expressions that imply very different behavior outside the observed regime, making numerical fit and syntactic complexity insufficient measures of scientific credibility. Existing approaches largely focus on improving expressions, yet the search typically retains little beyond the resulting formula and score, losing the scientific record, such as motivations and probes, that inform what to try next. AutoSR preserves this record in a Research State, coupling each candidate equation with the reasoning, computational evidence, and independent review developed along its branch. Proposer–reviewer agents develop these states under progressive-widening Monte Carlo tree search (PW-MCTS), which allocates computation across competing investigations, while the accumulated research record is ultimately synthesized into a final report that explains the leading relation and the basis for its selection. Across nine selected challenges from two benchmark suites, AutoSR recovers algebraically equivalent relations in every case, including three cp3-bench problems that no published system recovers and six structurally diverse LSR-Transform problems. Overall, AutoSR extends symbolic regression from equation-level search toward automated scientific investigation, allowing scientific knowledge and accumulated evidence to shape both what is explored and how the resulting equation is justified.

We consider the Euclidean distance geometry problem (EDG): given a subset of the pairwise distances of an unknown cloud of $n$ points in $\mathbb{R}^\ell$, recover the point cloud up to rigid motions. When $n$ is large, a popular practical approach is to minimize a nonconvex quartic, known as the squared-stress or s-stress, over point clouds in $\mathbb{R}^k$, with $k$ potentially larger than $\ell$. It is a long-standing open problem to understand the optimization landscape of the s-stress when all pairwise distances are known (Malone and Trosset, 2000; Parhizkar, 2013). It was recently shown that the landscape is not benign when $k=\ell$, and it was conjectured that the landscape becomes benign as soon as $k\ge \ell+1$ (Song et al., 2025; Criscitiello et al., 2026). Here, we show that the complete-graph s-stress has a benign landscape whenever $k\ge 2(\ell+1)$, establishing the conjecture up to a factor of two. A key idea is to view second-order criticality as a containment of two ellipsoids; finding a descent direction then corresponds to finding a separating hyperplane that violates this containment. This dual perspective yields the stated landscape result, and also applies to any measurement operator whose inverse satisfies a simple frame condition.

We study the nonlinear stochastic Cahn-Hilliard equation per- turbed by additive colored noise. We show almost sure existence and regularity of solutions. We introduce spatial approximation by a standard finite element method and prove error estimates of optimal order on sets of probability arbitrarily close to 1. We also prove strong convergence without known rate.

In Vilenkin's group we present an algorithm for constructing a step scaling function with a given support and a constant on given cosets.

We propose and analyze a time-staggered numerical scheme for mean field game (MFG) systems with individual noise on finite graphs. Numerically solving such coupled forward–backward systems is delicate because the density evolves in the open probability simplex and the coefficients may degenerate at its boundary. The scheme preserves mass and satisfies a discrete fundamental identity compatible with the Lasry–Lions monotonicity argument, leading to uniqueness of the numerical solution. By establishing a timestep-uniform positive lower bound for the density and uniform bounds for the value variable, we prove first-order convergence for every interior discrete solution. For potential MFGs, we establish a variational characterization by identifying the scheme with the KKT system of a convex discrete action, yielding existence of the discrete solution and an optimization-based realization. The resulting optimization problem is solved by a feasible primal–dual Newton method in mass-preserving coordinates. Numerical experiments confirm the predicted convergence rate and illustrate topology-dependent transport and congestion-driven route choice.

Neural PDE training yields a finite checkpoint archive, yet its logged energy errors are inaccessible without the exact solution, while loss-based selection does not necessarily recover the logged energy oracle. For admissible neural approximations of symmetric coercive variational problems, we introduce a reference-free selection rule based on minimizing a computable conforming Riesz monitor. The exact residual-energy identity and conforming projection make the monitor an unconditional lower bound converging monotonically to each logged energy error under nested conforming refinement; under saturation, hierarchical enrichment yields a computable upper estimate and hence a lower-upper bracket. A key finding is that archive selection is order-sensitive: unresolved checkpoint-dependent components can reverse the oracle-non-oracle ranking at finite resolution, so checkpointwise recovery alone is insufficient. For finite archives, we prove uniform recovery, yielding convergence to the logged-oracle error and, without saturation, logged-oracle selection at sufficiently fine auxiliary resolution. Under saturation, the bracket gives a computable near-oracle bound and certifies unique logged-oracle selection upon interval separation. We also bound logging-resolution loss and certify oracle inclusion over prescribed comparison trajectories. The resulting criterion replaces inaccessible exact-error minimization by computable, training-independent post-training selection on the intrinsic energy-error scale, requiring only the computed candidates and the variational problem. Experiments on diffusion and elasticity, including a non-manufactured perforated plate, demonstrate energy-scale calibration, oracle-level selection, and modest post-processing cost.

Elastic wave propagation in fractured media is relevant to applications such as analysis of seismic waves and non-destructive characterization of materials. Understanding attenuation and scattering behavior arising from wave-fracture interaction is important for interpreting observations at both field and laboratory scales. This work presents a computational framework for elastic wave propagation in fractured media based on a mixed-dimensional discrete fracture-matrix representation. Fracture deformation is governed by four models of increasing complexity, ranging from widely used spring-based formulations to fracture contact mechanics models with friction, all incorporated within a unified computational framework. Many previous studies are often restricted to simplified wave fields, single fractures or subsets of the relevant fracture deformation mechanisms. In contrast, the proposed framework enables fully coupled simulation of elastic wave propagation with fracture deformation models that account for elastic normal deformation, frictional contact and fracture opening and closure. The elastic wave equation is discretized in space using the cell-centered finite volume method Multi-Point Stress Approximation with weak symmetry and in time using the Newmark method. The spatial discretization is locally conservative and applicable to general polyhedral grids, making it well suited for media containing fractures, material heterogeneities and anisotropy. The proposed framework is verified through numerical convergence analyses and is subsequently applied to wave propagation and fracture deformation in two- and three-dimensional media containing multiple intersecting fractures.

Torus computed tomography (Torus CT) reconstructs the Fourier coefficients of an image from X-ray data along closed geodesics on the flat torus. We extend Torus CT to experimental X-ray tomographic data by converting the data to a parallel-beam geometry and mapping the measured projection angles to the geodesic directions on the torus. In addition, we introduce two extensions of the Torus CT framework: the Star TCT method, which recovers additional Fourier coefficients beyond the standard frequency box using the same data, and torus backprojection (TBP), which sums the data over the measured directions. We prove that filtered TBP converges to the Tikhonov-regularized solution for general data up to an explicit error term that vanishes for data which is invariant along fixed geodesics. We also show that the direction-wise averaging of the zero mode is an unbiased estimator and its mean squared error decays like $N^{-2}$. We evaluate the methods on simulated phantoms and experimental X-ray data of a walnut using an updated implementation. The results show that torus-based reconstruction methods are applicable to experimental data, avoid the streak artifacts of filtered backprojection (FBP), and yield competitive reconstruction errors on both simulated and experimental data.

Statistical inverse problems have garnered significant attention in recent years due to the growing importance of statistical learning theory and functional analytic approaches in the fields of machine learning and artificial intelligence. In this paper, we investigate the stable approximation of the element $u^{\dagger}$ that satisfies the equation $Au = g$, where $A$ is a linear operator that maps a Banach space into an appropriate function space. The function $g$ is observed only through independently and identically distributed data points that are corrupted by noise and assumed to follow an unknown distribution $\rho$. We employ the Tikhonov regularization scheme, leveraging statistical learning techniques and the framework of reproducing kernel Banach spaces to estimate the solution. We establish convergence and derive the convergence rate of the estimated solution with respect to the true solution as the number of data points increases, with the rate expressed in probabilistic terms. The theoretical findings are further supported by numerical experiments that demonstrate the effectiveness of the proposed approach.

We present a Runge-Kutta method of weak order 2 for the numerical time integration of stochastic delay differential equations. This scheme extends the class of second order Runge-Kutta methods introduced by A. Rößler in [SIAM J. Numer. Anal., 47(3):1713-1738, 2009] for stochastic ordinary differential equations. The proposed integrator is applicable to equations with discrete commensurable delays and is particularly efficient for problems involving multiple noise terms. Experimental confirmation of the weak order 2 is provided and MATLAB codes are freely available.

We consider splitting methods for partial differential equations involving unbounded operators. For non-time-reversible dynamics, such as dissipative systems, negative splitting coefficients are generally not admissible because they require stepping backward in time, leading to an order barrier when all coefficients are required to be positive. We introduce generalized modified splitting methods to overcome this barrier. To analyze their convergence, we develop the corresponding multi-index series formalism, which provides a systematic framework for deriving order conditions. Using this formalism, we derive the order conditions and construct a generalized modified splitting method of order $6$. We also provide Python scripts that automate the generation and verification of order conditions, as well as the construction of new generalized modified splitting methods. Finally, we establish connections between the introduced multi-index series and related series formalisms from the literature, in particular, word series, Lie-Butcher series, and multi-index Butcher series.

In this paper, we introduce new approximation spaces and a locally adaptive refinement strategy for the hierarchical sparse-grid PIC (HSG-PIC) method to improve the bias while preserving the noise-reduction properties of sparse-grid methods. We first propose an energy-based approximation space, which optimizes the relation between the $\mathrm{H}^1$-norm error and the number of degrees of freedom, together with a family of generalized sparse-grid spaces that continuously connects classical sparse-grid and full-grid approximations. We then develop a locally adaptive approximation strategy based on hierarchical surpluses, combined with an efficient incremental refinement algorithm that avoids solving the Galerkin problem on the complete generalized space. Numerical experiments demonstrate that the proposed adaptive HSG-PIC method substantially improves the approximation of solutions with localized structures while maintaining the statistical advantages of sparse-grid discretizations. Compared with a standard full-grid PIC method, the adaptive approach achieves comparable or higher accuracy with a significantly reduced number of mesh nodes and particles. These results demonstrate the potential of adaptive sparse-grid PIC methods for efficient simulations of kinetic plasmas with complex solution structures.

For wave scattering by axisymmetric surfaces, the fast Fourier transform (FFT) method provides an effective tool to accelerate standard boundary integral equation (BIE) solvers. Surface BIEs can be decoupled into a series of curve integral equations on the generating curve, due to the convolution-like integral operators. The Fourier coefficients of the three-dimensional fundamental kernels can be rapidly computed through three-term recurrence relations based on Miller's algorithm. Such well-established techniques break down for nonaxisymmetric surfaces. This paper proposes a novel FFT-accelerated boundary integral method for wave scattering by smooth surfaces of arbitrary shapes. The Fourier coefficients of the singular kernels now satisfy higher-order recurrence relations. Although they can be solved with an optimal linear complexity by the standard Olver's algorithm, it turns out that a singularity swapping approach that rewrites each kernel as the product of a smooth function and an axisymmetric-related singular factor is realistically much faster. Consequently, Miller's algorithm together with the standard FFT convolution yields an ${\cal O}(M\log M)$ approach for evaluating the ${\cal O}(M)$ Fourier coefficients of the kernels, attaining exactly the same order of complexity for axisymmetric surfaces! With such FFT-based efficient procedures, we rewrite the surface BIEs in terms of ${\cal O}(M)$ curve integrals, which are proved to exhibit logarithmic singularities, discretize them by panel-based generalized Gaussian quadratures, and obtain spectrally accurate linear systems to approximate the wavefields. Extensive numerical experiments are carried out to demonstrate the effectiveness and spectral accuracy of the new approach.

We study resonance prediction in dispersive media, formulated as nonlinear spectral problems for volume integral operators. The main idea is to use asymptotic analysis not only as a baseline approximation, but also as a guide for constructing predictive correction models. We learn the residual between asymptotic and reference resonances using features suggested by the subwavelength expansion, including the logarithmic scales specific to two dimensions. The resulting corrections substantially improve single-resonator and dimer predictions, and symbolic regression produces compact formulas for the learned residual. The results show that asymptotic analysis can be used not only to approximate resonances, but also to design the feature space in which data-driven corrections become accurate, low-dimensional, and interpretable.

We develop and analyze a continuous/discontinuous Galerkin (CDG) method based on reconstructed surface Laplacians for the biharmonic equation on a smooth closed surface. Continuous mapped finite elements of degree $k\ge2$ are used on fitted parametric meshes of degree $r\ge1$, while a discontinuous degree-$(k-2)$ lifting corrects the broken Laplace-Beltrami operator for two-sided conormal-flux jumps. The resulting completed-square form is coercive on the mean-zero space for every fixed $\beta>0$, without requiring a sufficiently large penalty parameter. Under the standard geometric assumptions, we prove that the energy and reconstructed-Laplacian errors are $\mathcal O(h^{k-1}+h^r)$, and the $L^2$-error is $\mathcal O(h^{q_k}+h^{r+1})$, where $q_2=2$ and $q_k=k+1$ for $k\ge3$. Benchmark computations support these rates, while a surface Swift-Hohenberg experiment illustrates the extension of the method to nonlinear Laplacian-dominated models.

We study an inverse source problem for a linear non-autonomous parabolic equation with additive source term $f(t)+\zeta(t,x)$, where the unknown component depends only on time and is recovered from an integral observation of the solution. After reducing the problem to an equivalent linear inverse problem, we establish existence and uniqueness of the Tikhonov-regularized solution and derive a first-order optimality condition. We show that the forward operator admits a Volterra representation in time, yielding a weak-norm stability estimate in $H^{-1}(0,T)$ and, under an a priori $H^r(0,T)$ bound on the source, a conditional Hölder stability estimate in $L^2(0,T)$. Motivated by this Volterra structure, we introduce an adjoint-free integral feedback method that reconstructs the source using only forward solves. We analyze the feedback iteration by establishing its well-definedness and fixed-point properties, convergence for exact data, and finite-iteration stability with respect to noisy data. We further show that, with an appropriate noise-dependent stopping rule, the method constitutes an iterative regularization scheme. Numerical experiments for smooth and piecewise constant sources, supplemented by temporal regularization and automatic parameter selection, demonstrate accurate and stable reconstructions in the presence of noise.