Papers

Variational regularization techniques are dominant in the field of mathematical imaging. The expressiveness of these methods heavily depends on the adopted regularization formulation which usually involves a number of parameters which have to be fine-tuned by the user. To circumvent the heavy computational costs associated with this parameter tuning, data-driven learning strategies have turned out to be a useful strategy. Selecting parameters by learning from data gives rise to a hierarchical (bilevel) optimization problem, which is in general computationally hard to solve. Our aim in this paper is to develop a data-efficient solution strategy for the bilevel learning problem. To that end, we propose a multilevel Monte Carlo derivative-free optimization method that couples Gaussian smoothing of the upper-level objective (parameter selection) with a hierarchy of lower-level (variational regularization) solution accuracies. Being agnostic towards the specific way the lower-level problem is solved, we derive explicit oracle complexity and computational cost statements for approximately solving the entire bilevel learning problem. Our statements are parameterized by the employed lower-level solution method, thus allowing for a direct cost accounting from the variational reconstruction method to the parameter selection problem. This enables us to give a universal end-to-end complexity analysis of bilevel learning, valid for the most important optimization routines in inverse problems, rather than for a specific choice of methods.

The rectangular multiparameter eigenvalue problem (rMEP) is a recently introduced form of an eigenvalue problem, in which we are looking for the combinations of parameters where the multivariate rectangular matrix loses full column rank. Existing numerical methods for rMEPs typically require substantial memory and are not suitable for problems with a large number of eigenvalues. We propose a novel numerical approach based on the homotopy method that has low memory requirements and can compute all eigenvalues of an rMEP. Moreover, the homotopy method is naturally parallelizable, providing a simple and efficient way to accelerate computations. An implementation, together with a collection of test examples, is available in a public repository.

It is a classical fact in spectral geometry that the eigenvalues of the Dirichlet Laplacian depend continuously on the geometry of the underlying domain, under suitable geometric restrictions. We show that this theoretical property is reflected in a concrete numerical situation: the convergence rate of iterative solvers for convection-diffusion problems remains stable under deformations of the domain. The connection is a superlinear estimate for preconditioned GMRES/GCR in which the relative residual errors are bounded by an explicit expression of the eigenvalues of the Laplacian, or of a symmetric positive elliptic operator in divergence form. We prove that this bound depends continuously on the domain, uniformly in the iteration index, under either regularity hypotheses on the boundary or a topological constraint in dimension two. We also establish corresponding statements on the discrete level. Two examples delimit the theory: domains bounded by a fractal curve, where continuity still holds, and perforated domains, where it fails. The theoretical results, including both examples with rough domains, are confirmed by numerical experiments.

We present a mass-weighted meshless formulation for the shifted incompressible lattice Boltzmann method (LBM) whose streaming step is exactly unitary, combining with the Carleman embedding. In traditional meshless LBM, on a scattered node set the semi-Lagrangian streaming step is an interpolation, which is not unitary; a generic block encoding of this matrix therefore requires a normalization proportional to its sparsity and logarithmically many ancillas. To overcome the non-unitarity, we instead express the interpolation as the exponential exp(a Dx + b Dy) of mimetic derivative operators. Endowing each node with its Voronoi mass and applying a mass-weighted similarity transform renders these derivative operators skew-Hermitian by enforcing moment constraints together with low-wave number trigonometric differentiation conditions, so that the interpolation becomes exactly unitary. Such unitarization of interpolation leave the local collision step unchanged up to a rescaling of its quadratic term, where the same weights define the amplitude encoding of the distribution functions. For quasi-uniform node sets, the mass-weighting diagonal remains within a unit-order normalization that degrades only mildly as the Carleman truncation order grows, thereby mitigating the cost of block encoding. Numerical experiments on the Taylor Green vortex confirm that the method reproduces the full meshless LBM, with the velocity error decreasing as the Carleman order increases.

To solve a family of large linear systems where the system matrix depends on the family parameter quadratically, we propose a new solution approach. Our method employs the derivatives of solutions with respect to the parameter. It turns out that the accumulated solution derivatives of increasing order, evaluated at a single specially chosen parameter value, provide successively better approximations to solution subspaces for the whole linear system family. We show that the derivatives can be computed by a recurrence, involving, at each step, solving a linear system with the same matrix. Constructing an orthonormal basis for the subspace spanned by the derivatives then leads to an efficient solver where just a single (sparse) LU factorization (or a single preconditioner setup in an iterative scheme) suffices to solve the whole linear system family. We call the new method ODB (orthonormal derivative basis method). As shown in numerical tests, compared to conventional solution methods, the ODB method allows to drastically reduce computational costs.

We establish the holomorphic dependence of the Helmholtz Calderón projector on closed $C^{1,1}$ surfaces and of the operator associated with the Maxwell electric field integral equation (EFIE) on closed and open surfaces of the same regularity. More precisely, after transporting the boundary integral operators (BIOs) to fixed reference spaces, we prove that the resulting maps admit bounded holomorphic extensions to complex neighborhoods of compact sets of admissible boundary transformations, a property referred to as shape holomorphy. A key step in our proof consists in extending the results in [Lanza de Cristoforis and Rossi, Analytic Methods of Analysis and Differential Equations: AMADE 2006 (2008), pp. 193–220] and [Dalla Riva, Luzzini, and Musolino, Inverse Problems 38 (2022), 055004] from Schauder to fractional Sobolev spaces. For the EFIE, the single layer representation and the surface Piola transformation give the result in Maxwell trace spaces, whereas the result for open surfaces or screens is obtained by introducing suitable restriction and extension by zero operators. As opposed to existing results, we do not rely on an affine-parametric representation of the boundary deformations, thus allowing for more general parametric shape variations. Finally, we establish shape holomorphy of the shape-to-solution map for the EFIE, subject to interior electric Maxwell nonresonance on closed surfaces but with no such condition on screens. Under suitable summability and stability assumptions, these shape holomorphy results provide an analytic basis for the dimension-robust convergence of interpolation and quadrature methods used in computational uncertainty quantification (UQ) of scattering problems with high-dimensional, shape-parametric uncertain geometry.

Deep learning provides a flexible mesh-free framework for solving elliptic interface problems. However, domain-specific neural approaches often rely on separate networks for individual subdomains, causing the model complexity to grow rapidly as the number of subdomains increases. To address this limitation, we propose a Domain-Specific Adaptation Network (DSAN), which combines a shared backbone with lightweight subdomain-specific adapters to provide local represen- representation flexibility without replicating a complete network for each subdomain. We further establish a stability and error analysis showing that the continuous physics-informed residual controls the solution error, and decompose the resulting error into approximation, optimization, and sampling components. Numerical experiments on two- and three-dimensional problems with multiple regular and irregular interfaces demonstrate that DSAN achieves accuracy comparable to or better than fully decomposed multi-network methods while requiring substantially fewer trainable parameters. More- over, the parameter savings become increasingly pronounced as the number of subdomains grows, demonstrating the scalability of the proposed domain-specific adaptation framework.

Implicit time integration is the standard way to overcome the tight CFL time step restriction when solving nonstationary PDEs numerically. However, within an adaptive time accurate integration, the time step size selection is governed not only by stability but also by accuracy requirements. Therefore, an eligible question to ask is: is (unconditionally) stable implicit time stepping really needed if the time step is restricted by accuracy considerations anyway? And how much do we then lose (in terms of the time step size) when switching from the expensive implicit time stepping to stabilized explicit schemes? An answer to this question is also important because such a switch is often done when solving large scale nonstationary PDEs on massively parallel computers. This paper attempts to answer these questions by numerical tests for two particular adaptive time accurate settings: (i)~nonstationary convection-diffusion equations and (ii)~nonlinear heat equations.

In this work, long-time Sobolev stability and error estimates with optimal $\varepsilon$-dependence are established for the Strang splitting method applied to the Allen–Cahn equation with periodic boundary conditions. By combining the modified energy dissipation law, and Hölder-Gagliardo-Nirenberg product estimates, the numerical solution is proven to admit uniform-in-time $H^r$-seminorm bounds $\mathcal{O}(\varepsilon^{-r}+1)$, consistent with the corresponding Sobolev bounds of the exact solution. A second-order error estimate $\mathcal{O}(\varepsilon^2\tau^2)$ is further established on any prescribed finite time interval, exhibiting quadratic dependence on $\varepsilon$. Numerical experiments are presented to verify the optimal $\varepsilon$-dependence of both the stability and error estimates.

We study randomized row-action methods for problems with linear observations, convex constraints, and non-Euclidean geometry. Our starting point is a unified formulation in which the data-fitting term acts in observation space, while the feasible set and the regularizer encode structure in variable space. It covers constrained least squares, nonnegative and simplex-constrained models, regularized data fitting, and equality-constrained Bregman selection problems. The structure of these problems leads naturally to two algorithmic branches. The first uses primal Bregman-proximal row steps for soft data-fitting problems. The second uses dual Bregman–Kaczmarz steps for problems with hard equality constraints. We first derive a common one-step inequality. We then give deterministic and stochastic convergence results for the soft branch under relative smoothness. For constrained least squares, the normalized row oracle has exact first- and second-moment identities. They yield an expected ergodic bound in the inconsistent case, an $O(1/K)$ bound in the interpolating case, and linear convergence of the actual iterates under a residual error bound. For equality-constrained selectors, we prove linear convergence in expected Bregman distance under a residual-to-Bregman error bound. We recall global sufficient conditions for this bound and derive local constants for the nonnegative quadratic geometry and the entropy geometry on the simplex. Numerical experiments illustrate the noisy stepsize trade-off and the local entropy selector theory, and show that entropic row-action updates can recover interpretable factors in a simplex-structured nonnegative matrix factorization problem.

Geometrically exact beams naturally arise in physical phenomena and engineering practice for modeling flexible structures, such as soft robotic manipulators and DNA structures. To the best of our knowledge, there appear to be no geometric integrators for stochastic geometrically exact beams in a variational formulation. In this paper, we introduce variational principles for the dynamics of stochastic geometrically exact beams, in both continuous- and discrete-time settings. By introducing perturbations via stochastic potentials and stochastic kinematics, we derive the equations of motion from a variational viewpoint. The resulting stochastic geometric integrators simultaneously discretize time and space in a parallel manner, which provides a distributed computational framework for solving the stochastic dynamics of complex beam systems. In particular, the proposed integrators are shown to be a stochastic extension of the well-known deterministic staggered leapfrog scheme for the wave equation, and also a stochastic extension of Yee's scheme for spatially one-dimensional computational electromagnetism. Numerical results are shown to demonstrate the effectiveness of the proposed approaches.

We consider stable and effective restarts of the SS–RR and SS–RRR methods for computing the eigenpairs of a large Hermitian matrix $A\in\mathbb{C}^{n\times n}$ with the $n_{ev}$ eigenvalues in a given interval. Taking the CJ–SS–RR and CJ–SS–RRR methods for instance, we prove that the commonly used restart scheme and its orthonormalized variant, named the natural restart schemes, are not numerically stable, thereby lowering attainable accuracy of approximate eigenpairs and missing some of the $n_{ev}$ desired eigenpairs. To cure these deficiencies, we propose a new restart scheme, called the combination restart, that takes the $n\times \ell$ starting matrix to be the product of the matrix, which consists of the selected $n_{ev}$ approximate eigenvectors, and an $n_{ev}\times \ell$ Gaussian random matrix. We prove the remarkable mathematical and numerical superiority of the combination restart over the natural restarts, then present a hybrid restart scheme that fully exploits the advantages of the natural restart and combination restart, and finally develop CJ–SS–RR and CJ–SS–RRR algorithms with the hybrid restart. Numerical experiments confirm that the new algorithms are numerically backward stable to the working precision and more efficient than the CJ–SS–RRR algorithm with the natural restart. The combination restart and hybrid restart schemes can be straightforwardly extended to the SS–RR and SS–RRR methods for general standard and generalized eigenvalue problems.

Conventional control approaches to optimization problems under uncertainty entail a practical compromise: deterministic controls offer robustness at the expense of pathwise accuracy while stochastic controls can adapt to individual realizations but suffer from the curse of dimensionality. To achieve both robustness and adaptivity, a neural operator-based approach is proposed in this paper for distributed elliptic optimal control problems with high-dimensional random inputs. A data-driven training algorithm is first established to learn the mapping from coefficient realizations directly to their optimal pathwise controls. To bypass the burden of generating training data pairs, a physics-informed training algorithm is then proposed and recast as weakly supervised learning using pseudo label targets produced by standard adjoint-based control updates. Numerical experiments with rough coefficients are conducted to validate the effectiveness and efficiency of our robust-pathwise controls, which outperforms the benchmark robust deterministic control on unseen coefficient realizations owing to the generalization capability and fast inference of trained models.

Many applications in science and engineering require solving sparse Hermitian eigenvalue problems that are so large that even algorithms whose runtime scales linearly with the matrix dimension are unacceptably slow. Responding to this challenge, there is a growing body of work on sparsified iterative methods and submatrix selection methods that often perform well in practice. However, theoretical understanding of these methods is limited. In this paper, we take a step toward better understanding sublinear-time eigenvalue computation by analyzing the sparsified power method. We prove that the method approximates a leading eigenpair to a prescribed accuracy in a number of iterations comparable to the standard power method. The required number of retained nonzeros depends explicitly on the target accuracy, the initial overlap with the leading eigenvector, the compressibility of the leading eigenvector, and the relative spectral gap. When these quantities and the column sparsity are independent of dimension, the runtime is independent of dimension. In contrast, we show that a comparable worst-case guarantee cannot hold for a prototypical sparsified submatrix selection method. Submatrix selection is the approach of choice in many quantum chemistry applications. However, we show that, even with substantial initial eigenvector overlap, a single iteration of the method can produce a vector that is orthogonal to the leading eigenvector, after which the method cannot recover. Nevertheless, experiments on matrices from the SuiteSparse matrix collection and quantum chemistry problems illustrate the practical effectiveness of submatrix selection methods in a variety of domains despite these worst-case limitations.

We consider an elliptic Cauchy problem for the Poisson equation in which both Dirichlet and Neumann boundary data are prescribed only on a portion of the boundary, while no boundary conditions are available on the remaining part of the boundary. Such problems are severely ill-posed. In this work, we study the convergence of an iterative method based on the alternating imposition of Dirichlet and Neumann boundary conditions. We formulate the scheme within a Hilbert space framework, where a yet another alternation becomes apparent: each iteration can now be interpreted as an orthogonal projection onto an appropriate closed subspace. This perspective provides a transparent geometric interpretation of the algorithm and yields a concise convergence proof based on von Neumann's classical result on alternating (orthogonal) projections. In addition to elucidating the convergence of the method, this alternating projections perspective also offers a general rate of convergence in terms of the Friedrichs angle (between two closed subspaces). We point out by means of an example that such convergence rate estimate degenerates for the elliptic Cauchy problem under consideration.

In this work, we develop a general memory-augmented phase-field (MAPF) framework based on Onsager's variational principle, in which memory effects are incorporated through a general relaxation kernel while preserving the thermodynamic consistency of the underlying system. Rather than treating fractional models as fundamental formulations, we show that they can be interpreted as reduced descriptions associated with specific power-law memory kernels. The proposed MAPF equations provide a unified modeling framework for non-Markovian phase-field dynamics. To facilitate analysis and computation, we further introduce an equivalent internal-variable representation that transforms the original nonlocal-in-time dynamics into an extended Markovian gradient system. This formulation preserves key structural properties, including energy dissipation and mass conservation, and enables the development of efficient structure-preserving numerical schemes. In particular, we construct an energy-stable scalar auxiliary variable approach combined with a fast rank-one implementation that exploits the relaxation-spectrum-induced coupling structure, resulting in a computationally efficient algorithm with reduced storage requirements. The proposed framework establishes a systematic connection between memory-dependent phase-field modeling, internal-variable thermodynamics, and structure-preserving numerical computation. Numerical experiments demonstrate the accuracy, stability, and efficiency of the proposed methods for phase-field dynamics with long-time memory effects.

We establish inverse inequalities in Gaussian Sobolev norms for shallow $\mathrm{ReLU}^k$ trial spaces. The estimates identify the sharp regularity threshold and quantify the dependence on parameter separation and bias magnitude. Based on these estimates, we prove convergence of weighted least squares discretizations of linear PDEs on $\mathbb R^d$ using shallow networks with prescribed inner parameters and finitely many residual samples in an expanding ball. Under suitable stability and regularity assumptions, we obtain an error bound of order $h^{k+(d+1)/2-r}\delta_\pm^{-(S-r)}$ in the Gaussian Sobolev norm of order $S$, up to logarithmic factors, where $r$ is the PDE order, $h$ is the parameter fill distance, and $\delta_\pm$ is the antipodal separation. The convergence results cover deterministic and random sampling and extend to fixed bounded regions by restriction. We construct parameter sets with $\delta_\pm\asymp h$, yielding a convergence rate in terms of the network width that matches a lower bound up to logarithmic factors. Numerical experiments support the theoretical results and provide practical guidance on balancing accuracy and stability.

We develop and analyze a symmetric two-level preconditioner for the method of reflections for Laplace problems in perforated domains. Using constant boundary modes, we define a coarse correction through a symmetric coercive dual formulation. Combining this coarse correction with parallel reflection corrections gives a symmetric positive definite preconditioner for solving the dual system by the conjugate gradient method. Under a geometric separation condition, we prove a condition number bound independent of the number and sizes of the holes and the diameter of the ambient domain. We also derive geometry-dependent convergence estimates for the one-level methods that explain the need for a coarse correction.

In this work, we develop and analyze a splitting mixed finite element method, combined with a fully implicit Euler time discretization, for a stochastic Keller–Segel system with general polynomial logistic growth and multiplicative noise. The proposed splitting formulation allows the use of standard finite element spaces without requiring the LBB inf–sup condition associated with traditional mixed finite element methods and, at each time step, avoids the need for solving a coupled nonlinear system. We first establish the well-posedness of the variational solution and derive new high-moment, higher-order regularity, and exponential stability estimates. These estimates provide the key analytical foundation for the strong convergence analysis of the numerical method. A principal contribution of this work is the derivation of error estimates directly in full expectation, without the use of localization arguments. By combining the new stability and regularity estimates with a stochastic Gronwall framework, we overcome the difficulties arising from the interaction between the nonlinear chemotactic coupling and the multiplicative noise. As a result, we establish strong convergence in high-moment norms and obtain optimal convergence rates with respect to both the temporal and spatial discretizations. In particular, the proposed method achieves the optimal rates of order $k^{\frac{1}{2}}$ in time and $h$ in space in the relevant $L^2$- and $H^1$-based error norms. These results provide, to our knowledge, the first optimal strong convergence analysis in full expectation for a splitting mixed finite element approximation of this class of stochastic Keller–Segel systems.

We extend the study of Hamilton-Jacobi equations on networks from a single equation to a weakly coupled system of two equations on a star graph, and we develop a numerical scheme for it. Instead of prescribing the boundary values and imposing a junction condition, we compute the solution on each branch from its outer vertex toward the junction and determine the boundary values for which these branchwise solutions agree at the junction. On each branch, we replace the right-hand sides by their cell averages and solve the resulting problems exactly on every cell. This gives an iteration whose limits we compute in closed form. From these limits, the values of the discrete solution at the junction are affine functions of the boundary values on the branch, with a $2 \times 2$ matrix whose rows sum to one and whose determinant factors into explicit contributions of the cells. We use these formulas to characterize the boundary values that give solutions continuous at the junction and monotone on every branch: when the matrices are nonsingular, the admissible boundary values of one branch form a convex polygon, bounded by conditions on the difference of the two components and by positivity conditions, and the boundary values on the other branches follow from them. We also treat the singular case, and we illustrate the results with numerical experiments.

The modified Anderson-Bjorck (modAB) method was first proposed in 2023 by Ganchovski and Traykov for solving nonlinear equations in structural mechanics. It was improved further in 2026 by introducing better switching criteria with account to function symmetry, and bisection fallback based on measuring the bracket reduction rate. Although the algorithm demonstrated superior performance across a great variety of functions, it lacked formal mathematical analysis and proof of convergence. In addition, some discrepancies were noted by L. Tomov in the latest version that need to be addressed: premature application of the Anderson-Bjorck (AB) correction on the first AB step, caused by retaining the endpoint-update record from the preceding bisection step, and using the corrected endpoint ordinates in symmetry factor calculation instead of the true residuals. The purpose of the present paper is to provide rigorous mathematical proof of global convergence in exact arithmetic, remove the above discrepancies, and provide floating-point safeguards for interpolation, switching and stopping, with explicit treatment of infinite and NaN values, underflow, overflow and round-off errors. These safeguards do not certify the accuracy of the supplied function values. Additional performance improvement was introduced by allowing for early bisection fallback when AB fails to make progress on reducing either the bracket width or the residual, while retaining a finite bound on consecutive exceptions to the prescribed width condition.

We present several new and surprising results on the theory and practice of Arbitrary Precision Approximate (APA) matrix multiplication, which refers to algorithms that output merely approximations to matrix multiplication even in exact arithmetic. While APA algorithms have been around for nearly half a century, and despite the fact that they have potentially lower (asymptotic) cost than many fast exact matrix multiplication schemes, they have not been studied as closely, to the extent that relatively little is known about their numerical properties. In this paper, we shine a spotlight on these algorithms by presenting three important results. First, we provide a new search algorithm for APA routines that is easy to use, automatically generates encoding matrices with functional entries, and can match or beat the best known APA algorithms for several sizes. In particular, we improve the upper bound on the border rank of $\langle 2,4,5 \rangle$ matrix multiplication. Second, we combine this insight with alternative bases to minimize the number of additions in an APA algorithm, thereby improving its leading coefficient. Finally, we prove the same kind of stability that characterizes fast exact matrix multiplication for an algorithm of Bini [Calcolo 1980] that computes an exact product of two matrices using a modified APA scheme (with or without alternative bases).

Gaussian process solvers for linear partial differential equations (PDEs) return posteriors rather than point estimates, but they condition only on equality information: the PDE residual, boundary data, and measurements. Many solutions of interest also obey inequality constraints, including nonnegativity, monotonicity, and absence of arbitrage, which equality conditioning cannot enforce. We prove that conditioning a Gaussian process prior jointly on linear equality observations and a finite-dimensional inequality event is equivalent to conditioning sequentially in terms of Borel measures on the solution space. The physics-informed posterior is therefore a valid prior for the inequality update at virtual points. The combined law is non-Gaussian yet computable through truncated-Gaussian moments, and neither the PDE conditioning nor the inequality update requires modification. We evaluated the construction on a Poisson benchmark, option-price surfaces, and air-quality fields. Under incomplete boundary data, a monotonicity constraint lowered the error of a two-dimensional Poisson reconstruction by $63\%$. On a real S&P 500 option chain thinned to $20$ retained quotes, no-arbitrage constraints reduced held-out extrapolation error by $22\%$ and cut the calendar-arbitrage rate from $10\%$ to $1\%$. On four real air-quality snapshots, nonnegativity reduced a Gaussian moment-based sub-zero-tail proxy by $74\%$ at the sparsest monitor count. Changes in mean accuracy and constraint diagnostics differed across data regimes. The probability of the noisy virtual constraint event under the equality posterior measures agreement with the specified bounds before the update. The construction turns admissibility from a post-processing step into conditioning information for probabilistic PDE solvers.

We consider an inexact version of the Schwarz domain decomposition method in $3$D to solve a Laplace problem in a domain that is the union of two overlapping spheres. We solve each subproblem approximately by projecting its boundary data $L^2$-orthogonally onto the space of spherical harmonics of degree at most $N$. The model problem arises from the ddCOSMO solver used for implicit solvation in computational chemistry. We analyze the Schwarz-spherical domain decomposition method by establishing a new variant of the maximum principle.

We consider the Helmholtz equation in a Dirichlet angular sector (which may be locally perturbed) covered on one side by a thin homogeneous layer. First, we propose an asymptotic expansion of the solution at any order and for any angle when the thickness of the layer tends to 0, using a matched asymptotic expansion. This consists in introducing different asymptotic expansions in three subdomains – the vicinity of the corner, the layer and the rest of the domain – that are linked through matching conditions. The presence of the corner makes these matching conditions delicate to derive because the fields have complex singular behaviors. We overcome this difficulty using a novel system of algebraic tools, including formal series and operators on these series. Moreover, we justify the asymptotic expansion with error estimates and provide details on the form of the expansion depending on the angle. We numerically validate this expansion up to order 5. This requires computing the behavior at infinity of some corner profiles. To do this, we introduce a Galerkin-type method enriched by complements at infinity. These complements are known functions that appear in the behavior at infinity of the corner profiles with coefficients that become unknowns of the system. Finally, we propose effective models of order 3, resp. 5 where, far from the corner, the layer is replaced with an effective Robin, resp. Ventcel, boundary condition and, near the corner, an artificial boundary is introduced and a Dirichlet-to-Neumann (DtN) type condition is imposed. The DtN map relies on a modal decomposition of the truncated asymptotic expansion at order 3, resp. 5. The modes take into account the singular behavior of the expansion terms, the corner profiles and some algebraic matching coefficients. We prove the well-posedness of our effective models, as well as error estimates, and provide numerical results.

An inverse spectral problem for a $(1+1)$-dimensional wave equation with wave speed $c(x)$ is considered. The spectral data consist of the eigenvalues of the Neumann–Dirichlet problem for the operator $-c^2(x)\partial_x^2$ together with their derivatives with respect to a Robin parameter at the boundary. We propose a two-step recovery procedure. In the first step, the Neumann-to-Dirichlet map $\Lambda$ is reconstructed from the spectral data using the results of arXiv/0406582. In the second step, the reconstructed $\Lambda$ is used within the regularization strategy for the boundary control method presented in arXiv:1509.04478. As an improvement, the regularization of the control problems employs a weighted combination of $L^2$ and $H^1$ penalty terms, and we prove a uniform convergence result for the family of optimal controls, thereby allowing the regularization parameters to be tuned simultaneously. We also prove that the forward spectral map, which assigns to a wave speed its spectral data, is Lipschitz continuous, and that the inverse spectral map underlying the reconstruction procedure is Hölder continuous on a subclass of wave speeds. The method is implemented numerically and validated on both smooth and piecewise smooth wave speed profiles. Numerical experiments demonstrate that the recovery of the wave speed exhibits high accuracy and robustness.

We present a fully discrete energy conserving discretisation of the rotating shallow water equations with arbitrary order in space and time. A space time finite element discretisation is used, built from compatible finite element methods in space and continuous Petrov-Galerkin schemes in time. Projected auxiliary variables for Hamiltonian derivatives are necessary for exact energy conservation. We provide an iterative solver approach based on monolithic patch smoothers, which allows us to efficiently converge a Newton method sufficiently that energy errors due to the solver are very small. We extend this construction to the energy conserving Stochastic Forcing by Lie Transport formulation. Then we demonstrate the convergence, conservation and robustness in numerical tests.

The proximal Galerkin algorithm is a structure-preserving solver for variational inequalities that results in discretization-independent nonlinear iteration counts. However, in large three-dimensional simulations, the arising Newton saddle point systems cannot be feasibly solved by direct methods. To leverage an iterative linear solver, we propose a robust and scalable block preconditioner. Instead of using the usual dual Schur complement that acts on the latent variable, we design a preconditioner for the primal Schur complement. This leads to a block preconditioner that requires two approximate inverses of the latent block and only one of the primal Schur complement. The latent block represents an algebraic equation which, after a small perturbation, can be cheaply inverted. We then introduce two choices of approximations for the primal Schur complement: a PDE-operator-based and an algebraic approach. Either approximation can be efficiently inverted by geometric multigrid. We investigate the effectiveness of these strategies when proximal Galerkin is applied to obstacle, Signorini, and strain-constrained elasticity problems and also compare with a monolithic multigrid approach. Choices of standard relaxation methods, in the block preconditioner route, lead to robust iteration counts with respect to the mesh size and discretization degree. Moreover, in the context of the Signorini, specialized relaxation and prolongation also result in effective solvers for nearly-incompressible elastic materials. The linear solves are parameter-robust in 2D, with only mild growth in 3D. This allows for the solution of a three-dimensional contact problem, with a Poisson's ratio of 0.499, discretized with over ten million degrees of freedom.

Spherical harmonics have been studied extensively and are used in many fields of application. We consider the probably earliest application of spherical harmonics, namely for a representation of the solution of the Laplace equation on the $d$-dimensional unit sphere where the boundary data are represented in a spherical harmonics basis. We consider general $d \geq 2$. For the solution $u$ of the Laplace equation the fundamental maximum principle holds. If one approximates the boundary data in the finite dimensional space spanned by all spherical harmonics of degree at most $N$, the resulting discrete solution, denoted by $u_N$, does in general not satisfy the maximum principle. In this paper we derive sharp estimates for $r_N^{\rm max}$ such that the maximum principle holds for $u_N(x)$ for all $x \in \Bbb{R}^d$ with $\|x\| \leq r_N^{\rm max}$. We show that for $d=3$ one has $r_N^{\rm max} \gtrsim 1- 3 N^{-1} \ln\big( 2^\frac13 N \big) $. Similar estimates are derived for $d \neq 3$. We also study the sharpness of these estimates.

In this paper we propose a novel regularisation strategy for structured linear inverse problems of the form $Ax=F \circ G x = z$. Instead of using the operator $A$ or a complex scheme involving both $F$ and $G$ simultaneously, combining already existing and simpler regularisations for $F$ and $G$ allows the use of established and optimised programs for high-speed computations. As classical parameter choice rules fail in the scenario when no source conditions are known because intermediate error bounds may be unavailable, we provide both a-priori and a-posteriori rules with a convergence speed analysis given source conditions. We test the proposed method on the pre-clinical medical imaging technique of magnetic particle imaging with different choices of regularisation and on a convolution problem for an analysis of the convergence speed.