Papers

In this paper, we propose an unsupervised, two-step model-informed deep learning framework for solving the inverse source problem. In the first step, boundary measurement data are transformed into an imaging function that encodes information about the shape and location of the unknown source. Leveraging this information, we derive a Fourier-based model equation in the spectral domain that relates the imaging function to the Fourier coefficients of the source function. In the second step, this model equation is incorporated into the training of a model-informed neural network to recover the parameters of interest. The proposed approach enables fast and accurate reconstruction of the source function while maintaining robustness to noise. The effectiveness of the method is demonstrated through numerical experiments in both two- and three-dimensional settings. For the two-dimensional case, we further compare our approach with the traditional least-squares method to validate its computational efficiency and reconstruction accuracy.

On unstructured grids, the tangent of hyperbola for interface capturing (THINC) method with quadratic surface representation and Gaussian quadrature (THINC/QQ) requires topology-dependent numerical quadrature for cell and face integrals and Newton iteration for the surface constant. This study introduces a quadratic moment-matched sigmoid (QMMS) approximation that evaluates these quantities from the mean and variance of the quadratic field while retaining the hyperbolic-tangent reconstruction profile and the boundary variation diminishing (BVD) formulation. A slope-matched Gaussian cumulative distribution function and a moment-matched Gaussian model provide closed form approximations for the surface constant and face averages, eliminating Newton iteration and topology-dependent runtime quadrature from these evaluations. Single-cell tests and nine two- and three-dimensional benchmarks assess reconstruction accuracy, flow-field differences, and computational cost. Across the tested cases, the principal transported and shock structures remain closely aligned with those obtained using THINC/QQ, with local differences concentrated mainly near discontinuities, contact regions, and developed shear layers. Reconstruction-stage profiling gives mean THINC/QQ-to-QMMS time ratios of 2.12 in two dimensions and 1.85 in three dimensions. QMMS therefore reduces reconstruction-stage cost while replacing topology-dependent runtime integral evaluation with a common moment-based formulation

This paper analyzes finite element approximations for Maxwell's equations with time-dependent coefficients and nonlinear E-H coupling in both the total current density and the Silver-Müller-type boundary condition. The nonlinear coupling and the time-varying material parameters prevent the use of standard analytical tools, including the discrete compactness property and the Minty-Browder argument, making the numerical analysis particularly challenging. As the main novelty, we establish two primary contributions: Uniform convergence on every finite time interval $[0,T]$ and exponential stability on the infinite time horizon $[0,\infty)$. The uniform convergence is established via a Cauchy-type argument within a time-dependent Hilbert space framework that entirely bypasses Minty's trick and the discrete compactness argument. Furthermore, under a positive homogeneity assumption on the nonlinearities, the exponential stability is analyzed by means of a nonlinear sup-max problem arising from ratios of consecutive discrete energies. With this strategy, we ultimately prove a fully discrete unconditional exponential decay result under a specific sufficient condition. In the absence of this condition, a concrete counterexample shows that the energy decay fails to hold. Numerical experiments are provided to validate the theoretical findings.

The long-run statistics of chaotic dissipative partial differential equations are described by invariant measures, and when the equation carries a parameter these measures form a family whose members can differ in kind, from steady states to sustained spatiotemporal chaos. We ask whether a single generative model can cover such a family. We train a single conditional diffusion model once, on states of the Kuramoto-Sivashinsky equation at 160 values of its hyperviscosity, and evaluate it against held-out data at interleaved parameter values not seen in training, using a classifier two-sample test supported by spectral, geometric and tail statistics. Over most of the family, agreement improves as the dynamics become more chaotic. We identify a structural cause: the weakly chaotic members have invariant measures concentrated on low-dimensional sets, which a model with full support output cannot represent, and the target approaches full dimensionality as chaos increases. Restricting training to the strongly chaotic members reduces the classifier's excess over chance on those members by roughly 40%, and this restricted conditional model matches or exceeds models trained at single parameter values. On those members, generated samples reproduce energy spectra to within a resampling floor, with a small systematic deficit at the largest scales, and match the energy balance of held-out data.

A clear link is established between conventional numerical methods and neural network approximations for solving one-dimensional scalar conservation laws. The focus is on the construction of an appropriate flux term in the case of convex flux functions for improving the classical schemes. The first neural network developed here is able to rediscover Godunov's method, while the second one emulates the behavior of a second-order slope-limiter function. In this way, by merging them, second-order reconstruction-based schemes can be developed. The networks presented here employ a minimal number of parameters, significantly reducing the complexity compared to previous approaches. These networks can also be linked consecutively to get a deep one corresponding to multiple time steps. Training them with an appropriate loss leads to stable schemes, improving even the classical methods without increasing their complexity.

In this paper, we investigate how adaptive time-integration strategies can be effectively combined with parallel-in-time numerical methods for solving systems of ordinary differential equations. Our focus is particularly on their influence on tolerance proportionality. We examine various grid-refinement strategies within the multigrid reduction-in-time (MGRIT) framework. Our results show that a simple adjustment to the original refinement factor can substantially improve computational stability and reliability. Through numerical experiments on standard test problems using the XBraid library, we demonstrate that parallel-in-time solutions closely match their sequential counterparts. Moreover, with the use of multiple processors, computing time can be significantly reduced.

We compute $\E[(S_T-K)^+]$ by Monte Carlo for a scalar stochastic-volatility model with a fast mean-reverting factor of time scale $\eps$, for $\eps$ ranging from $1$ down to $10^{-3}$. A conditional (mixing) estimator gives finite variance, whereas the direct estimator has infinite variance for this model. The volatility factor is simulated with its exact Ornstein–Uhlenbeck transition. As $\eps\to0$ the price converges, at rate $O(\eps)$, to the Black–Scholes price with the averaged volatility $\bar\sigma$, and the implied-volatility smile flattens to $\bar\sigma$. Finally, we test a martingale control variate built on the Black–Scholes delta with volatility $\bar\sigma$. If the martingale is driven by the true volatility $\sigma(Y_t)$, the variance is reduced by a factor that grows like $1/\eps$, about $160$ at $\eps=10^{-3}$. If it is driven by the constant $\bar\sigma$, the variance is essentially not reduced. Last, we calibrate the model with $\rho\neq0$ to S&P~500 implied volatilities by full simulation of $S$. At the same number of paths the control variate reduces the variance by a factor $1.7$–$16$ (median $5.6$) and gives more accurate calibrated parameters; at equal CPU time it pays off only if the delta is rebalanced on a coarser grid than the time step. For a basket of four indices (state dimension $8$) the gain is larger (median $8$) and the relative cost smaller, so that the control variate is about $4$–$5$ times more efficient than plain Monte Carlo at equal CPU time.

This work considers strong and weak stochastic Galerkin approximations of random domain problems for the elliptic random partial differential equation (PDE). The random domain problems are realized by sufficiently regular random domain mapping, allowing a transformation on deterministic reference domain. A traditional numerical method for solving the resulting high-dimensional coupled stochastic Galerkin systems is replaced by deep learning techniques. We compare a physics-informed neural network approach, based on the strong stochastic Galerkin residual, with a Deep Ritz approach, based on the weak stochastic Galerkin system, reformulated as Ritz energy minimization. The neural networks serve as surrogates for the deterministic spectral coefficients of the respective stochastic Galerkin solution. The stochastic parameters are not used as neural network input parameters and the stochastic domain is not sampled during training. The efficiency of the methods is demonstrated on a randomly stretched interval and a randomly deformed annulus in two spatial dimensions.

We present a fully GPU-oriented SVD pipeline based on polar decomposition, motivated by iterative methods that rely solely on matrix multiplications, such as the Newton-Schulz iteration. We show that this approach enables up to a $2\times$ speedup compared to standard implementations. Furthermore, we derive a numerically stable backward pass for the polar decomposition and leverage it to obtain a fully differentiable SVD. Our methods are released as open-source implementations in both PyTorch and JAX: https://github.com/fallnlove/cans_svd.

Bayesian state estimation for high-dimensional nonlinear dynamical systems entails a fundamental tension between statistical fidelity and computational tractability, as particle weights can collapse, while Gaussian ensemble updates can miss non-Gaussian posterior structure. Score-based diffusion filters offer a sampling-based alternative, but existing training-free score filters often rely on heuristic likelihood corrections, which can compromise posterior accuracy by neglecting uncertainty about the system state associated with each noisy reverse particle. To address these issues, we propose AECSF, a training-free adaptive ensemble conditional score filter. AECSF constructs an analytically tractable score estimator from the conditional Tweedie identity, which recasts noisy posterior score estimation as estimating the conditional mean of the system state given a noisy reverse particle and the observation. To estimate these conditional means efficiently, AECSF employs a shared adaptive weighted proposal ensemble, while particle-specific conditional weights yield an estimate for each noisy reverse particle without separate proposal sampling. The proposal ensemble is updated using reverse-particle information within the same reverse-diffusion run to improve conditional-mean estimation. Theoretically, we characterize when a fixed weighted proposal measure yields the exact noisy posterior score. Under stated assumptions, we establish a bound relating conditional-mean estimation errors to reverse-sampling endpoint error. Numerical experiments demonstrate that AECSF improves the accuracy of posterior sampling and nonlinear filtering in high-dimensional problems with limited forecast ensembles.

Finite-difference discretizations of Laplace operators yield discrete Laplacians whose spectral properties are closely tied to the stability, convergence, and physical fidelity of numerical solvers. While symmetric discretizations are well understood, many high-order finite-difference schemes produce nonsymmetric matrices for which rigorous spectral analysis is overlooked. Surprisingly, we prove that two fourth-order schemes with Dirichlet boundary conditions yield nonsymmetric discrete Laplacians whose spectra are purely real and strictly positive. Computer-assisted computations show that this property does not hold for higher-order extensions of these schemes in general.

We develop a structure-preserving locally optimal block preconditioned conjugate gradient(LOBPCG) method for computing several of the largest positive generalized singular values of a matrix pair \((A, B)\), where \(B\) has full column rank, together with the associated generalized singular vectors. The method is derived from the Hermitian-definite Jordan–Wielandt pencil and therefore avoids explicitly forming the normal-equation matrices. By exploiting the block structure and spectral symmetry of this pencil, we construct separate search subspaces for the left and right generalized singular vectors. The resulting Rayleigh–Ritz procedure requires only a singular value decomposition of a small projected cross matrix. We adapt the improved Hetmaniuk–Lehoucq (IHL) trick to construct stable conjugate search directions in the two component spaces. We also establish an explicit correspondence between the coefficient matrices produced by the componentwise and augmented Rayleigh–Ritz procedures. This correspondence permits the componentwise IHL construction to be used in an augmented, structure-preserving LOBPCG formulation without explicitly assembling the augmented basis. Numerical experiments demonstrate the computational advantages of the proposed method over unstructured and Gram-matrix-based LOBPCG implementations.

Geometrically parametrised PDEs arise in shape optimisation, where solving a PDE over many deformed domains at manageable cost is a central challenge. Free-form Deformation (FFD) parametrises shape variations through a lattice of control points, but the resulting parameter space is typically high-dimensional, with most directions barely affecting the quantity of interest (QoI). We address this with a two-stage reduction pipeline. An Active Subspace (AS)-guided structured control-point selection identifies the low-dimensional subspace of FFD weights driving the QoI. A Sinkhorn Autoencoder (SAE) is then trained on the reduced weights, learning their distribution in a compact latent space by minimising the Wasserstein distance between aggregated posterior and prior. We compare three non-intrusive Reduced-order Models (ROMs), mapping full FFD weights, AS-reduced weights, or SAE latent codes to the QoI, each using Random Forests, Gaussian Processes, and K-Nearest Neighbours. We test the methodology on a non-linear heat equation with Gaussian random Robin boundary data on a deformed Stanford Bunny, and on drag prediction for the DTC Hull bulb under FFD deformations. In both cases, the generative model captures the bulk of the QoI distribution over a comparable support, thinning the shoulder opposite the design objective, while its extremes in the design-relevant direction match the FFD sample within numerical and sampling uncertainty. ROMs learned in the SAE latent space yield smaller test errors than those on full FFD weights in most regressor-metric combinations, whereas AS-reduced weights give no comparable gain. Only the latent representation attains a positive predictive coefficient $Q^2$.

Solution of convection-dominated problems using stabilized methods often requires to specify stabilization parameters whose optimal choice is not known but which considerably influence the quality of the approximate solution. In a recent work of the authors, an approach for computing these parameters by machine learning was proposed and applied to the streamline upwind/Petrov-Galerkin (SUPG) method for convection-diffusion equations. In the present paper, this approach is studied numerically for unstructured meshes and extended to meshes obtained by anisotropic mesh adaptation.

We study the problem of implementing multivariable matrix-valued functions on a quantum computer and propose two quantum algorithms for multivariable matrix eigenvalue transformations acting on tuples of pairwise commuting matrices. The first algorithm is based on multivariable contour integrals and applies to arbitrary holomorphic functions on $\mathbb{C}^n$. The second algorithm is based on high-dimensional Fourier transforms and designed for smooth functions of commuting Hermitian matrices. For both algorithms, we discuss and analyze the complexity of their quantum implementation leveraging quantum singular value transformation, compression gadgets, and linear combination of unitaries. We additionally study their hybrid quantum–classical variants that reduce the number of required ancilla qubits to logarithmic dependence on the number of variables, at the cost of increased query complexity. As an application, we show that multivariable matrix polynomials can be implemented with no additional explicit degree dependence, depending only on global properties of the polynomials.

Operator learning of periodic solutions requires deciding how simulation data should be recorded and represented. A natural choice is to integrate long enough for transients to decay and record a window wide enough to contain at least one full period of all trajectories. We find that these conservative choices can make the resulting trajectories difficult to learn, even when the underlying periodic orbits vary regularly with system parameters. Unaligned trajectories generalize poorly even within the training distribution. Phase alignment substantially improves in-distribution generalization, but models trained on a fixed physical-time window still have large errors on trajectories with periods outside the training range. We explain both failures through a common mechanism: frequency differences accumulate over time, so the target phase varies rapidly with the parameters. Predictors that cannot track this variation incur a population MSE floor in both settings; for fixed window prediction, we also derive a per-sample lower bound. We then study one of the simplest representations that escape these floors: learning an aligned, normalized waveform and its period separately. We establish regularity of the decoupled targets under ODE assumptions and show experimentally that this approach avoids both failures in ODE systems and a PDE case study.

We develop entropy-stable isothermal-wall conditions for the compressible Navier–Stokes equations at continuous and semi-discrete levels. For a stationary, impermeable, no-slip wall at constant positive temperature, adding total-energy density divided by wall temperature to the canonical convex entropy cancels the wall contribution without suppressing heat transfer. This affine modification preserves the entropy Hessian, relative entropy, and volume entropy production. It is proportional to total ballistic-energy density and yields a conditional $L_2$-type bound based on the initial field under the stated assumptions. Two wall formulations use a common discontinuous Galerkin discretization with the summation-by-parts property and simultaneous approximation terms. Both give zero adapted wall contribution without a penalty and a nonpositive contribution with it. The complete-residual formulation penalizes the full wall-data residual and balances fluid energy through its complete numerical wall flux. The Fourier-energy-exact formulation combines an entropy-neutral skew correction with a momentum-projected penalty so that the complete numerical outward energy flux equals the one-sided numerical Fourier heat flux at every wall node and spatial resolution. Pointwise tests verify the local identities. Manufactured-solution studies with entropy-stable interior coupling demonstrate convergence rates consistent with order $p_s+1$ in the discrete $L_2$ norm, where $p_s$ is the solution polynomial degree. Thermal-relaxation calculations verify the semi-discrete energy and adapted-entropy balances. Three-dimensional calculations demonstrate balance closure in a closed domain, convergence of curved-wall traces, and applicability to a complex supersonic separated flow.

We establish that finite expressions form a symbolic representation class capable of overcoming the curse of dimensionality for several classes of high-dimensional partial differential equations. For semilinear heat equations, we show how finite expression approximations of the terminal condition and nonlinearity can be propagated through the multilevel Picard framework to produce randomized pointwise approximations with prescribed root-mean-square accuracy. For semilinear Kolmogorov equations with Laplacian diffusion and zero drift, diagonal Black–Scholes equations, and the Laplace Dirichlet problem on a half-space, we construct deterministic finite expression approximations with arbitrarily small spatial $L^p$ error. Under suitable growth and regularity assumptions, the evaluation cost of the resulting approximants is bounded polynomially in the dimension and the reciprocal accuracy. A key ingredient in the nonlinear setting is the construction of finite expressions that approximate the PDE solution while preserving the growth and Lipschitz structures required by the stochastic solution theory. More broadly, our results show that dimension-robust approximation of high-dimensional PDEs is not restricted to conventional neural-network architectures: structured symbolic expressions generated from a fixed dictionary can achieve comparable polynomial-complexity guarantees. This provides a rigorous foundation for finite expression methods as a representation paradigm for high-dimensional scientific computing.

We develop a high-order and positivity-preserving kinetic solver for the multi-species Vlasov-Ampère system that is built within a semi-Lagrangian discontinuous Galerkin (SLDG) framework. The method is designed for efficiency and accuracy across a broad range of conditions, including stiff, multi-scale regimes with strong ion-electron scale separation. Phase space is discretized via high-order discontinuous Galerkin finite elements; time stepping is handled via high-order operator splitting and a series of one-dimensional, unconditionally stable, semi-Lagrangian updates. To improve both accuracy and efficiency, each plasma species is represented on a different phase-space mesh. On every element and for each species, local positivity-preserving limiters are applied to prevent spurious undershoots and rigorously ensure positivity of the cell averages at the next time step. The accuracy and robustness of the resulting method are verified on several standard single- and two-species test cases, including manufactured solutions, two-stream instability, weak and strong Landau damping, ion-acoustic wave propagation, ion-acoustic shock formation, and KEEN wave dynamics.

In this article, we investigate the stochastic Camassa–Holm equation driven by additive noise, which models nonlinear shallow-water waves under random external forcing. We establish that, almost surely, any finite-time breakdown occurs through wave breaking: the wave amplitude and the positive part of the slope remain bounded, whereas the minimum slope tends to $-\infty$. Under an additional moment assumption, we derive an initial-slope condition ensuring wave breaking with positive probability. We also formulate the equation as a stochastic Hamiltonian PDE, establish its stochastic multi-symplectic conservation law, and derive the evolution law for the averaged $H^1$ energy. Combining the Fourier–Galerkin method, symplectic Runge–Kutta integration, and exact solution of the linear stochastic subsystem, we propose a fully discrete method which preserves an integrated discrete stochastic multi-symplectic conservation law and the linear growth of the discrete averaged \(H^1\) energy. Under suitable regularity and resolution assumptions, we present that the fully discrete dynamics retain the slope-steepening mechanism responsible for wave breaking, providing a discrete counterpart of the continuous behavior.

We analyse a discontinuous Galerkin method for two-dimensional scalar conservation laws on curved domains. Small intersections of a Cartesian grid with the domain are combined by rectangular cell merging. The local space is the restriction of $\mathbb Q_k$ to each physical merged cell. An entropy-conservative flux is supplemented by face viscosity and, for odd degrees, a skew derivative-jump correction on unmerged interior squares. We prove exact semidiscrete $L^2$ dissipation and an optimal error estimate for every fixed $k\ge1$. For smooth solutions with homogeneous trace on the entire boundary, the semi-discrete method satisfies an $O(h^{k+1})$ error bound. The constants are independent of the original cut-cell area fractions. The proof combines a physical-cell $L^2$ projection and cancellation of regular-face projection errors. A separate semidiscrete result covers linear transport with upwind inflow data. Computations for Burgers' equation on a disk and linear advection on a smooth star exhibit orders close to $k+1$ for $k=1,2,3$.

Pair arithmetic omits the renormalization stage performed at the end of each operation of multi-component multiple-precision arithmetic, reducing the operation count and, by eliminating conditional branches, easing SIMD vectorization. Renormalization cannot be omitted entirely, however, because iterative solvers may then fail to converge. We propose RenormBF-pair, a branch-free fixed-trip-count renormalizer parameterized by the word count \(K\) and the input length \(n\): a VecSum sweep, a tail fold, and \(r\) rounds of a FastTwoSum chain, costing \(6(n-1)+(n-K)+3r(K-1)\) flops. Over \(200{,}000\) trials per condition, the smallest round count producing no non-overlap violation was \(r=1\) for \(K \leq 3\) and \(r=2\) for \(K=4\), whereas the same cost spent on a VecSum\(K\)-style TwoSum chain fails for \(K=4\). We delimit with explicit counterexamples both the input condition under which the sum relation holds and the fact that the scheme does not, in general, guarantee a correctly rounded lowest component. Integrating it into CG and BiCGStab over pair arithmetic in binary64 and binary32, we sweep \(240\) paired runs including real matrices screened from SuiteSparse. The proposed scheme and VecSum\(K\) returned bitwise identical words wherever both succeeded, so they are indistinguishable for CG in binary64; differences appear only where VecSum\(K\) breaks, namely binary32 and BiCGStab. There the \(22\)–\(16\) win–loss record is not significant (sign test, \(p=0.42\)); the clear bias is structural, non-overlap violations dropping from \(66\) rows to \(37\). The cost is \(+3.4\%\) for the whole CG solver.

In this paper, an arbitrary Lagrangian-Eulerian (ALE)-based finite element method (FEM) is developed and studied in a monolithic framework for a class of two-phase Stokes flow problems with moving interfaces and jump coefficients, where the mixed finite element approximation to Stokes moving interface problems is established and analyzed in both semi- and fully discrete schemes based on the ALE formulation. The key analytical technique involves a specific H^1-projection associated with the ALE-induced mesh motion due to the evolving interface. Optimal convergence properties of the proposed H^1-projection and its ALE temporal derivative are proved in both H^1 and L^2 norms, with which optimal error estimates are obtained for both semi- and fully discrete mixed finite element approximations to the studied Stokes moving interface problem in both H1 and L2 norms as well. Numerical experiments are carried out to validate all derived theoretical results. The developed analytical approach can be extended to Taylor-Hood and MINI mixed elements, as well as to more general two-phase flow problems.

We introduce Distance-Residual Physics-Informed Neural Networks (DR-PINNs), a physics-informed learning framework for approximating solutions of ordinary and partial differential inclusions (DIs), governing laws in which a differential operator is constrained to lie in a set-valued map rather than equaling a prescribed function. The method replaces the classical pointwise PDE/ODE residual by the squared distance from the differential operator to the admissible set. This distance vanishes exactly when the inclusion is satisfied and measures the infimal correction needed for the operator to enter the admissible set. For a fixed closed convex admissible set, the squared distance is differentiable with respect to the operator value, with gradient given by the metric projection. When the admissible set also depends on the network state, that dependence is included through the chain rule. The framework encompasses ordinary and partial DIs with set-valued reaction terms. For both settings we prove consistency: under the stated closedness, measurability, convexity, and growth assumptions, any sequence of candidates satisfying the initial (and, in the parabolic case, boundary) conditions whose continuous distance-residual functional tends to zero admits a subsequence converging to an exact solution of the target inclusion. These are conditional statements for the continuous distance-residual functional; they do not cover the finite-collocation training loss, the behavior of the optimizer, or convergence rates. For admissible sets given as convex hulls of finitely many vertices, projection onto the set reduces to a small convex quadratic program, making the loss efficiently computable inside the training loop. Numerical experiments demonstrate high accuracy on the differential-inclusion benchmarks.

The classical method of Kaczmarz, introduced in 1937, is a textbook iterative method for solving linear systems $A x = b.$ Despite its widespread use, particularly in the context of solving inverse problems where it is commonly included in software packages within Matlab, Python, and Julia, the precise characterization of convergence has long been deemed difficult to obtain. While different bounds on convergence rates have been established, they are largely unsatisfying as they cannot explain the classical, cyclic-update method's efficient convergence in practice. In this work, we obtain a tight characterization of convergence of the classical Kaczmarz method, establishing both linear and sublinear convergence bounds. The worst-case tight (i.e., exactly attained by some instances in the considered family) bounds are expressed in terms of a fixed matrix that depends only on $A,$ but are not fully interpretable in terms of the matrix spectrum and row correlations, which had been observed to have an impact on convergence. We thus provide relaxations of these bounds that are fully expressible in terms of matrix row norms, row correlations, rank, and extremal positive singular values. The provided relaxed bounds explain one-cycle convergence in special cases where the matrix rows are all either parallel or orthogonal to each other. We further argue that the dependence on different parameters appearing in the bounds is necessary and within a small constant factor of the best attainable in the worst case. Finally, our bounds explain why the classical cyclic update is faster than the randomized one when matrix rows are weakly correlated, which is often observed in inverse problems where the cyclic method is used.

We study finite element approximations for an optimal control problem with pointwise tracking governed by a semilinear elliptic PDE. The control variable enters the state equation as a reaction coefficient, resulting in a bilinear control, and the cost functional includes point evaluations of the state variable. We propose two discretization schemes: a semidiscrete scheme, where the control variable is not discretized, and a fully discrete scheme, where the control variable is discretized using piecewise constant functions. In two and three dimensions, for both schemes, we establish convergence of discrete solutions and prove a priori error bounds that behave as $\mathcal{O}(h|\log h|)$ in the $L^2(\Omega)$-norm for approximating a locally optimal control. In two dimensions, these bounds are improved to $\mathcal{O}(h)$ for the fully discrete scheme and $\mathcal{O}(h^{2}|\log h|^2)$ for the semidiscrete scheme. We conclude with numerical experiments that demonstrate the performance of the proposed schemes.

We study density of states (DOS) estimation for large sparse real symmetric matrices by fitting the cumulative density of states (CDOS) and differentiating the fit. We use a stochastic Lanczos method to supply Ritz values and weights, which we average across Lanczos runs to construct CDOS midpoint data. Monotone piecewise-cubic interpolation of these data yields a smooth CDOS approximation that can be easily differentiated to yield a nonnegative, normalized, piecewise-quadratic DOS without requiring a Gaussian smoothing with a fixed bandwidth in its construction. We also consider approximating CDOS by using Gaussian-process regression (GPR) with single- and double-Gaussian covariance kernels, yielding an approximation with uncertainty information. DOS approximation is then obtained by analytic derivative of the GPR-based CDOS. Experiments on matrices from diverse scientific applications are performed to compare these estimators with Gaussian-broadened stochastic Lanczos and Jackson-damped kernel polynomial approximations. At a common Gaussian validation resolution, we use error measures that assess local discrepancies, integrated errors, and agreement in the overall spectral shape. The results demonstrate that the midpoint spline can produce accurate approximations to the DOS with relatively few matrix–vector products, supporting the use of cumulative interpolation as a practical DOS estimation.

We study a fixed-feature solver, referred to as Spectral-ELM, in which Chebyshev–Gauss–Lobatto (CGL) differentiation matrices are applied to the nodal values of an extreme learning machine (ELM) trial function. Curved domains are treated using a mapped multi-patch formulation. Each curved quadrilateral patch is represented as the image of a reference square, physical derivatives are computed from the associated metric terms, and adjacent patches are coupled by enforcing continuity of the solution and its normal flux. Local QR orthogonalization is used to reduce near-linear dependence among the discrete features. The numerical study includes comparisons with direct CGL collocation, an analytically differentiated ELM, a global Random Feature Method, and a TransNet sampling strategy. It also includes a matched KdV test, a small three-dimensional example, and a Poisson problem on a curved domain containing a hole. For a smooth elliptic problem on a square, direct CGL collocation gives the smallest error. When applied to the same fixed-feature trial space, spectral and analytic differentiation produce comparable errors. For the curved-domain problem, mapped direct CGL collocation and mapped Spectral-ELM achieve similar accuracy with comparable numbers of unknowns. The present formulation is particularly well suited to smooth PDEs in one to three dimensions, where tensor-product grids enable accurate high-order discretizations, while sparse-grid and dimension-adaptive strategies offer promising directions for extension to higher-dimensional problems.

The Koopman operator is a well-established framework for lifting nonlinear dynamical systems to an infinite-dimensional space where dynamics are linear. Extended Dynamic Mode Decomposition (EDMD) is a widely used method in data-driven dynamics as it provides a Galerkin approximation of the Koopman operator on the finite-dimensional span of a prescribed dictionary of observable functions. In particular, EDMD only requires samples from the dynamical system, describing the Koopman operator without knowledge of the underlying dynamics. While many studies analyze asymptotic convergence results for EDMD in terms of data, the question of more practical finite-data convergence results remains incomplete in general settings. In this work, we prove bounds for the finite-data EDMD Galerkin approximation of the Koopman operator by a rate proportional to the inverse of the root of the number of samples (the Monte Carlo rate). These results apply to multiple classes of problems, i.e., for discrete or continuous dynamics, including stochastic systems. In certain more structured setting, we demonstrate super-polynomial rates are achievable both theoretically and computationally. Finally, we derive recovery guarantees for an EDMD variant in the undersampled regime inspired by compressed sensing techniques.

The electrical properties of biological tissues depend on their composition and microstructure and determine their response to applied electric fields. In bone tissue, these properties are closely related to its composition, porosity, and microstructural organization. However, simplified models that link bone microstructure (trabecular and cortical) to its effective electrical properties while accounting for the electrical anisotropy of cortical bone remain scarce. In this work, a multiscale approach was developed, combining finite element simulations with Bruggeman-based homogenization models. Trabecular bone was modeled using microCT derived geometries, and its effective conductivity was analyzed as a function of BV/TV and free water content. Cortical bone was represented by a simplified three-dimensional network of Haversian canals (axial = z) and Volkmann canals (transverse), incorporating their directional organization through an anisotropic Bruggeman formulation. The models captured the conductivity range reported experimentally for trabecular bone (approximately 25-200 mS/m) and reproduced the anisotropic behavior observed in the cortical FEM simulations, with higher conductivity in the axial direction. In the generated cortical geometries, 65.5% of the total porosity was associated with axial Haversian canals. At the macroscale, the tibia model with lower trabecular and cortical BV/TV and reduced cortical thickness showed an approximately 38-39% reduction in impedance relative to the reference model. These results indicate that microstructural variations in bone tissue are reflected in its effective electrical properties and in the macroscopic impedance of the model, providing a physical and computational basis for future investigations into the sensitivity of electrical measurements to changes in bone microstructure.