Numerical Analysis
See recent articles
Showing new listings for Friday, 18 April 2025
- [1] arXiv:2504.12567 [pdf, html, other]
-
Title: The existence of explicit symplectic integrators for general nonseparable Hamiltonian systemsSubjects: Numerical Analysis (math.NA)
The existence of explicit symplectic integrators for general nonseparable Hamiltonian systems is an open and important problem in both numerical analysis and computing in science and engineering, as explicit integrators are usually more efficient than the implicit integrators of the same order of accuracy. Up to now, all responses to this problem are negative. That is, there exist explicit symplectic integrators only for some special nonseparable Hamiltonian systems, whereas the universal design involving explicit symplectic integrators for general nonseparable Hamiltonian systems has not yet been studied sufficiently. In this paper, we present a constructive proof for the existence of explicit symplectic integrators for general nonseparable Hamiltonian systems via finding explicit symplectic mappings under which the special submanifold of the extended phase space is invariant. It turns out that the proposed explicit integrators are symplectic in both the extended phase space and the original phase space. Moreover, on the basis of the global modified Hamiltonians of the proposed integrators, the backward error analysis is made via a parameter relaxation and restriction technique to show the linear growth of global errors and the near-preservation of first integrals. In particular, the effective estimated time interval is nearly the same as classical implicit symplectic integrators when applied to (near-) integrable Hamiltonian systems. Numerical experiments with a completely integrable nonseparable Hamiltonian and a nonintegrable nonseparable Hamiltonian illustrate the good long-term behavior and high efficiency of the explicit symplectic integrators proposed and analyzed in this paper.
- [2] arXiv:2504.12631 [pdf, html, other]
-
Title: Geometry-preserving Numerical Scheme for Riemannian Stochastic Differential EquationsSubjects: Numerical Analysis (math.NA)
Stochastic differential equations (SDEs) on Riemannian manifolds have numerous applications in system identification and control. However, geometry-preserving numerical methods for simulating Riemannian SDEs remain relatively underdeveloped. In this paper, we propose the Exponential Euler-Maruyama (Exp-EM) scheme for approximating solutions of SDEs on Riemannian manifolds. The Exp-EM scheme is both geometry-preserving and computationally tractable. We establish a strong convergence rate of $\mathcal{O}(\delta^{\frac{1 - \epsilon}{2}})$ for the Exp-EM scheme, which extends previous results obtained for specific manifolds to a more general setting. Numerical simulations are provided to illustrate our theoretical findings.
- [3] arXiv:2504.12650 [pdf, html, other]
-
Title: Tangent Space Parametrization for Stochastic Differential Equations on SO(n)Subjects: Numerical Analysis (math.NA)
In this paper, we study the numerical simulation of stochastic differential equations (SDEs) on the special orthogonal Lie group $\text{SO}(n)$. We propose a geometry-preserving numerical scheme based on the stochastic tangent space parametrization (S-TaSP) method for state-dependent multiplicative SDEs on $\text{SO}(n)$. The convergence analysis of the S-TaSP scheme establishes a strong convergence order of $\mathcal{O}(\delta^{\frac{1-\epsilon}{2}})$, which matches the convergence order of the previous stochastic Lie Euler-Maruyama scheme while avoiding the computational cost of the exponential map. Numerical simulation illustrates the theoretical results.
- [4] arXiv:2504.12713 [pdf, html, other]
-
Title: Efficient Primal-dual Forward-backward Splitting Method for Wasserstein-like Gradient Flows with General Nonlinear MobilitiesComments: 47pages, 12 figuresSubjects: Numerical Analysis (math.NA); Optimization and Control (math.OC)
We construct an efficient primal-dual forward-backward (PDFB) splitting method for computing a class of minimizing movement schemes with nonlinear mobility transport distances, and apply it to computing Wasserstein-like gradient flows. This approach introduces a novel saddle point formulation for the minimizing movement schemes, leveraging a support function form from the Benamou-Brenier dynamical formulation of optimal transport. The resulting framework allows for flexible computation of Wasserstein-like gradient flows by solving the corresponding saddle point problem at the fully discrete level, and can be easily extended to handle general nonlinear mobilities. We also provide a detailed convergence analysis of the PDFB splitting method, along with practical remarks on its implementation and application. The effectiveness of the method is demonstrated through several challenging numerical examples.
- [5] arXiv:2504.12892 [pdf, html, other]
-
Title: Manifold-valued function approximation from multiple tangent spacesComments: 25 pages, 7 figuresSubjects: Numerical Analysis (math.NA)
Approximating a manifold-valued function from samples of input-output pairs consists of modeling the relationship between an input from a vector space and an output on a Riemannian manifold. We propose a function approximation method that leverages and unifies two prior techniques: (i) approximating a pullback to the tangent space, and (ii) the Riemannian moving least squares method. The core idea of the new scheme is to combine pullbacks to multiple tangent spaces with a weighted Fréchet mean. The effectiveness of this approach is illustrated with numerical experiments on model problems from parametric model order reduction.
- [6] arXiv:2504.12938 [pdf, html, other]
-
Title: Optimal analysis of penalized lowest-order mixed FEMs for the Stokes-Darcy modelSubjects: Numerical Analysis (math.NA)
This paper is concerned with non-uniform fully-mixed FEMs for dynamic coupled Stokes-Darcy model with the well-known Beavers-Joseph-Saffman (BJS) interface condition. In particular, a decoupled algorithm with the lowest-order mixed non-uniform FE approximations (MINI for the Stokes equation and RT0-DG0 for the Darcy equation) and the classical Nitsche-type penalty is studied. The method with the combined approximation of different orders is commonly used in practical simulations. However, the optimal error analysis of methods with non-uniform approximations for the coupled Stokes-Darcy flow model has remained challenging, although the analysis for uniform approximations has been well done. The key question is how the lower-order approximation to the Darcy flow influences the accuracy of the Stokes solution through the interface condition. In this paper, we prove that the decoupled algorithm provides a truly optimal convergence rate in L^2-norm in spatial direction: O(h^2) for Stokes velocity and O(h) for Darcy flow in the coupled Stokes-Darcy model. This implies that the lower-order approximation to the Darcy flow does not pollute the accuracy of numerical velocity for Stokes flow. The analysis presented in this paper is based on a well-designed Stokes-Darcy Ritz projection and given for a dynamic coupled model. The optimal error estimate holds for more general combined approximations and more general coupled models, including the corresponding model of steady-state Stokes-Darcy flows and the model of coupled dynamic Stokes and steady-state Darcy flows. Numerical results confirm our theoretical analysis and show that the decoupled algorithm is efficient.
- [7] arXiv:2504.13036 [pdf, other]
-
Title: A generalized energy-based modeling framework with application to field/circuit coupled problemsSubjects: Numerical Analysis (math.NA)
This paper presents a generalized energy-based modeling framework extending recent formulations tailored for differential-algebraic equations. The proposed structure, inspired by the port-Hamiltonian formalism, ensures passivity, preserves the power balance, and facilitates the consistent interconnection of subsystems. A particular focus is put on low-frequency power applications in electrical engineering. Stranded, solid, and foil conductor models are investigated in the context of the eddy current problem. Each conductor model is shown to fit into the generalized energy-based structure, which allows their structure-preserving coupling with electrical circuits described by modified nodal analysis. Theoretical developments are validated through a numerical simulation of an oscillator circuit, demonstrating energy conservation in lossless scenarios and controlled dissipation when eddy currents are present.
New submissions (showing 7 of 7 entries)
- [8] arXiv:2504.12836 (cross-list from math.AP) [pdf, html, other]
-
Title: Inverse iteration method for higher eigenvalues of the $p$-LaplacianComments: 29 pages, 5 figuresSubjects: Analysis of PDEs (math.AP); Numerical Analysis (math.NA); Spectral Theory (math.SP)
We propose a characterization of a $p$-Laplace higher eigenvalue based on the inverse iteration method with balancing the Rayleigh quotients of the positive and negative parts of solutions to consecutive $p$-Poisson equations. The approach relies on the second eigenvalue's minimax properties, but the actual limiting eigenvalue depends on the choice of initial function. The well-posedness and convergence of the iterative scheme are proved. Moreover, we provide the corresponding numerical computations. As auxiliary results, which also have an independent interest, we provide several properties of certain $p$-Poisson problems.
- [9] arXiv:2504.12949 (cross-list from cs.LG) [pdf, html, other]
-
Title: RL-PINNs: Reinforcement Learning-Driven Adaptive Sampling for Efficient Training of PINNsSubjects: Machine Learning (cs.LG); Numerical Analysis (math.NA)
Physics-Informed Neural Networks (PINNs) have emerged as a powerful framework for solving partial differential equations (PDEs). However, their performance heavily relies on the strategy used to select training points. Conventional adaptive sampling methods, such as residual-based refinement, often require multi-round sampling and repeated retraining of PINNs, leading to computational inefficiency due to redundant points and costly gradient computations-particularly in high-dimensional or high-order derivative scenarios. To address these limitations, we propose RL-PINNs, a reinforcement learning(RL)-driven adaptive sampling framework that enables efficient training with only a single round of sampling. Our approach formulates adaptive sampling as a Markov decision process, where an RL agent dynamically selects optimal training points by maximizing a long-term utility metric. Critically, we replace gradient-dependent residual metrics with a computationally efficient function variation as the reward signal, eliminating the overhead of derivative calculations. Furthermore, we employ a delayed reward mechanism to prioritize long-term training stability over short-term gains. Extensive experiments across diverse PDE benchmarks, including low-regular, nonlinear, high-dimensional, and high-order problems, demonstrate that RL-PINNs significantly outperforms existing residual-driven adaptive methods in accuracy. Notably, RL-PINNs achieve this with negligible sampling overhead, making them scalable to high-dimensional and high-order problems.
- [10] arXiv:2504.12981 (cross-list from physics.med-ph) [pdf, html, other]
-
Title: Efficient Chebyshev Reconstruction for the Anisotropic Equilibrium Model in Magnetic Particle ImagingComments: This work has been submitted to the IEEE for possible publicationSubjects: Medical Physics (physics.med-ph); Image and Video Processing (eess.IV); Numerical Analysis (math.NA)
Magnetic Particle Imaging (MPI) is a tomographic imaging modality capable of real-time, high-sensitivity mapping of superparamagnetic iron oxide nanoparticles. Model-based image reconstruction provides an alternative to conventional methods that rely on a measured system matrix, eliminating the need for laborious calibration measurements. Nevertheless, model-based approaches must account for the complexities of the imaging chain to maintain high image quality. A recently proposed direct reconstruction method leverages weighted Chebyshev polynomials in the frequency domain, removing the need for a simulated system matrix. However, the underlying model neglects key physical effects, such as nanoparticle anisotropy, leading to distortions in reconstructed images. To mitigate these artifacts, an adapted direct Chebyshev reconstruction (DCR) method incorporates a spatially variant deconvolution step, significantly improving reconstruction accuracy at the cost of increased computational demands. In this work, we evaluate the adapted DCR on six experimental phantoms, demonstrating enhanced reconstruction quality in real measurements and achieving image fidelity comparable to or exceeding that of simulated system matrix reconstruction. Furthermore, we introduce an efficient approximation for the spatially variable deconvolution, reducing both runtime and memory consumption while maintaining accuracy. This method achieves computational complexity of O(N log N ), making it particularly beneficial for high-resolution and three-dimensional imaging. Our results highlight the potential of the adapted DCR approach for improving model-based MPI reconstruction in practical applications.
Cross submissions (showing 3 of 3 entries)
- [11] arXiv:2311.01451 (replaced) [pdf, other]
-
Title: Randomized Strong Recursive Skeletonization: Simultaneous Compression and LU Factorization of Hierarchical Matrices using Matrix-Vector ProductsSubjects: Numerical Analysis (math.NA)
The hierarchical matrix framework partitions matrices into subblocks that are either small or of low numerical rank, enabling linear storage complexity and efficient matrix-vector multiplication. This work focuses on the $\mathcal{H}^2$-matrix format, whose defining feature is the nested basis property which allows basis matrices to be reused across different levels of the hierarchy. While $\mathcal{H}^2$-matrices support fast Cholesky and LU factorizations, implementing these methods is challenging -- especially for 3D PDE discretizations -- due to the complexity of nested recursions and recompressions. Moreover, compressing $\mathcal{H}^2$-matrices becomes particularly difficult when only matrix-vector multiplication operations are available.
This paper introduces an algorithm that simultaneously compresses and factorizes a general $\mathcal{H}^{2}$-matrix, using only the action of the matrix and its adjoint on vectors. The number of required matrix-vector products is independent of the matrix size, depending only on the problem geometry and a rank parameter that captures low-rank interactions between well-separated boxes. The resulting LU factorization is invertible and can serve as an approximate direct solver, with its accuracy influenced by the spectral properties of the matrix.
To achieve competitive sample complexity, the method uses dense Gaussian test matrices without explicitly encoding structured sparsity patterns. Samples are drawn only once at the start of the algorithm; as the factorization proceeds, structure is dynamically introduced into the test matrices through efficient linear algebraic operations. Numerical experiments demonstrate the algorithm's robustness to indefiniteness and ill-conditioning, as well as its efficiency in terms of sample cost for challenging problems arising from both integral and differential equations in 2D and 3D. - [12] arXiv:2404.16324 (replaced) [pdf, html, other]
-
Title: Improved impedance inversion by the iterated graph LaplacianSubjects: Numerical Analysis (math.NA); Machine Learning (cs.LG); Signal Processing (eess.SP)
We introduce a data-adaptive inversion method that integrates classical or deep learning-based approaches with iterative graph Laplacian regularization, specifically targeting acoustic impedance inversion - a critical task in seismic exploration. Our method initiates from an impedance estimate derived using either traditional inversion techniques or neural network-based methods. This initial estimate guides the construction of a graph Laplacian operator, effectively capturing structural characteristics of the impedance profile. Utilizing a Tikhonov-inspired variational framework with this graph-informed prior, our approach iteratively updates and refines the impedance estimate while continuously recalibrating the graph Laplacian. This iterative refinement shows rapid convergence, increased accuracy, and enhanced robustness to noise compared to initial reconstructions alone. Extensive validation performed on synthetic and real seismic datasets across varying noise levels confirms the effectiveness of our method. Performance evaluations include four initial inversion methods: two classical techniques and two neural networks - previously established in the literature.
- [13] arXiv:2410.21817 (replaced) [pdf, html, other]
-
Title: Backward error analysis of stochastic Poisson integratorsSubjects: Numerical Analysis (math.NA)
We address our attention to the numerical time discretization of stochastic Poisson systems via Poisson integrators. The aim of the investigation regards the backward error analysis of such integrators to reveal their ability of being structure-preserving, for long times of integration. In particular, we first provide stochastic modified equations suitable for such integrators and then we rigorously study them to prove accurate estimates on the long-term numerical error along the dynamics generated by stochastic Poisson integrators, with reference to the preservation of the random Hamiltonian conserved along the exact flow of the approximating Wong-Zakai Poisson system. Finally, selected numerical experiments confirm the effectiveness of the theoretical analysis.
- [14] arXiv:2503.19185 (replaced) [pdf, html, other]
-
Title: Least Squares with Equality constraints Extreme Learning Machines for the resolution of PDEsSubjects: Numerical Analysis (math.NA)
In this paper, we investigate the use of single hidden-layer neural networks as a family of ansatz functions for the resolution of partial differential equations (PDEs). In particular, we train the network via Extreme Learning Machines (ELMs) on the residual of the equation collocated on -- eventually randomly chosen -- points. Because the approximation is done directly in the formulation, such a method falls into the framework of Physically Informed Neural Networks (PINNs) and has been named PIELM. Since its first introduction, the method has been refined variously, and one successful variant is the Extreme Theory of Functional Connections (XTFC). However, XTFC strongly takes advantage of the description of the domain as a tensor product. Our aim is to extend XTFC to domains with general shapes. The novelty of the procedure proposed in the present paper is related to the treatment of boundary conditions via constrained imposition, so that our method is named Least Squares with Equality constraints ELM (LSEELM). An in-depth analysis and comparison with the cited methods is performed, again with the analysis of the convergence of the method in various scenarios. We show the efficiency of the procedure both in terms of computational cost and in terms of overall accuracy.
- [15] arXiv:2504.06998 (replaced) [pdf, html, other]
-
Title: A Krylov projection algorithm for large symmetric matrices with dense spectraComments: Block Lanczos, Quadrature, Transfer function, Kreĭn-Nudelman, Hermite-PadéSubjects: Numerical Analysis (math.NA)
We consider the approximation of $B^T (A+sI)^{-1} B$ for large s.p.d. $A\in\mathbb{R}^{n\times n}$ with dense spectrum and $B\in\mathbb{R}^{n\times p}$, $p\ll n$. We target the computations of Multiple-Input Multiple-Output (MIMO) transfer functions for large-scale discretizations of problems with continuous spectral measures, such as linear time-invariant (LTI) PDEs on unbounded domains. Traditional Krylov methods, such as the Lanczos or CG algorithm, are known to be optimal for the computation of $(A+sI)^{-1}B$ with real positive $s$, resulting in an adaptation to the distinctively discrete and nonuniform spectra. However, the adaptation is damped for matrices with dense spectra. It was demonstrated in [Zimmerling, Druskin, Simoncini, Journal of Scientific Computing 103(1), 5 (2025)] that averaging Gauß and Gauß-Radau quadratures computed using the block-Lanczos method significantly reduces approximation errors for such problems. Here, we introduce an adaptive Kreĭn-Nudelman extension to the (block) Lanczos recursions, allowing further acceleration at negligible $o(n)$ cost. Similar to the Gauß-Radau quadrature, a low-rank modification is applied to the (block) Lanczos matrix. However, unlike the Gauß-Radau quadrature, this modification depends on $\sqrt{s}$ and can be considered in the framework of the Hermite-Padé approximants, which are known to be efficient for problems with branch-cuts, that can be good approximations to dense spectral intervals. Numerical results for large-scale discretizations of heat-diffusion and quasi-magnetostatic Maxwell's operators in unbounded domains confirm the efficiency of the proposed approach.
- [16] arXiv:2504.11968 (replaced) [pdf, other]
-
Title: Dynamical reweighting for estimation of fluctuation formulasSubjects: Numerical Analysis (math.NA)
We propose a variance reduction method for calculating transport coefficients in molecular dynamics using an importance sampling method via Girsanov's theorem applied to Green--Kubo's formula. We optimize the magnitude of the perturbation applied to the reference dynamics by means of a scalar parameter~$\alpha$ and propose an asymptotic analysis to fully characterize the long-time behavior in order to evaluate the possible variance reduction. Theoretical results corroborated by numerical results show that this method allows for some reduction in variance, although rather modest in most situations.
- [17] arXiv:2410.03802 (replaced) [pdf, html, other]
-
Title: Mesh-Informed Reduced Order Models for Aneurysm Rupture Risk PredictionGiuseppe Alessio D'Inverno, Saeid Moradizadeh, Sajad Salavatidezfouli, Pasquale Claudio Africa, Gianluigi RozzaSubjects: Medical Physics (physics.med-ph); Machine Learning (cs.LG); Numerical Analysis (math.NA)
The complexity of the cardiovascular system needs to be accurately reproduced in order to promptly acknowledge health conditions; to this aim, advanced multifidelity and multiphysics numerical models are crucial. On one side, Full Order Models (FOMs) deliver accurate hemodynamic assessments, but their high computational demands hinder their real-time clinical application. In contrast, Reduced Order Models (ROMs) provide more efficient yet accurate solutions, essential for personalized healthcare and timely clinical decision-making. In this work, we explore the application of computational fluid dynamics (CFD) in cardiovascular medicine by integrating FOMs with ROMs for predicting the risk of aortic aneurysm growth and rupture. Wall Shear Stress (WSS) and the Oscillatory Shear Index (OSI), sampled at different growth stages of the thoracic aortic aneurysm, are predicted by means of Graph Neural Networks (GNNs). GNNs exploit the natural graph structure of the mesh obtained by the Finite Volume (FV) discretization, taking into account the spatial local information, regardless of the dimension of the input graph. Our experimental validation framework yields promising results, confirming our method as a valid alternative that overcomes the curse of dimensionality.
- [18] arXiv:2503.04649 (replaced) [pdf, html, other]
-
Title: Transferable Foundation Models for Geometric Tasks on Point Cloud Representations: Geometric Neural OperatorsSubjects: Machine Learning (cs.LG); Computer Vision and Pattern Recognition (cs.CV); Numerical Analysis (math.NA); Optimization and Control (math.OC)
We introduce methods for obtaining pretrained Geometric Neural Operators (GNPs) that can serve as basal foundation models for use in obtaining geometric features. These can be used within data processing pipelines for machine learning tasks and numerical methods. We show how our GNPs can be trained to learn robust latent representations for the differential geometry of point-clouds to provide estimates of metric, curvature, and other shape-related features. We demonstrate how our pre-trained GNPs can be used (i) to estimate the geometric properties of surfaces of arbitrary shape and topologies with robustness in the presence of noise, (ii) to approximate solutions of geometric partial differential equations (PDEs) on manifolds, and (iii) to solve equations for shape deformations such as curvature driven flows. We release codes and weights for using GNPs in the package geo_neural_op. This allows for incorporating our pre-trained GNPs as components for reuse within existing and new data processing pipelines. The GNPs also can be used as part of numerical solvers involving geometry or as part of methods for performing inference and other geometric tasks.
- [19] arXiv:2504.11619 (replaced) [pdf, html, other]
-
Title: Computing the Tropical Abel--Jacobi Transform and Tropical Distances for Metric GraphsComments: 51 pages, 9 figuresSubjects: Algebraic Geometry (math.AG); Metric Geometry (math.MG); Numerical Analysis (math.NA)
Metric graphs are important models for capturing the structure of complex data across various domains. While much effort has been devoted to extracting geometric and topological features from graph data, computational aspects of metric graphs as abstract tropical curves remains unexplored. In this paper, we present the first computational and machine learning-driven study of metric graphs from the perspective of tropical algebraic geometry. Specifically, we study the tropical Abel--Jacobi transform, a vectorization of points on a metric graph via the tropical Abel--Jacobi map into its associated flat torus, the tropical Jacobian. We develop algorithms to compute this transform and investigate how the resulting embeddings depend on different combinatorial models of the same metric graph.
Once embedded, we compute pairwise distances between points in the tropical Jacobian under two natural metrics: the tropical polarization distance and the Foster--Zhang distance. Computing these distances are generally NP-hard as they turn out to be linked to classical lattice problems in computational complexity, however, we identify a class of metric graphs where fast and explicit computations are feasible. For the general case, we propose practical algorithms for both exact and approximate distance matrix computations using lattice basis reduction and mixed-integer programming solvers. Our work lays the groundwork for future applications of tropical geometry and the tropical Abel--Jacobi transform in machine learning and data analysis.