Numerical Analysis
See recent articles
Showing new listings for Friday, 11 October 2024
- [1] arXiv:2410.07254 [pdf, html, other]
-
Title: Uniform accuracy of implicit-explicit Runge-Kutta methods for linear hyperbolic relaxation systemsComments: arXiv admin note: text overlap with arXiv:2306.08742 by other authorsSubjects: Numerical Analysis (math.NA); Computational Physics (physics.comp-ph)
In this paper, we study the uniform accuracy of implicit-explicit (IMEX) Runge-Kutta (RK) schemes for general linear hyperbolic relaxation systems satisfying the structural stability condition proposed in \cite{yong_singular_1999}. We establish the uniform stability and accuracy of a class of IMEX-RK schemes with spatial discretization using a Fourier spectral method. Our results demonstrate that the accuracy of the fully discretized schemes is independent of the relaxation time across all regimes. Numerical experiments on applications in traffic flows and kinetic theory verify our theoretical analysis.
- [2] arXiv:2410.07375 [pdf, html, other]
-
Title: Boundary-value problems of functional differential equations with state-dependent delaysComments: 44 pages, 2 figuresSubjects: Numerical Analysis (math.NA)
We prove convergence of piecewise polynomial collocation methods applied to periodic boundary value problems for functional differential equations with state-dependent delays. The state dependence of the delays leads to nonlinearities that are not locally Lipschitz continuous preventing the direct application of general abstract discretization theoretic frameworks. We employ a weaker form of differentiability, which we call mild differentiability, to prove that a locally unique solution of the functional differential equation is approximated by the solution of the discretized problem with the expected order.
An additional difficulty is that linearizations required for solving the discretized nonlinear problem with Newton iterations are not well defined. We show that Newton iterations still converge if one uses the linearization in regularized solutions. The Newton iterations' asymptotic convergence ratio is limited by the numerical discretization error. Thus, Newton iterations should show better convergence for approximations on finer meshes. - [3] arXiv:2410.07465 [pdf, html, other]
-
Title: Preconditioning Low Rank Generalized Minimal Residual Method (GMRES) for Implicit Discretizations of Matrix Differential EquationsSubjects: Numerical Analysis (math.NA)
This work proposes a new class of preconditioners for the low rank Generalized Minimal Residual Method (GMRES) for multiterm matrix equations arising from implicit timestepping of linear matrix differential equations. We are interested in computing low rank solutions to matrix equations, e.g. arising from spatial discretization of stiff partial differential equations (PDEs). The low rank GMRES method is a particular class of Krylov subspace method where the iteration is performed on the low rank factors of the solution. Such methods can exploit the low rank property of the solution to save on computational and storage cost.
Of critical importance for the efficiency and applicability of the low rank GMRES method is the availability of an effective low rank preconditioner that operates directly on the low rank factors of the solution and that can limit the iteration count and the maximal Krylov rank.
The preconditioner we propose here is based on the basis update and Galerkin (BUG) method, resulting from the dynamic low rank approximation. It is a nonlinear preconditioner for the low rank GMRES scheme that naturally operates on the low rank factors. Extensive numerical tests show that this new preconditioner is highly efficient in limiting iteration count and maximal Krylov rank. We show that the preconditioner performs well for general diffusion equations including highly challenging problems, e.g. high contrast, anisotropic equations. Further, it compares favorably with the state of the art exponential sum preconditioner. We also propose a hybrid BUG - exponential sum preconditioner based on alternating between the two preconditioners. - [4] arXiv:2410.07468 [pdf, html, other]
-
Title: Fast Real Evaluation Through Sound Mixed-Precision TuningSubjects: Numerical Analysis (math.NA); Mathematical Software (cs.MS)
Evaluating a real-valued expression to high precision is a key building block in computational mathematics, physics, and numerics. A typical implementation uses a uniform precision for each operation, and doubles that precision until the real result can be bounded to some sufficiently narrow interval. However, this is wasteful: usually only a few operations really need to be performed at high precision, and the bulk of the expression could use much lower precision. Uniform precision can also waste iterations discovering the necessary precision and then still overestimate by up to a factor of two. We propose to instead use mixed-precision interval arithmetic to evaluate real-valued expressions. A key challenge is deriving the mixed-precision assignment both soundly and quickly. To do so, we introduce a sound variation of error Taylor series and condition numbers, specialized to interval arithmetic, that can be evaluated with minimal overhead thanks to an "exponent trick". Our implementation, Reval, achieves a speed-up of 1.25x compared to the state-of-the-art Sollya tool, with the speed-up increasing to 2.99x on the most difficult input points. An examination of the precisions used with and without precision tuning shows that the speed-up results come from quickly assigning lower precisions for the majority of operations.
- [5] arXiv:2410.07703 [pdf, html, other]
-
Title: Time-domain direct sampling method for inverse electromagnetic scattering with a single incident sourceSubjects: Numerical Analysis (math.NA)
In this paper, we consider an inverse electromagnetic medium scattering problem of reconstructing unknown objects from time-dependent boundary measurements. A novel time-domain direct sampling method is developed for determining the locations of unknown scatterers by using only a single incident source. Notably, our method imposes no restrictions on the the waveform of the incident wave. Based on the Fourier-Laplace transform, we first establish the connection between the frequency-domain and the time-domain direct sampling method. Furthermore, we elucidate the mathematical mechanism of the imaging functional through the properties of modified Bessel functions. Theoretical justifications and stability analyses are provided to demonstrate the effectiveness of the proposed method. Finally, several numerical experiments are presented to illustrate the feasibility of our approach.
- [6] arXiv:2410.07723 [pdf, other]
-
Title: High-order discretized ACMS method for the simulation of finite-size two-dimensional photonic crystalsSubjects: Numerical Analysis (math.NA)
The computational complexity and efficiency of the approximate mode component synthesis (ACMS) method is investigated for the two-dimensional heterogeneous Helmholtz equations, aiming at the simulation of large but finite-size photonic crystals. The ACMS method is a Galerkin method that relies on a non-overlapping domain decomposition and special basis functions defined based on the domain decomposition. While, in previous works, the ACMS method was realized using first-order finite elements, we use an underlying hp-finite element method. We study the accuracy of the ACMS method for different wavenumbers, domain decompositions, and discretization parameters. Moreover, the computational complexity of the method is investigated theoretically and compared with computing times for an implementation based on the open source software package NGSolve. The numerical results indicate that, for relevant wavenumber regimes, the size of the resulting linear systems for the ACMS method remains moderate, such that sparse direct solvers are a reasonable choice. Moreover, the ACMS method exhibits only a weak dependence on the selected domain decomposition, allowing for greater flexibility in its choice. Additionally, the numerical results show that the error of the ACMS method achieves the predicted convergence rate for increasing wavenumbers. Finally, to display the versatility of the implementation, the results of simulations of large but finite-size photonic crystals with defects are presented.
- [7] arXiv:2410.08042 [pdf, other]
-
Title: {\varphi}-FD : A well-conditioned finite difference method inspired by {\varphi}-FEM for general geometries on elliptic PDEsSubjects: Numerical Analysis (math.NA)
This paper presents a new finite difference method, called {\varphi}-FD, inspired by the {\phi}-FEM approach for solving elliptic partial differential equations (PDEs) on general geometries. The proposed method uses Cartesian grids, ensuring simplicity in implementation. Moreover, contrary to the previous finite difference scheme on non-rectangular domain, the associated matrix is well-conditioned. The use of a level-set function for the geometry description makes this approach relatively flexible. We prove the quasi-optimal convergence rates in several norms and the fact that the matrix is well-conditioned. Additionally, the paper explores the use of multigrid techniques to further accelerate the computation. Finally, numerical experiments in both 2D and 3D validate the performance of the {\varphi}-FD method compared to standard finite element methods and the Shortley-Weller approach.
New submissions (showing 7 of 7 entries)
- [8] arXiv:2410.07605 (cross-list from cs.CV) [pdf, html, other]
-
Title: A Variational Bayesian Inference Theory of Elasticity and Its Mixed Probabilistic Finite Element Method for Inverse Deformation Solutions in Any DimensionSubjects: Computer Vision and Pattern Recognition (cs.CV); Machine Learning (cs.LG); Numerical Analysis (math.NA)
In this work, we have developed a variational Bayesian inference theory of elasticity, which is accomplished by using a mixed Variational Bayesian inference Finite Element Method (VBI-FEM) that can be used to solve the inverse deformation problems of continua. In the proposed variational Bayesian inference theory of continuum mechanics, the elastic strain energy is used as a prior in a Bayesian inference network, which can intelligently recover the detailed continuum deformation mappings with only given the information on the deformed and undeformed continuum body shapes without knowing the interior deformation and the precise actual boundary conditions, both traction as well as displacement boundary conditions, and the actual material constitutive relation. Moreover, we have implemented the related finite element formulation in a computational probabilistic mechanics framework. To numerically solve mixed variational problem, we developed an operator splitting or staggered algorithm that consists of the finite element (FE) step and the Bayesian learning (BL) step as an analogue of the well-known the Expectation-Maximization (EM) algorithm. By solving the mixed probabilistic Galerkin variational problem, we demonstrated that the proposed method is able to inversely predict continuum deformation mappings with strong discontinuity or fracture without knowing the external load conditions. The proposed method provides a robust machine intelligent solution for the long-sought-after inverse problem solution, which has been a major challenge in structure failure forensic pattern analysis in past several decades. The proposed method may become a promising artificial intelligence-based inverse method for solving general partial differential equations.
- [9] arXiv:2410.07646 (cross-list from math.OC) [pdf, html, other]
-
Title: R-Adaptive Mesh Optimization to Enhance Finite Element Basis CompressionComments: 22 pagesSubjects: Optimization and Control (math.OC); Numerical Analysis (math.NA)
Modern computing systems are capable of exascale calculations, which are revolutionizing the development and application of high-fidelity numerical models in computational science and engineering. While these systems continue to grow in processing power, the available system memory has not increased commensurately, and electrical power consumption continues to grow. A predominant approach to limit the memory usage in large-scale applications is to exploit the abundant processing power and continually recompute many low-level simulation quantities, rather than storing them. However, this approach can adversely impact the throughput of the simulation and diminish the benefits of modern computing architectures. We present two novel contributions to reduce the memory burden while maintaining performance in simulations based on finite element discretizations. The first contribution develops dictionary-based data compression schemes that detect and exploit the structure of the discretization, due to redundancies across the finite element mesh. These schemes are shown to reduce memory requirements by more than 99 percent on meshes with large numbers of nearly identical mesh cells. For applications where this structure does not exist, our second contribution leverages a recently developed augmented Lagrangian sequential quadratic programming algorithm to enable r-adaptive mesh optimization, with the goal of enhancing redundancies in the mesh. Numerical results demonstrate the effectiveness of the proposed methods to detect, exploit and enhance mesh structure on examples inspired by large-scale applications.
- [10] arXiv:2410.07776 (cross-list from math.AP) [pdf, html, other]
-
Title: Median filter method for mean curvature flow using a random Jacobi algorithmComments: 56 pages, 7 figuresSubjects: Analysis of PDEs (math.AP); Numerical Analysis (math.NA)
We present an efficient scheme for level set mean curvature flow using a domain discretization and median filters. For this scheme, we show convergence in $L^\infty$-norm under mild assumptions on the number of points in the discretization. In addition, we strengthen the weak convergence result for the MBO thresholding scheme applied to data clustering of Lelmi and one of the authors. This is done through a strong convergence of the discretized heat flow in the optimal regime. Different boundary conditions are also discussed.
Cross submissions (showing 3 of 3 entries)
- [11] arXiv:2308.02467 (replaced) [pdf, html, other]
-
Title: Element learning: a systematic approach of accelerating finite element-type methods via machine learning, with applications to radiative transferSubjects: Numerical Analysis (math.NA)
In this paper, we propose a systematic approach for accelerating finite element-type methods by machine learning for the numerical solution of partial differential equations (PDEs). The main idea is to use a neural network to learn the solution map of the PDEs and to do so in an element-wise fashion. This map takes input of the element geometry and the PDEs' parameters on that element, and gives output of two operators -- (1) the in2out operator for inter-element communication, and (2) the in2sol operator (Green's function) for element-wise solution recovery. A significant advantage of this approach is that, once trained, this network can be used for the numerical solution of the PDE for any domain geometry and any parameter distribution without retraining. Also, the training is significantly simpler since it is done on the element level instead on the entire domain. We call this approach element learning. This method is closely related to hybridizbale discontinuous Galerkin (HDG) methods in the sense that the local solvers of HDG are replaced by machine learning approaches. Numerical tests are presented for an example PDE, the radiative transfer equation, in a variety of scenarios with idealized or realistic cloud fields, with smooth or sharp gradient in the cloud boundary transition. Under a fixed accuracy level of $10^{-3}$ in the relative $L^2$ error, and polynomial degree $p=6$ in each element, we observe an approximately 5 to 10 times speed-up by element learning compared to a classical finite element-type method.
- [12] arXiv:2310.14290 (replaced) [pdf, html, other]
-
Title: Data-driven Morozov regularization of inverse problemsSubjects: Numerical Analysis (math.NA)
The solution of inverse problems is crucial in various fields such as medicine, biology, and engineering, where one seeks to find a solution from noisy observations. These problems often exhibit non-uniqueness and ill-posedness, resulting in instability under noise with standard methods. To address this, regularization techniques have been developed to balance data fitting and prior information. Recently, data-driven variational regularization methods have emerged, mainly analyzed within the framework of Tikhonov regularization, termed Network Tikhonov (NETT). This paper introduces Morozov regularization combined with a learned regularizer, termed DD-Morozov regularization. Our approach employs neural networks to define non-convex regularizers tailored to training data, enabling a convergence analysis in the non-convex context with noise-dependent regularizers. We also propose a refined training strategy that improves adaptation to ill-posed problems compared to NETT's original strategy, which primarily focuses on addressing non-uniqueness. We present numerical results in the context of attenuation correction in photoacoustic tomography, comparing DD-Morozov regularization with NETT using the same trained regularizer, both with and without an additional total variation regularizer.
- [13] arXiv:2401.04783 (replaced) [pdf, html, other]
-
Title: Hyperbolic Machine Learning Moment Closures for the BGK EquationsComments: 30 pages, 7 figuresSubjects: Numerical Analysis (math.NA); Machine Learning (cs.LG); Computational Physics (physics.comp-ph)
We introduce a hyperbolic closure for the Grad moment expansion of the Bhatnagar-Gross-Krook's (BGK) kinetic model using a neural network (NN) trained on BGK's moment data. This closure is motivated by the exact closure for the free streaming limit that we derived in our paper on closures in transport \cite{Huang2022-RTE1}. The exact closure relates the gradient of the highest moment to the gradient of four lower moments. As with our past work, the model presented here learns the gradient of the highest moment in terms of the coefficients of gradients for all lower ones. By necessity, this means that the resulting hyperbolic system is not conservative in the highest moment. For stability, the output layers of the NN are designed to enforce hyperbolicity and Galilean invariance. This ensures the model can be run outside of the training window of the NN. Unlike our previous work on radiation transport that dealt with linear models, the BGK model's nonlinearity demanded advanced training tools. These comprised an optimal learning rate discovery, one cycle training, batch normalization in each neural layer, and the use of the \texttt{AdamW} optimizer. To address the non-conservative structure of the hyperbolic model, we adopt the FORCE numerical method to achieve robust solutions. This results in a comprehensive computing model combining learned closures with methods for solving hyperbolic models. The proposed model can capture accurate moment solutions across a broad spectrum of Knudsen numbers. Our paper details the multi-scale model construction and is run on a range of test problems.
- [14] arXiv:2404.18505 (replaced) [pdf, other]
-
Title: R3MG: R-tree based agglomeration of polytopal grids with applications to multilevel methodsComments: 25 pages, 21 figures, 5 tablesSubjects: Numerical Analysis (math.NA)
We present a novel approach to perform agglomeration of polygonal and polyhedral grids based on spatial indices. Agglomeration strategies are a key ingredient in polytopal methods for PDEs as they are used to generate (hierarchies of) computational grids from an initial grid. Spatial indices are specialized data structures that significantly accelerate queries involving spatial relationships in arbitrary space dimensions. We show how the construction of the R-tree spatial database of an arbitrary fine grid offers a natural and efficient agglomeration strategy with the following characteristics: i) the process is fully automated, robust, and dimension-independent, ii) it automatically produces a balanced and nested hierarchy of agglomerates, and iii) the shape of the agglomerates is tightly close to the respective axis aligned bounding boxes. Moreover, the R-tree approach provides a full hierarchy of nested agglomerates which permits fast query and allows for efficient geometric multigrid methods to be applied also to those cases where a hierarchy of grids is not present at construction time. We present several examples based on polygonal discontinuous Galerkin methods, confirming the effectiveness of our approach in the context of challenging three-dimensional geometries and the design of geometric multigrid preconditioners.
- [15] arXiv:2407.05230 (replaced) [pdf, html, other]
-
Title: Matrix perturbation bounds via contour bootstrappingSubjects: Numerical Analysis (math.NA); Statistics Theory (math.ST)
Matrix perturbation bounds play an essential role in the design and analysis of spectral algorithms. In this paper, we introduce a new method to deduce matrix perturbation bounds, which we call "contour bootstrapping". As applications, we work out several new bounds for eigensubspace computation and low rank approximation. Next, we use these bounds to study utility problems in the area of differential privacy.
- [16] arXiv:2410.06850 (replaced) [pdf, html, other]
-
Title: A robust solver for large-scale heat transfer topology optimizationSubjects: Numerical Analysis (math.NA)
This paper presents a large-scale parallel solver, specifically designed to tackle the challenges of solving high-dimensional and high-contrast linear systems in heat transfer topology optimization. The solver incorporates an interpolation technique to accelerate convergence in high-resolution domains, along with a multiscale multigrid preconditioner to handle complex coefficient fields with significant contrast. All modules of the optimization solver are implemented on a high performance computing cluster by the PETSc numerical library. Through a series of numerical investigations, we demonstrate the effectiveness of our approach in enhancing convergence and robustness during the optimization process, particularly in high-contrast scenarios with resolutions up to $1024^3$. Our performance results indicate that the proposed preconditioner achieves over $2\times$ speedup against the default algebraic multigrid in PETSc for high-contrast cases.
- [17] arXiv:2302.02665 (replaced) [pdf, html, other]
-
Title: Adjoint Method in PDE-based Image CompressionSubjects: Analysis of PDEs (math.AP); Numerical Analysis (math.NA)
We consider a shape optimization based method for finding the best interpolation data in the compression of images with noise. The aim is to reconstruct missing regions by means of minimizing a data fitting term in an $L^p$-norm between original images and their reconstructed counterparts using linear diffusion PDE-based inpainting. Reformulating the problem as a constrained optimization over sets (shapes), we derive the topological asymptotic expansion of the considered shape functionals with respect to the insertion of small ball (a single pixel) using the adjoint method. Based on the achieved distributed topological shape derivatives, we propose a numerical approach to determine the optimal set and present numerical experiments showing, the efficiency of our method. Numerical computations are presented that confirm the usefulness of our theoretical findings for PDE-based image compression.
- [18] arXiv:2404.16015 (replaced) [pdf, html, other]
-
Title: Neural Operators Learn the Local Physics of MagnetohydrodynamicsComments: 48 pages, 24 figuresSubjects: Computational Physics (physics.comp-ph); Artificial Intelligence (cs.AI); Numerical Analysis (math.NA)
Magnetohydrodynamics (MHD) plays a pivotal role in describing the dynamics of plasma and conductive fluids, essential for understanding phenomena such as the structure and evolution of stars and galaxies, and in nuclear fusion for plasma motion through ideal MHD equations. Solving these hyperbolic PDEs requires sophisticated numerical methods, presenting computational challenges due to complex structures and high costs. Recent advances introduce neural operators like the Fourier Neural Operator (FNO) as surrogate models for traditional numerical analyses. This study explores a modified Flux Fourier neural operator model to approximate the numerical flux of ideal MHD, offering a novel approach that outperforms existing neural operator models by enabling continuous inference, generalization outside sampled distributions, and faster computation compared to classical numerical schemes.
- [19] arXiv:2409.09906 (replaced) [pdf, html, other]
-
Title: Variance-reduced first-order methods for deterministically constrained stochastic nonconvex optimization with strong convergence guaranteesComments: Significantly improves the previous complexity resultsSubjects: Optimization and Control (math.OC); Machine Learning (cs.LG); Numerical Analysis (math.NA); Machine Learning (stat.ML)
In this paper, we study a class of deterministically constrained stochastic optimization problems. Existing methods typically aim to find an $\epsilon$-stochastic stationary point, where the expected violations of both constraints and first-order stationarity are within a prescribed accuracy $\epsilon$. However, in many practical applications, it is crucial that the constraints be nearly satisfied with certainty, making such an $\epsilon$-stochastic stationary point potentially undesirable due to the risk of significant constraint violations. To address this issue, we propose single-loop variance-reduced stochastic first-order methods, where the stochastic gradient of the stochastic component is computed using either a truncated recursive momentum scheme or a truncated Polyak momentum scheme for variance reduction, while the gradient of the deterministic component is computed exactly. Under the error bound condition with a parameter $\theta \geq 1$ and other suitable assumptions, we establish that these methods respectively achieve a sample and first-order operation complexity of $\widetilde O(\epsilon^{-\max\{\theta+2, 2\theta\}})$ and $\widetilde O(\epsilon^{-\max\{4, 2\theta\}})$ for finding a stronger $\epsilon$-stochastic stationary point, where the constraint violation is within $\epsilon$ with certainty, and the expected violation of first-order stationarity is within $\epsilon$. For $\theta=1$, these complexities reduce to $\widetilde O(\epsilon^{-3})$ and $\widetilde O(\epsilon^{-4})$ respectively, which match, up to a logarithmic factor, the best-known complexities achieved by existing methods for finding an $\epsilon$-stochastic stationary point of unconstrained smooth stochastic optimization problems.