Papers

Relative periodic orbits (RPOs) are solutions of the three-body problem that are periodic in a uniformly rotating reference frame and, in general, quasi-periodic in inertial coordinates. We present a numerical procedure for computing and continuing one-parameter families of RPOs of the planar Newtonian three-body problem. The method exploits consecutive syzygies, understood here as configurations in which the three bodies are aligned and their velocities satisfy the corresponding symmetry conditions. Matching the positions and momenta at two consecutive alignments reduces the computation of RPOs to a low-dimensional nonlinear problem. Its solutions are then numerically continued, and linear stability is determined from the nontrivial eigenvalues of the rotated monodromy matrix after removing the neutral directions associated with conserved quantities and continuous symmetries. The procedure is applied to several mass distributions and initial configurations, producing families of Poincaré, Hill, and binary-type solutions. These families exhibit transitions from nearly circular to highly eccentric motion, changes of stability near resonances and turning points, and absolute periodic solutions when the rotation angle is a rational multiple of 2{\pi}. In the Hill families, the continuation connects satellite configurations with circumstellar motion as the smallest body loses its gravitational binding to the intermediate body. Circumbinary and circumstellar configurations are also obtained in the binary regime. The results illustrate the dynamical diversity of RPOs and provide coherent three-body motions that can be used as prescribed trajectories in restricted four-body models.

We study optimal sampling recovery in reproducing kernel Hilbert spaces (RKHS) in the uniform norm. For every RKHS with bounded kernel, we establish new comparisons between linear sampling widths and Gelfand widths that overcome the known square-root gap, without requiring a measure or a Christoffel-type condition. Our bounds rely on nested sampling designs obtained by kernel interpolation at (weak) P-greedy points. Under additional (polynomial) decay assumptions the decay rate of the Gelfand widths directly transfers to the sampling widths. With either a logarithmic oversampling or passing to the square root of the Gelfand widths we obtain a direct comparison (requiring no decay assumption) between them. This is particularly effective for super-polynomial decay, such as in Paley-Wiener spaces. Our results follow from representations of both widths in terms of kernel translates and yield, in the opposite direction, a new existence result for a sharp reduced basis selection. Numerical experiments for Legendre, mixed-Sobolev, and Paley-Wiener kernels illustrate our findings.

Reliable epidemic monitoring often requires inferring regional infection burden and transmission heterogeneity from noisy, spatially aggregated, and potentially sparse surveillance data. Agent-based models (ABMs) are attractive for this task because they represent individual behavior, contact heterogeneity, and localized interventions, but these same features make them difficult to calibrate online. We develop a generative AI-based data-assimilation (GenDA) framework for partially observed epidemic ABMs that estimates both the epidemic state and a heterogeneous parameter field while respecting the gap between observable macrostates and latent agent-level microstates. GenDA combines a training-free, score-based generative update for macrostate correction with a direct parameter update based on macrostate discrepancies, followed by a macro-micro reassignment step that restores consistency with the ABM. In controlled and geographically explicit synthetic experiments, the framework recovers regional epidemic burden, dominant hotspot structures, and effective transmission heterogeneity from aggregated observations, while improving post-assimilation forecasts relative to state-only assimilation.

We consider convex optimization with nonlinear inequality constraints and develop a primal-dual multiplier framework that is consistent in continuous and discrete time. We first propose continuous-time dynamics with Nesterov-type vanishing damping $\alpha/t$, together with compatible extrapolations of the dual variable and the nonlinear constraint mapping. Under convexity assumptions and $\alpha\geq3$, we establish $\mathcal O(t^{-2})$ convergence rates for both nonlinear feasibility and the objective residual. In the noncritical regime $\alpha>3$, with an admissible choice of the extrapolation parameter, we further prove that the entire primal-dual trajectory converges to a KKT pair and sharpen both continuous-time estimates to $o(t^{-2})$. We then derive an inexact accelerated primal-dual algorithm through a compatible discretization of a perturbed version of the dynamic. For composite convex objectives, a weighted summability condition on the primal inexactness yields $\mathcal O(k^{-2})$ rates for feasibility and the objective residual. In the corresponding noncritical regime, the discrete primal-dual sequence converges to a KKT pair and both residual estimates improve to $o(k^{-2})$. Thus the continuous and discrete results exhibit matching accelerated rates and matching asymptotic improvements. To the best of our knowledge, this is the first Nesterov-type primal-dual multiplier framework for convex optimization with nonlinear inequality constraints.

We apply Anderson acceleration to the improved Arrow–Hurwicz (IAH) method for the finite element solution of the steady-state incompressible Navier–Stokes equations. The IAH scheme avoids saddle-point solves by decoupling the velocity and pressure updates, but can require prohibitively many iterations, particularly at high Reynolds numbers. To place the acceleration on a rigorous footing, we reformulate the IAH iteration as a nonlinear fixed-point operator $G$ for the grad-div augmented discrete formulation induced by the scheme and establish its well-definedness, Lipschitz continuity, and Fréchet differentiability, thereby verifying the required smoothness conditions locally near the fixed point. Numerical experiments on problems with known analytical solutions, lid-driven cavity flow up to $Re = 15{,}000$, and channel flow over a full step demonstrate that the resulting Anderson-accelerated improved Arrow–Hurwicz algorithm substantially reduces iteration counts and CPU time while retaining the reported manufactured-solution convergence rates and centerline-velocity agreement.

We study the singular limit of a family of linear degenerate interface transmission problems arising from a regularization procedure in the newly proposed Two-Parameter Diffuse Domain Method (DDM2p). For $\alpha>0$, the regularized problem admits a strictly convex variational formulation on $H^{1}(\Omega)$. In the limit $\alpha\to0$, the problem degenerates to a weakly coupled interface system with a nonstandard energy structure. To characterize the limit, we introduce a closed Hilbert subspace $\mathcal{H}\subset H^{1}(\Omega)$, defined through an auxiliary Helmholtz problem on an annular subdomain $\Omega_2\subset \Omega$, and identify the limiting energy functional $\mathcal{E}_{0}$ on $\mathcal{H}$. We prove that the regularized energies $\mathcal{E}_{\alpha}$ $\Gamma$-converge to $\mathcal{E}_{0}$ in the strong $L^{2}(\Omega)$ topology, using the standard framework. Consequently, minimizers of $\mathcal{E}_{\alpha}$ converge to the unique minimizer of $\mathcal{E}_{0}$, which is shown to be equivalent to the solution of the limiting interface problem. We further prove strong convergence $u_{\alpha}\to u_{0}$ in $H^{1}(\Omega)$ and establish an $O(\alpha)$ convergence rate. Numerical experiments in one spatial dimension confirm the predicted first-order convergence rate and suggest that this rate is sharp.

Isogeometric Analysis (IgA) uses the same spline functions to represent the computational domain and to approximate the solution. This allows exact geometry descriptions, but the resulting mass and stiffness matrices are expensive to assemble and to store, especially in three dimensions. We present a projection-based low-rank approach for assembling the mass and stiffness tensors of orientation-preserving tensor-product B-spline geometries. For the mass tensor, we exploit the polynomial structure of the determinant of the Jacobian of the geometry map and represent it in reduced spline product spaces by univariate coefficient transfer operators; with exact quadrature and without truncation, the resulting low-rank tensor is an exact reformulation of the standard Galerkin mass tensor. For the stiffness tensor, we split the rational weight function into a polynomial numerator, again represented in reduced spline product spaces, and the reciprocal determinant, which is in general not a spline function and is therefore approximated by an $L^2$-projection onto a tensor-product spline space. Both constructions are carried out entirely in the tensor-train (TT) format, with the projection system solved by the alternating minimal energy (AMEn) method, so that full high-order coefficient tensors are never formed and the multidimensional integrals reduce to univariate integrals and contracted products. The method is implemented in MATLAB using GeoPDEs and the TT-Toolbox. Numerical experiments show that it is competitive with full assembly and with the interpolation-based low-rank method, and that it applies in two situations in which interpolation is problematic: a singular interpolation system and nearly singular geometries. The construction is restricted to orientation-preserving tensor-product B-spline geometries and does not cover NURBS.

Let A denote the infinite Gram matrix associated with the standard L2 inner product of asymmetrically-weighted (AW) Hermite functions. We derive an explicit representation of its entries and its Cholesky factorization. We further show that this factorization admits a natural interpretation on a scaled Bargmann-Fock basis. An explicit formula for the inverse of A is also obtained. We then consider the corresponding finite Gram matrix and analyze its asymptotic property, as well as that of its Schur complement. The analysis is motivated by numerical methods for plasma physics, in particular Galerkin spectral methods applied to the Vlasov-Poisson (VP) system. As an application, we demonstrate how the derived Gram matrix formulas and asymptotic results can be exploited in the analysis and implementation of a Galerkin spectral method for the VP system.

This work introduces low-complexity, multiplierless approximations for the 32-point discrete Fourier transform. The proposed methods are obtained by minimizing the Frobenius error compared against the DFT matrix over a set of trivial multipliers. A row-wise, symmetry-constrained parameterization is employed to reduce the search space size, rendering the task computationally tractable. The resulting approximations could outperform the reference method in the literature according to energy-based error measurements. A sparse matrix factorization is provided for efficient computation; the arithmetic costs are 152 real additions and 34 bit-shifts only.

A symmetry framework is developed for second-order linear elliptic equations subject to Robin boundary conditions on bounded domains. Orthogonal transformations of the domain are represented through left group actions on scalar, vector, and second-order tensor fields. The corresponding transformation rules for the gradient, divergence, diffusion flux, and conormal boundary term are derived in detail. Based on these relations, symmetry groups are introduced for the principal coefficient tensor, the first-order and zeroth-order coefficients, the differential operator, the Robin coefficient, and the boundary operator. The symmetry properties of the volume and boundary source terms are then incorporated into a common symmetry group for the complete boundary value problem. Under the assumption of unique solvability, it is proved that every element of this common group is also a symmetry of the solution. For reflection symmetries, homogeneous generalized Neumann conditions are obtained on artificial symmetry boundaries, leading to an exact reduction of the computational domain. A corresponding finite element formulation is presented, and numerical results are provided to verify the theoretical symmetry properties and the validity of the resulting domain-reduction strategy. Three examples illustrate radial reduction from multiple dimensions to one dimension, reflection-based domain reduction, and a variable-coefficient problem in which all coefficients and source terms are nonzero.

We develop a new class of inexact block triangular preconditioners for double saddle-point systems arising from PDE-constrained optimization. The proposed preconditioners are constructed through matrix factorization techniques while preserving the inherent block structure of the original systems. A comprehensive spectral analysis of the preconditioned matrices is provided, yielding explicit bounds for both real and nonreal eigenvalues. To enable efficient construction of the inexact preconditioners, randomized strategies are introduced to select the required subblocks. We establish high-probability bounds for the expected approximation error, with the error estimates explicitly characterized in terms of the eigenvalues of the associated matrices. Numerical experiments demonstrate the effectiveness, robustness, and scalability of the proposed preconditioners, and validate the efficiency of the randomized construction strategies.

OptiXDE is a matrix-free spectral operator framework for differential equations on uniform grids and embedded domains. Inspired by angular-spectrum propagation in Fourier optics, it maps transform-diagonal spatial operators to analytical modal multipliers and composes them with physical-space operators for nonlinearities, geometry and boundary enforcement. A common transform–operator–inverse-transform backbone is demonstrated across transient diffusion, periodic and embedded-domain Poisson problems, the cubic nonlinear Schr"odinger equation, viscous Burgers dynamics, the two-dimensional Allen–Cahn equation and incompressible flows from the Taylor–Green vortex to embedded-cylinder vortex shedding. Transform-compatible linear problems are recovered near the floating-point limit, whereas errors on the singular L-shaped domain remain localized near the re-entrant corner and regularized interface. Nonlinear benchmarks recover second-order temporal convergence and the expected conservative or dissipative behavior, while incompressibility remains near round-off level during long-time vortex shedding. The matrix-free updates require \(\mathcal{O}(N\log N)\) work and \(\mathcal{O}(N)\) memory. Device-resident transform workloads reach \(94.9\times\) GPU acceleration, and the complete embedded-cylinder solver achieves a \(42.1\times\) CPU–GPU speedup under matched numerical settings. These results establish OptiXDE as a deterministic and extensible operator-centric framework for structured and embedded-domain differential equations.

The current paper considers the problem of recovering a vector-function on $[-1,1]$ from a limited number of coefficients of its expansion into a series of eigenfunctions of a one-dimensional Dirac system. The Krylov–Lanczos–Eckhoff–Gottlieb acceleration method is examined in the situation when the boundary values it requires have to be computed from the generalized Fourier coefficients themselves. This leads to a $2q\times 2q$ linear system whose matrix is a block Vandermonde matrix; its determinant and inverse are computed explicitly, and the asymptotic $L_2$-error constant of the method is found, paralleling the classical trigonometric case.

In this paper, we propose a novel objective function-free trust-region method designed to solve optimization problems over decentralized networks. Unlike traditional approaches that often rely on stepsize tuning, our framework employs a function-free trust-region procedure that enables adaptive selection of the step length. Our approach accommodates first- and second-order models and eliminates the need to share local function values and gradients among agents, thereby enhancing privacy and computational efficiency. On the theoretical side, we establish provable iteration complexity guarantees that, for some variants, match those established for classical centralized trust-region methods. Numerical evaluations demonstrate that our approach achieves a favorable trade-off between performance and efficiency, requiring only moderate communication overhead compared to state-of-the-art methods in the literature.

We establish sharp mixed spectral Barron regularity for eigenfunctions of molecular Coulomb Hamiltonians. The mixed norm is a Fourier $L^1$ norm with one isotropic weight and coordinate-product weights, and therefore detects regularity invisible to the isotropic Barron scale. For a nonempty set $I$ of electron indices on which the wave function is antisymmetric, we derive an explicit admissible region for the isotropic order $s$ and the coordinate orders $\alpha,\beta$. This region is optimal as a uniform statement over the class of clamped-nuclei Coulomb Hamiltonians. For fixed-spin components with two occupied spin blocks, it reduces to $s+\alpha+\beta<1$; in the fully spin-polarized class it reduces to $s+\alpha<1$. In particular, if $\mathcal I_\sigma$ denotes the family of occupied same-spin blocks determined by $\sigma$, then every fixed-spin spatial component $\psi_\sigma$ satisfies, for every $0\leq\alpha<1$, \[ \left(\sum_{I\in\mathcal I_\sigma}\prod_{i\in I}\langle\xi_i\rangle^\alpha\right)\widehat{\psi_\sigma}\in L^1(\mathbb{R}^{3N}). \] For a fully spin-polarized state, $\mathcal I_\sigma=\{\{1,\ldots,N\}\}$.

We assess computed lift coefficient data for flow past a circular cylinder to evaluate their suitability for practical applications. Specifically, we consider lift coefficient data for a circular cylinder over Reynolds numbers from 120 to 8000. The results are obtained from two-dimensional finite element simulations of the incompressible Navier-Stokes equations using pressure robust discretizations. We compare the computed lift coefficients with published experimental and numerical results, finding good agreement in some cases but significant disagreement in others. Because lift fluctuations are central to vortex-induced vibration concepts, these data therefore provide input for the analysis and preliminary design of bladeless turbines.

The Schottky–Klein prime function is a standard tool for boundary-value problems on multiply connected circular domains. Because this function is represented as an infinite product over a Schottky group, numerical evaluation requires truncation to finitely many factors. The standard word-length truncation grows exponentially in cost and becomes inefficient when the boundary circles nearly touch one another or the unit circle. To address this difficulty, we assign to each group element a cross-ratio potential measuring the size of its contribution, and retain only terms below a prescribed threshold. We establish uniform closed-form bounds on the change in this potential when prepending Schottky-group generators, and from these bounds we derive an efficient enumeration algorithm. The resulting relative error decays exponentially with the threshold at a rate determined by the Hausdorff dimension of the limit set of the Schottky group. Numerical experiments demonstrate that the proposed formulation achieves substantial computational speedups over word-length truncation in challenging geometric configurations.

Functional tissue units (FTUs), including tertiary lymphoid structures (TLSs), blood vessels, and glands, encode localized immune, vascular, and epithelial organization in histopathology. Accurate quantification of these structures is important for studying tissue architecture and disease-associated tissue organization. However, FTUs are frequently sparse, heterogeneous, and surrounded by large amounts of morphologically similar background tissue, making automated segmentation in whole-slide images (WSIs) challenging. We therefore developed FTU-Seek, a pathology foundation model-guided framework that treats morphology-aware negative-patch selection as a key component of sparse FTU segmentation. FTU-Seek uses frozen multi-depth features from the UNI pathology foundation model to train a patch-level classifier that distinguishes FTU-containing from FTU-absent tissue. Target-absent patches are subsequently ranked according to their predicted target-containing probabilities, and the highest-scoring hard negatives are selected through a static Top$K$ strategy to construct compact segmentation training sets. The framework was evaluated using five-fold cross-validation and internal test cohorts across TLS, blood-vessel, and gland segmentation tasks, with an additional independent 30-WSI held-out cohort for TLS. Positive-only, all-tissue, random-negative, and matched random Top$K$ sampling strategies served as comparators. Segmentation-derived phenotypes were further explored in external TCGA cohorts.

In the application of incomplete Cholesky preconditioning to the incomplete Cholesky-conjugate gradient (ICCG) method, forward and backward substitutions exhibit sequential dependencies that constitute a major bottleneck for parallelization in multicore environments. To alleviate this bottleneck, the algebraic block multi-coloring (ABMC) method achieves both parallelism and data locality through block-wise coloring. However, ABMC requires the number of blocks to be specified as an input parameter in advance. This study evaluates the Leiden method as an alternative blocking approach for parallel preconditioning in the ICCG method. As a community detection technique that maximizes a quality function for graph partitioning, the Leiden method automatically generates blocks that reflect the matrix structure without requiring the number of blocks a priori. We partition the adjacency graphs of sparse matrices using the Leiden method and utilize the resulting blocks for parallel preconditioning. We implement the Leiden method using modularity and the constant Potts model as quality functions and compare its performance with that of the ABMC method in terms of the number of iterations, execution time, and L2 cache efficiency across eight symmetric positive definite matrices. The experimental results demonstrate that the Leiden method with the constant Potts model achieves performance comparable to that of the ABMC method configured with an optimized number of blocks.

A bottleneck in the theory of kernel methods in machine learning is the storage requirement. To ameliorate this, a standard trick is to replace the kernel with an explicit feature map. Perhaps, the most well known example is the Gaussian kernel which can be expressed in terms of the Fourier features. Analytically, the kernel $K$ can be expressed in terms of an integral expression that involves a possibly asymmetric kernel $G$ representing the feature map. Numerically, one needs to approximate this integral by a suitable numerical integration scheme, typically Monte Carlo. In this paper, we demonstrate that the eigenvalues of $K$ are approximated much better by the eigenvalues of the kernel obtained by discretizing the integral using suitable quadrature formulas instead. We illustrate this fact in the case of the Gaussian kernel and neural tangent kernels on the unit sphere of a four-dimensional Euclidean space corresponding to the sigmoid and ReLU activation functions.

The complexity of time-domain simulation of modern power systems has increased significantly because converter-based resources introduce control dynamics that must be simulated alongside slower system-level and fast electromagnetic dynamics. The resulting wide range of timescales may force classical time-domain solvers to use small timesteps, complicate the solution of nonlinear equations at each timestep, and reduce solver reliability under strongly nonlinear and multi-timescale transient conditions. This paper introduces RAPTOR, a first-of-its-kind time-domain simulation framework for power system dynamic simulations. At its core, RAPTOR introduces a new integration technique that represents the unknown trajectory of hybrid differential-algebraic equations (DAEs) over a time interval using a physics-informed random-projection neural network (PIRPNN) built from fixed Gaussian radial basis functions (RBFs). A nonlinear solver, e.g., Newton-Raphson, determines how to combine the fixed RBFs so that the resulting trajectory satisfies the governing equations. By combining RBFs with broad and localized shapes, RAPTOR can represent complex multi-timescale behavior over comparatively long time intervals. This enables larger effective simulation advances, fewer overall nonlinear solver iterations, and faster integration over long time horizons. Numerical studies on stiff RMS models, IEEE 9-, 14-, 39-, 57-, and 118-bus RMS simulation benchmarks, and EMT test cases show that RAPTOR accurately finds solutions with substantially fewer sequential simulation advances while offering superior accuracy-runtime performance. With speedups of over 10x in certain cases, these results showcase RAPTOR's potential to outperform long-standing and widely used integration methods, such as the Radau method and the trapezoidal rule.

Time-harmonic Maxwell scattering by a perfectly conducting obstacle is governed by the electric field integral equation (EFIE), whose natural energy space is $\boldsymbol H^{-1/2}(\mathrm{div}_\Gamma,\Gamma)$. For countably parametrized surface deformations, both the operator and its energy space depend on the geometry. We establish operator-valued shape holomorphy of the EFIE operator on a fixed reference space and derive dimension-independent sparse approximation rates. The main analytical difficulty is a fractional mapping property absent from existing $L^2$-based operator-valued shape-holomorphy theory for weakly singular kernels: the Maxwell graph norm requires uniform holomorphy of the complex-deformed scalar single layer as an operator from $H^{-1/2}$ to $H^{1/2}$. We prove this one-order smoothing by realizing the Laplace principal part as the trace of a uniformly sectorial complex-coefficient divergence-form problem on a fixed real ambient domain. Combined with a surface contravariant Piola transformation, which identifies the geometry-dependent Maxwell trace spaces and cancels the surface Jacobians in the pulled-back EFIE exactly, this yields $(\mathbf b,p,\varepsilon)$-holomorphy of the EFIE operator family for $\ell^p$-summable shape deformations, $0<p<1$. Under pointwise exclusion of interior electric resonances, parameter-uniform invertibility is derived rather than assumed, and the surface current and far field inherit the same parametric regularity. Their Legendre coefficients are $\ell^p$-summable, yielding dimension-independent best $N$-term approximation rates. The operator-level results also transfer to Galerkin boundary element operators on a single reference mesh, with constants independent of the discretization dimension, and provide reusable surrogates for multiple incident fields and bounded linear observables.

Unlike the classical Heston model, Fourier pricing under the rough Heston model requires solving a fractional Riccati equation at every quadrature point. Since the required resolution varies with model parameters and quadrature point, a single uniform time discretization can be inefficient. We develop single- and multilevel Gauss-Laguerre quadrature methods that balance the time discretization and Fourier quadrature errors. Both methods scale the laguerre weight to the estimated Fourier integrand decay. The single-level method allocates a prescribed tolerance between the two errors. The multilevel method splits the integrand into a level-zero term and level differences, selecting quadrature points separately at each level. Suppose that the Fourier integrand discretization error is $O(\Delta t^p)$, that evaluating the characteristic function once costs $O(\Delta t^{-\beta})$, and that the algebraic Gauss-Laguerre quadrature error is $O(N^{-s_{SL}/2})$, where $s_{{SL}}$ is the smoothness index. Under this estimate and assumptions on the regularity and decay of level differences, we prove that the proposed single-level method requires $O(\epsilon^{-(\beta/p+2/s_{{SL}})})$ computational work to achieve accuracy $\epsilon$, whereas the proposed multilevel method requires $O(\epsilon^{-\beta/p})$ computational work. We also study root-exponential Gauss-Laguerre error models for practical multilevel quadrature allocation. Numerical experiments support the observed fractional Riccati and Fourier integrand convergence rates and root-exponential quadrature behavior, and show substantial reductions in quadrature cost from the proposed scaling. The multilevel method provides clear computational savings over the single-level method. We further benchmark the multilevel fractional Riccati method against the BL2 Markovian approximation and report lower total CPU time in the tested configurations.

Generative models for implied volatility surfaces must produce outputs that satisfy static no-arbitrage constraints. We study these constraints in latent space. For a fixed generator, we assign each latent code a scalar margin determined by the no-arbitrage conditions of the generated surface. The codes with nonnegative margin form the admissible latent set. We establish conditions under which strictly admissible codes remain admissible under small perturbations and the boundary of the admissible set is characterized by zero margin. For regular boundary components, we formulate a level-set equation whose local dynamics are directed toward the zero-margin set. The analysis treats the generator as a map from latent variables to surfaces and is therefore not restricted to a particular architecture. It applies to variational autoencoders, generative adversarial networks, and other generative models with a deterministic realization map. Numerical tests recover known boundaries in analytic examples. Experiments with a variational autoencoder trained on Heston surfaces show that similar reconstruction errors can correspond to different admissible regions and that the latent prior may be concentrated inside such a region. The computed boundary can also be used to modify latent codes that generate violating surfaces.

Metropolitan regions are simultaneously pursuing several interventions to improve mobility, accessibility, and energy efficiency, necessitating integrated tools to understand how these policies interact to affect travel behavior, energy use, and infrastructure needs. This paper evaluates the combined impacts of electrification, freight demand management, road pricing, parking reform, and transit expansion on the Chicago metropolitan transportation system in 2050, using a business-as-usual (BAU) scenario as the baseline. We employ POLARIS, a large-scale agent-based modeling framework calibrated to 2019 conditions, to simulate nine policy scenarios for the seven-county northeastern Illinois region. The framework co-simulates activity-based passenger demand, endogenous freight generation, multimodal traffic assignment, and transit operations, with charging infrastructure and freight operations optimized for each case. Our findings reveal that under the high electrification scenario, total fuel mass declines by 68% while total charging energy increases by approximately 4-8x from BAU, resulting in a peak power demand near 4 GW concentrated in the urban core. Furthermore, freight management policies reduce freight VMT by increasing trip frequency but shortening distances, smart road pricing most effectively reduces auto VMT, and transit expansion boosts ridership by 18% relative to BAU. By presenting the first integrated, agent-based scenario framework for Chicago that jointly evaluates these interventions, this study provides actionable insights for regional transportation planning, grid infrastructure investment, and emissions reduction, highlighting the value of targeted charger upgrades and coordinated policy bundles.

We present a meshfree numerical solver for the incompressible Navier-Stokes equations on oriented curved surfaces that are represented by surface point clouds. On curved surfaces, numerical challenges pertaining to stiffness and pressure-velocity coupling are exacerbated. Moreover, vector calculus on curved surfaces differs from its Euclidean counterpart. The presented method operates on surface point clouds in an Eulerian frame of reference without requiring a computational grid or mesh. It achieves consistent approximation in space and time with high order of accuracy; we demonstrate up to order six. The incompressibility constraint is locally imposed as a weak artificial compressibility approximation, avoiding global matrix inversion. We show that the method provides consistent and convergent approximations of surface vector fields and differential operators. We study the relationship between error, spatial resolution, and artificial Mach number and characterize the frequency spectrum of the artificial oscillations. We provide numerical solutions of the incompressible Navier-Stokes equations on symmetric surfaces, such as the sphere and torus, and on parametric and non-parametric asymmetric surfaces. Since the proposed method works directly on unstructured surface point clouds, it provides a promising approach for simulations on image-derived geometries, such as in biological morphogenesis from microscopy videos.

We consider a fractional diffusion problem with temporal nonlocality acting on the diffusive flux. A coercive space–time variational formulation in Bochner-valued fractional Sobolev spaces is derived and the existence, uniqueness, and regularity of solutions are established. We further develop a conforming tensor-product Galerkin discretization and prove quasi-optimal error estimates in the anisotropic energy norm and improved convergence rates in weaker norms using duality arguments. In contrast to some space-time formulations for classical diffusion, the method preserves the causal structure of the evolution problem and leads to a time-stepping procedure with memory terms. On uniform time grids, the discrete history operator has a lower-triangular Toeplitz structure which enables an efficient implementation using fast recursive convolution techniques.

Quantum signal processing (QSP) provides a simple and efficient framework for implementing polynomial transformations using quantum circuits. Its classical design stage leads to a constrained minimax approximation problem: find a polynomial of prescribed parity that approximates a target function uniformly on a fitting set while remaining bounded in magnitude by one on the domain $[0,1]$, which can be viewed as a semi-infinite constraint. Discretization converts the problem into a linear program, but feasibility at a set of finitely many sampled points does not ensure feasibility on the whole domain, especially when an optimal approximant reaches the boundary of the feasible set. We investigate two approaches to address this difficulty. A Remez exchange method combined with active-set constraint enforcement is efficient on many tested instances, but its stability depends on the target and problem geometry. We then introduce nonlinear Fourier retraction, which uses QSP completion and phase synthesis to turn a nearly feasible polynomial into phase factors for a feasible QSP polynomial without increasing the degree. Across representative problems, retraction largely preserves approximation accuracy and remains effective on instances where the Remez heuristic is unstable. The resulting workflow connects classical minimax approximation and semi-infinite optimization with nonlinear Fourier analysis, and is implemented in the qsppack software package.

Markovian open quantum systems are well modeled by the Lindblad Master Equation (ME) $\frac{\mathrm{d}}{\mathrm{d} t} \rho_t = \mathcal{L} \rho_t$, where $\mathcal{L}$ is a linear (super-)operator and $\rho_t$ is the system state, a positive matrix. When designing or characterizing a quantum system, one is usually interested in the steady state $\rho_\infty$ (such that $\mathcal{L} \rho_\infty = 0$), the first few excited states, and trajectories $t\mapsto \rho_t$. In finite dimension, $\rho_t$ is an $n\times n$ matrix, $\mathcal{L}$ thus typically costs $n^4$ to store explicitly as a dense matrix, and $O(n^6)$ to diagonalize or invert exactly, making standard linear algebraic techniques expensive for large systems. However, $\mathcal{L}$ usually costs only $O(n^3)$ to apply. This makes iterative methods appealing, but they do not work without a good preconditioner. In this article, our main observation is that a part of the Lindblad equation, corresponding to the so-called no-jump evolution $\mathcal{S}$, can be inverted efficiently. Using this inverse map, we introduce an auxiliary completely positive trace-preserving (CPTP) map $\Phi$ whose fixed point is directly related to $\rho_\infty$, all the other eigenvalues having smaller magnitude. The map $\Phi$ is thus well suited to iterative methods, and $\rho_\infty$ can be found in a few Arnoldi iterations. Using the same inverse map $\mathcal{S}^{-1}$ as preconditioner, we compute the low-lying spectrum efficiently via shift-invert Arnoldi, and, as a proof of concept, build an implicit time integrator that is competitive on stiff systems in the low-precision regime. For the steady-state and low excited states problems, our methods scale like $O(n^3)$ per iteration and offer state-of-the-art performance on CPU and GPU.

In this paper, we consider structure-preserving model reduction of port-Hamiltonian (pH) systems which extend classical Hamiltonian systems with dissipation and an input-output port. These pH systems are often used in multi-physics systems, as the interconnection of one or more \pH systems results again in a pH system. If particularly the system matrices associated with the interconnection and/or dissipation of a pH system are state-dependent, then the evaluation of standard reduced-order models (ROMs) may depend on the dimension of the original full-order model, resulting in high computational costs. To circumvent these high costs, we propose to use structure-preserving neural networks. In particular, we perform two steps: (1) we use the generalized manifold Galerkin projection to project the pH system onto the reduced space; then (2) we train a neural network to learn the map from the reduced-order state to the reduced-order interconnection and dissipation system matrices. To ensure that the resulting ROM is again a pH system, the architecture of the neural network is chosen such that the skew-symmetry and positive semi-definiteness of the reduced-order systems matrices are maintained. In a numerical example, we consider a nonlinear mass-spring-damper system with state-dependent system matrices. The numerical results show that the proposed method achieves a significant computational speed-up compared to the original \ROM with comparable accuracy.