Numerical Analysis
See recent articles
Showing new listings for Friday, 25 July 2025
- [1] arXiv:2507.17916 [pdf, html, other]
-
Title: Recovery Thresholding Hyperinterpolations in Signal ProcessingComments: 16 pages, 3 figures, 5 tablesSubjects: Numerical Analysis (math.NA)
This paper introduces recovery thresholding hyperinterpolations, a novel class of methods for sparse signal reconstruction in the presence of noise. We develop a framework that integrates thresholding operators--including hard thresholding, springback, and Newton thresholding--directly into the hyperinterpolation structure to maintain sparsity during signal recovery. Our approach leverages Newton's method to minimize one-dimensional nonconvex functions, which we then extend to solve multivariable nonconvex regularization problems. The proposed methods demonstrate robust performance in reconstructing signals corrupted by both Gaussian and impulse noise. Through numerical experiments, we validate the effectiveness of these recovery thresholding hyperinterpolations for signal reconstruction and function denoising applications, showing their advantages over traditional approaches in preserving signal sparsity while achieving accurate recovery.
- [2] arXiv:2507.17928 [pdf, html, other]
-
Title: A novel finite element method for simulating surface plasmon polaritons on complex graphene sheetsSubjects: Numerical Analysis (math.NA)
Surface plasmon polaritons (SPPs) are generated on the graphene surface, and provide a window into the nano-optical and electrodynamic response of their host material and its dielectric environment. An accurate simulation of SPPs presents several unique challenges, since SPPs often occur at complex interfaces between materials of different dielectric constants and appropriate boundary conditions at the graphene interfaces are crucial. Here we develop a simplified graphene model and propose a new finite element method accordingly. Stability for the continuous model is established, and extensive numerical results are presented to demonstrate that the new model can capture the SPPs very well for various complex graphene sheets.
- [3] arXiv:2507.18235 [pdf, other]
-
Title: A stabilized Two-Step Formulation of Maxwell's Equations in the time-domainComments: 6 pages, 9 figuresSubjects: Numerical Analysis (math.NA); Computational Engineering, Finance, and Science (cs.CE)
Simulating electromagnetic fields across broad frequency ranges is challenging due to numerical instabilities at low frequencies. This work extends a stabilized two-step formulation of Maxwell's equations to the time-domain. Using a Galerkin discretization in space, we apply two different time-discretization schemes that are tailored to the first- and second-order in time partial differential equations of the two-step solution procedure used here. To address the low-frequency instability, we incorporate a generalized tree-cotree gauge that removes the singularity of the curl-curl operator, ensuring robustness even in the static limit. Numerical results on academic and application-oriented 3D problems confirm stability, accuracy, and the method's applicability to nonlinear, temperature-dependent materials.
- [4] arXiv:2507.18282 [pdf, other]
-
Title: EigenWave: An Optimal O(N) Method for Computing Eigenvalues and Eigenvectors by Time-Filtering the Wave EquationSubjects: Numerical Analysis (math.NA)
An algorithm named EigenWave is described to compute eigenvalues and eigenvectors of elliptic boundary value problems. The algorithm, based on the recently developed WaveHoltz scheme, solves a related time-dependent wave equation as part of an iteration. At each iteration, the solution to the wave equation is filtered in time. As the iteration progresses, the filtered solution generally contains relatively larger and larger proportions of eigenmodes whose eigenvalues are near a chosen target frequency (target eigenvalue). The ability to choose an arbitrary target frequency enables the computation of eigenvalues anywhere in the spectrum, without the need to invert an indefinite matrix, as is common with other approaches. Furthermore, the iteration can be embedded within a matrix-free Arnoldi algorithm, which enables the efficient computation of multiple eigenpairs near the target frequency. For efficiency, the time-dependent wave equation can be solved with implicit time-stepping and only about $10$ time-steps per-period are needed, independent of the mesh spacing. When the (definite) implicit time-stepping equations are solved with a multigrid algorithm, the cost of the resulting EigenWave scheme scales linearly with the number of grid points $N$ as the mesh is refined, giving an optimal $O(N)$ algorithm. The approach is demonstrated by finding eigenpairs of the Laplacian in complex geometry using overset grids. Results in two and three space dimensions are presented using second-order and fourth-order accurate approximations.
- [5] arXiv:2507.18318 [pdf, other]
-
Title: Parametric design and adaptive sizing of lattice structures for 3d additive manufacturingJorge Manuel Mercado-Colmenero, Daniel Diaz - Perete, Miguel Angel Rubio- Paramio, Cristina Martin-DonateJournal-ref: Manuel Mercado-Colmenero, J., D\'iaz-Perete, D., \'Angel Rubio-Paramio, M., & Martin-Do\~nate, C. (2025). PARAMETRIC DESIGN AND ADAPTIVE SIZING OF LATTICE STRUCTURES FOR 3D ADDITIVE MANUFACTURING. DYNA-Ingenier\'ia e Industria, 100(1)Subjects: Numerical Analysis (math.NA)
The present research is developed into the realm of industrial design engineering and additive manufacturing by introducing a parametric design model and adaptive mechanical analysis for a new lattice structure, with a focus on 3D additive manufacturing of complex parts. Focusing on the land-scape of complex parts additive manufacturing, this research proposes geometric parameterization, mechanical adaptive sizing, and numerical validation of a novel lattice structure to optimize the final printed part volume and mass, as well as its structural rigidity. The topology of the lattice structures exhibited pyramidal geometry. Complete parameterization of the lattice structure ensures that the known geometric parameters adjust to defined restrictions, enabling dynamic adaptability based on its load states and boundary conditions, thereby enhancing its mechanical performance. The core methodology integrates analytical automation with mechanical analysis by employing a model based in two-dimensional beam elements. The dimensioning of the lattice structure is analyzed using rigidity models of its sub-elements, providing an evaluation of its global structural behavior after applying the superposition principle. Numerical validation was performed to validate the proposed analytical model. This step ensures that the analytical model defined for dimensioning the lattice structure adjusts to its real mechanical behavior and allows its validation. The present manuscript aims to advance additive manufacturing methodologies by offering a systematic and adaptive approach to lattice structure design. Parametric and adaptive techniques foster new industrial design engineering methods, enabling the dynamic tailoring of lattice structures to meet their mechanical demands and enhance their overall efficiency and performance.
- [6] arXiv:2507.18443 [pdf, html, other]
-
Title: On MAP estimates and source conditions for drift identification in SDEsSubjects: Numerical Analysis (math.NA)
We consider the inverse problem of identifying the drift in an SDE from $n$ observations of its solution at $M+1$ distinct time points. We derive a corresponding MAP estimate, we prove differentiability properties as well as a so-called tangential cone condition for the forward operator, and we review the existing theory for related problems, which under a slightly stronger tangential cone condition would additionally yield convergence rates for the MAP estimate as $n\to\infty$. Numerical simulations in 1D indicate that such convergence rates indeed hold true.
- [7] arXiv:2507.18466 [pdf, html, other]
-
Title: Solution of Least Squares Problems with Randomized Preconditioned Normal EquationsSubjects: Numerical Analysis (math.NA)
We consider the solution of full column-rank least squares problems by means of normal equations that are preconditioned, symmetrically or non-symmetrically, with a randomized preconditioner. With an effective preconditioner, the solutions from the preconditioned normal equations are almost as accurate as those from the QR-based Matlab backslash (mldivide) command -- even for highly illconditioned matrices. This means the accuracy of the preconditioned normal equations depends on the residual of the original least squares problem. We present non-intuitive but realistic perturbation bounds for the relative error in the computed solutions and show that, with an effective preconditioner, these bounds are essentially equal to the perturbation bound for the original least squares problem. Probabilitistic condition number bounds corroborate the effectiveness of the randomized preconditioner computed with small amounts of sampling.
- [8] arXiv:2507.18491 [pdf, html, other]
-
Title: Fast Multipole Method for Maxwell's Equations in Layered MediaComments: Submitted to SIAM Journal on Scientific ComputingSubjects: Numerical Analysis (math.NA)
We present a fast multipole method (FMM) for solving Maxwell's equations in three-dimensional (3-D) layered media, based on the magnetic vector potential $\boldsymbol A$ under the Lorenz gauge, to derive the layered dyadic Green's function. The dyadic Green's function is represented using three scalar Helmholtz layered Green's functions, with all interface-induced reaction field components expressed through a unified integral representation. By introducing equivalent polarization images for sources and effective locations for targets to reflect the actual transmission distance of different reaction field components, multiple expansions (MEs) and local expansions (LEs) are derived for the far-field governed by actual transmission distance. To further enhance computational efficiency and numerical stability, we employ a Chebyshev polynomial expansion of the associated Legendre functions to speed up the calculation of multipole-to-local (M2L) expansion translations. Finally, leveraging the FMM framework of the Helmholtz equation in 3-D layered media, we develop a FMM for the dyadic Green's function of Maxwell's equations in layered media. Numerical experiments demonstrate the $\mathcal O(N\log N)$-complexity of the resulting FMM method, and rapid convergence for interactions of low-frequency electromagnetic wave sources in 3-D layered media.
New submissions (showing 8 of 8 entries)
- [9] arXiv:2507.17815 (cross-list from hep-th) [pdf, html, other]
-
Title: Analytic Regression of Feynman Integrals from High-Precision Numerical SamplingSubjects: High Energy Physics - Theory (hep-th); High Energy Physics - Phenomenology (hep-ph); Numerical Analysis (math.NA)
In mathematics or theoretical physics one is often interested in obtaining an exact analytic description of some data which can be produced, in principle, to arbitrary accuracy. For example, one might like to know the exact analytical form of a definite integral. Such problems are not well-suited to numerical symbolic regression, since typical numerical methods lead only to approximations. However, if one has some sense of the function space in which the analytic result should lie, it is possible to deduce the exact answer by judiciously sampling the data at a sufficient number of points with sufficient precision. We demonstrate how this can be done for the computation of Feynman integrals. We show that by combining high-precision numerical integration with analytic knowledge of the function space one can often deduce the exact answer using lattice reduction. A number of examples are given as well as an exploration of the trade-offs between number of datapoints, number of functional predicates, precision of the data, and compute. This method provides a bottom-up approach that neatly complements the top-down Landau-bootstrap approach of trying to constrain the exact answer using the analytic structure alone. Although we focus on the application to Feynman integrals, the techniques presented here are more general and could apply to a wide range of problems where an exact answer is needed and the function space is sufficiently well understood.
- [10] arXiv:2507.17887 (cross-list from cs.LG) [pdf, html, other]
-
Title: Fourier Neural Operators for Non-Markovian Processes:Approximation Theorems and ExperimentsSubjects: Machine Learning (cs.LG); Numerical Analysis (math.NA)
This paper introduces an operator-based neural network, the mirror-padded Fourier neural operator (MFNO), designed to learn the dynamics of stochastic systems. MFNO extends the standard Fourier neural operator (FNO) by incorporating mirror padding, enabling it to handle non-periodic inputs. We rigorously prove that MFNOs can approximate solutions of path-dependent stochastic differential equations and Lipschitz transformations of fractional Brownian motions to an arbitrary degree of accuracy. Our theoretical analysis builds on Wong--Zakai type theorems and various approximation techniques. Empirically, the MFNO exhibits strong resolution generalization--a property rarely seen in standard architectures such as LSTMs, TCNs, and DeepONet. Furthermore, our model achieves performance that is comparable or superior to these baselines while offering significantly faster sample path generation than classical numerical schemes.
- [11] arXiv:2507.18083 (cross-list from astro-ph.GA) [pdf, other]
-
Title: Numerical Study of Bar Suppression in Galaxy Models Due to Disc HeatingComments: Accepted for publication in the Galaxies journal. 34 pages (including references),12 figures and 2 tablesJournal-ref: Galaxies 2025, 13(2), 45Subjects: Astrophysics of Galaxies (astro-ph.GA); Numerical Analysis (math.NA)
The process of bar formation, evolution and destruction is still a controversial topic regarding galaxy dynamics. Numerical simulations show that these phenomena strongly depend on physical and numerical parameters. In this work, we study the combined influence of the softening parameter, $\epsilon$ and disc mass fraction, $m_{\mathrm{d}}$ on the formation and evolution of bars in isolated disc-halo models via $N$-body simulations with different particle resolutions. Previous studies indicate that the bar strength depends on $m_{\mathrm{d}}$ as $\propto m_{\mathrm{d}}^{-1}$, which is seen as a delay in bar formation. However, the distorsion parameter, $\eta$, which measures the bar's momentum through time, shows that an increase in $m_{\mathrm{d}}$ does not always induce a delay in bar formation. This suggests that $\epsilon$ interact to either enhance or weaken the bar. Moreover, numerical heating dominates in models with small softening values, creating highly accelerated particles at the centre of discs, regardless of $m_{\mathrm{d}}$ or resolution. These enhanced particle accelerations produce chaotic orbits for $\epsilon \leq 5\,$pc, resulting in bar suppression due to collisional dynamics in the centre. In our high resolution models ($N \approx 10^{7}$), small softening values are incapable of reproducing the bar instability. The role of disc mass is as follows: increasing $m_{\mathrm{d}}$ for moderate $\epsilon$ ($\geq 10\,$pc) reduces the amount of drift in the acceleration profile, without affecting the bar's behaviour. Models with lower $m_{\mathrm{d}}$ values coupled with small softening values, have an excess of highly accelerated particles, introducing unwanted effects into otherwise reliable simulations. Finally, we show that the evolution of the disc's vertical acceleration profile is a reliable indicator of numerical heating introduced by $\epsilon$ and the bar.
- [12] arXiv:2507.18346 (cross-list from cs.LG) [pdf, html, other]
-
Title: Low-rank adaptive physics-informed HyperDeepONets for solving differential equationsComments: 14 pages, 6 figures, 5 tablesSubjects: Machine Learning (cs.LG); Numerical Analysis (math.NA)
HyperDeepONets were introduced in Lee, Cho and Hwang [ICLR, 2023] as an alternative architecture for operator learning, in which a hypernetwork generates the weights for the trunk net of a DeepONet. While this improves expressivity, it incurs high memory and computational costs due to the large number of output parameters required. In this work we introduce, in the physics-informed machine learning setting, a variation, PI-LoRA-HyperDeepONets, which leverage low-rank adaptation (LoRA) to reduce complexity by decomposing the hypernetwork's output layer weight matrix into two smaller low-rank matrices. This reduces the number of trainable parameters while introducing an extra regularization of the trunk networks' weights. Through extensive experiments on both ordinary and partial differential equations we show that PI-LoRA-HyperDeepONets achieve up to 70\% reduction in parameters and consistently outperform regular HyperDeepONets in terms of predictive accuracy and generalization.
- [13] arXiv:2507.18434 (cross-list from math.CO) [pdf, html, other]
-
Title: Guessing sequences of eigenvectors for LMPs defining spectrahedral relaxations of Eulerian rigidly convex setsComments: 51 pages. Preprint extracted from a selection, rewrite and recombination of several sections and chapters from my PhD thesis. For more possible lines of research in these related directions, we direct the interested reader to arXiv:2503.04628 and to arXiv:2507.03800Subjects: Combinatorics (math.CO); Numerical Analysis (math.NA); Optimization and Control (math.OC)
Stable multivariate Eulerian polynomials were introduced by Brändén. Particularizing some variables, it is possible to extract real zero multivariate Eulerian polynomials from them. These real zero multivariate Eulerian polynomials can be fed into constructions of spectrahedral relaxations providing therefore approximations to the (Eulerian) rigidly convex sets defined by these polynomials. The accuracy of these approximations is measured through the behaviour in the diagonal, where the usual univariate Eulerian polynomials sit. In particular, in this sense, the accuracy of the global spectrahedral approximation produced by the spectrahedral relaxation can be measured in terms of bounds for the extreme roots of univariate Eulerian polynomials. The bounds thus obtained beat the previous bounds found in the literature. However, the bound explicitly studied and obtained before beat the previously known bounds by a quantity going to $0$ when $n$ goes to infinity. Here we use numerical experiments to construct a sequence of vectors providing a (linearized) bound whose difference with the previous known bounds is a growing exponential function (going therefore fast to infinity when $n$ grows). This allows us to establish a better (diagonal) measure of accuracy for the spectrahedral relaxation of the Eulerian rigidly convex sets. In particular, we will achieve this by linearizing through the sequence of vectors $\{(y,(-2^{m-i})_{i=3}^{m},(0,\frac{1}{2}),(1)_{i=1}^{m})\in\mathbb{R}^{n+1}\}_{n=1}^{\infty}$ for even $n=2m$.
- [14] arXiv:2507.18573 (cross-list from math.DG) [pdf, html, other]
-
Title: Jacobi Hamiltonian IntegratorsComments: 32 pagesSubjects: Differential Geometry (math.DG); Mathematical Physics (math-ph); Numerical Analysis (math.NA); Symplectic Geometry (math.SG)
We develop a method of constructing structure-preserving integrators for Hamiltonian systems in Jacobi manifolds. Hamiltonian mechanics, rooted in symplectic and Poisson geometry, has long provided a foundation for modelling conservative systems in classical physics. Jacobi manifolds, generalizing both contact and Poisson manifolds, extend this theory and are suitable for incorporating time-dependent, dissipative and thermodynamic phenomena.
Building on recent advances in geometric integrators - specifically Poisson Hamiltonian Integrators (PHI), which preserve key features of Poisson systems - we propose a construction of Jacobi Hamiltonian Integrators. Our approach explores the correspondence between Jacobi and homogeneous Poisson manifolds, with the aim of extending the PHI techniques while ensuring preservation of the homogeneity structure.
This work develops the theoretical tools required for this generalization and outlines a numerical integration technique compatible with Jacobi dynamics. By focusing on the homogeneous Poisson perspective rather than on direct contact realizations, we provide a clear pathway for structure-preserving integration of time-dependent and dissipative systems within the Jacobi framework.
Cross submissions (showing 6 of 6 entries)
- [15] arXiv:1808.02391 (replaced) [pdf, other]
-
Title: Energy-preserving continuous-stage partitioned Runge-Kutta methodsComments: The paper needs to be improved.Subjects: Numerical Analysis (math.NA)
In this paper, we present continuous-stage partitioned Runge-Kutta (csPRK) methods for energy-preserving integration of Hamiltonian systems. A sufficient condition for the energy preservation of the csPRK methods is derived. It is shown that the presented condition contains the existing condition for energy-preserving continuous-stage Runge-Kutta methods as a special case. A noticeable and interesting result is that when we use the simplifying assumptions of order conditions and the normalized shifted Legendre polynomials for constructing high-order energy-preserving csPRK methods, both the Butcher "weight" coefficients $B_\tau$ and $\widehat{B}_\tau$ must be equal to $1$. As illustrative examples, new energy-preserving integrators are acquired by virtue of the presented condition, and for the sake of verifying our theoretical results, some numerical experiments are reported.
- [16] arXiv:1808.08451 (replaced) [pdf, other]
-
Title: Energy-preserving continuous-stage Runge-Kutta-Nyström methodsComments: The paper needs to be improved.Subjects: Numerical Analysis (math.NA)
Many practical problems can be described by second-order system $\ddot{q}=-M\nabla U(q)$, in which people give special emphasis to some invariants with explicit physical meaning, such as energy, momentum, angular momentum, etc. However, conventional numerical integrators for such systems will fail to preserve any of these quantities which may lead to qualitatively incorrect numerical solutions. This paper is concerned with the development of energy-preserving continuous-stage Runge-Kutta-Nyström (csRKN) methods for solving second-order systems. Sufficient conditions for csRKN methods to be energy-preserving are presented and it is proved that all the energy-preserving csRKN methods satisfying these sufficient conditions can be essentially induced by energy-preserving continuous-stage partitioned Runge-Kutta methods. Some illustrative examples are given and relevant numerical results are reported.
- [17] arXiv:2407.13066 (replaced) [pdf, other]
-
Title: Fast And Scalable FFT-Based GPU-Accelerated Algorithms for Block-Triangular Toeplitz Matrices With Application to Linear Inverse Problems Governed by Autonomous Dynamical SystemsSubjects: Numerical Analysis (math.NA)
We present an efficient and scalable algorithm for performing matrix-vector multiplications ("matvecs") for block Toeplitz matrices. Such matrices, which are shift-invariant with respect to their blocks, arise in the context of solving inverse problems governed by autonomous systems, and time-invariant systems in particular. In this article, we consider inverse problems that infer unknown parameters from observational data of a linear time-invariant dynamical system given in the form of partial differential equations (PDEs). Matrix-free Newton-conjugate-gradient methods are often the gold standard for solving these inverse problems, but they require numerous actions of the Hessian on a vector. Matrix-free adjoint-based Hessian matvecs require solution of a pair of linearized forward/adjoint PDE solves per Hessian action, which may be prohibitive for large-scale inverse problems. Time invariance of the forward PDE problem leads to a block Toeplitz structure of the discretized parameter-to-observable (p2o) map defining the mapping from inputs (parameters) to outputs (observables) of the PDEs. This block Toeplitz structure enables us to exploit two key properties: (1) compact storage of the p2o map and its adjoint; and (2) efficient fast Fourier transform (FFT)-based Hessian matvecs. The proposed algorithm is mapped onto large multi-GPU clusters and achieves more than 80 percent of peak bandwidth on NVIDIA A100 GPUs. Excellent weak scaling is shown for up to 48 A100 GPUs. For the targeted problems, the implementation executes Hessian matvecs within fractions of a second, which is orders of magnitude faster than can be achieved by conventional matrix-free Hessian matvecs via forward/adjoint PDE solves.
- [18] arXiv:2409.03872 (replaced) [pdf, html, other]
-
Title: Continuous data assimilation for hydrodynamics: consistent discretization and application to moment recoveryComments: Updated version in Journal of Computational PhysicsSubjects: Numerical Analysis (math.NA)
Motivated by the challenge of moment recovery in hydrodynamic approximation in kinetic theory, we propose a data-driven approach for the hydrodynamic models. Inspired by continuous data assimilation, our method introduces a relaxation-based nudging system coupled with a novel discretization technique. This approach facilitates the simultaneous recovery of both the force term and a high-resolution solution from sparsely observed data. To address potential numerical artifacts, we use kernel regression to fit the observed data. We also analyze the convergence of the proposed nudging system under both full and partial data scenarios. When applied to moment systems, the source term involves the derivative of higher-order moments, our approach serves as a crucial step for data preparation in machine-learning based moment closure models. Multiple numerical experiments demonstrate the effectiveness of our algorithm, and we discuss its potential extension to high-dimensional systems.
- [19] arXiv:2410.02540 (replaced) [pdf, html, other]
-
Title: $hp$-error analysis of mixed-order hybrid high-order methods for elliptic problems on simplicial meshesSubjects: Numerical Analysis (math.NA)
We present both $hp$-a priori and $hp$-a posteriori error analysis of a mixed-order hybrid high-order (HHO) method to approximate second-order elliptic problems on simplicial meshes. Our main result on the $hp$-a priori error analysis is a $\frac12$-order $p$-suboptimal error estimate. This result is, to our knowledge, the first of this kind for hybrid nonconforming methods and matches the state-of-the-art for other nonconforming methods (as discontinuous Galerkin methods) with general (mixed Dirichlet/Neumann) boundary conditions. Our second main result is a residual-based $hp$-a posteriori upper error bound, comprising residual, normal flux jump, tangential jump, and stabilization estimators (plus data oscillation terms). The first three terms are $p$-optimal and only the latter is $\frac12$-order $p$-suboptimal. This result is, to our knowledge, the first $hp$-a posteriori error estimate for HHO methods. A novel approach based on the partition-of-unity provided by hat basis functions and on local Helmholtz decompositions on vertex stars is devised to estimate the nonconforming error. Finally, we establish local lower error bounds. Remarkably, the normal flux jump estimator is only $\frac12$-order $p$-suboptimal, as it can be bounded by the stabilization owing to the local conservation property of HHO methods. Numerical examples illustrate the theory.
- [20] arXiv:2503.11034 (replaced) [pdf, html, other]
-
Title: Inverse scattering for Schrödinger equation in the frequency domain via data-driven reduced order modelingSubjects: Numerical Analysis (math.NA)
In this paper we develop a numerical method for solving an inverse scattering problem of estimating the scattering potential in a Schrödinger equation from frequency domain measurements based on reduced order models (ROM). The ROM is a projection of Schrödinger operator onto a subspace spanned by its solution snapshots at certain wavenumbers. Provided the measurements are performed at these wavenumbers, the ROM can be constructed in a data-driven manner from the measurements on a surface surrounding the scatterers. Once the ROM is computed, the scattering potential can be estimated using non-linear optimization that minimizes the ROM misfit. Such an approach typically outperforms the conventional methods based on data misfit minimization. We develop two variants of ROM-based algorithms for inverse scattering and test them on a synthetic example in two spatial dimensions.
- [21] arXiv:2504.21603 (replaced) [pdf, html, other]
-
Title: Flow Through Porous Media: A Hopf-Cole Transformation Approach for Modeling Pressure-Dependent ViscositySubjects: Numerical Analysis (math.NA); Mathematical Physics (math-ph)
Most organic liquids exhibit a pressure-dependent viscosity, making it crucial to consider this behavior in applications where pressures significantly exceed ambient conditions (e.g., geological carbon sequestration). Mathematical models describing flow through porous media while accounting for viscosity-pressure dependence are nonlinear (e.g., the Barus model). This nonlinearity complicates mathematical analysis and makes numerical solutions more time-intensive and prone to convergence issues. In this paper, we demonstrate that the Hopf-Cole transformation, originally developed for Burgers' equation, can recast the governing equations -- describing flow through porous media with pressure-dependent viscosity -- into a linear form. The transformed equations, resembling Darcy's equations in the transformed variables, enable (a) systematic mathematical analysis to establish uniqueness and maximum principles, (b) the derivation of a mechanics-based principle, and (c) the development of efficient numerical solutions using solvers optimized for Darcy equations. Notably, many properties of the linear Darcy equations naturally extend to nonlinear models that depend on pressure. For example, solutions to these nonlinear models adhere to a reciprocal relation analogous to that observed in Darcy's equations.
- [22] arXiv:2507.14562 (replaced) [pdf, html, other]
-
Title: 1/2 order convergence rate of Euler-type methods for time-changed stochastic differential equations with super-linearly growing drift and diffusion coefficientsSubjects: Numerical Analysis (math.NA)
This paper investigates the convergence rates of two Euler-type methods for a class of time-changed stochastic differential equations with super-linearly growing drift and diffusion coefficients. Building upon existing research, we adapt the backward Euler method to time-changed stochastic differential equations where both coefficients exhibit super-linear growth and introduce an explicit counterpart, the projected Euler method. It is shown that both methods achieve the optimal strong convergence rate of order 1/2 in the mean-square sense for this class of equations. Numerical simulations confirm the theoretical findings
- [23] arXiv:2507.15103 (replaced) [pdf, html, other]
-
Title: Analysis of fully discrete Crank-Nicolson finite element methods for a stochastic Keller-Segel chemotaxis system with gradient-type multiplicative noiseSubjects: Numerical Analysis (math.NA)
We develop and analyze numerical methods for a stochastic Keller-Segel system perturbed by Stratonovich noise, which models chemotactic behavior under randomly fluctuating environmental conditions. The proposed fully discrete scheme couples a Crank-Nicolson time discretization with a splitting mixed finite element method in space. We rigorously prove the stability of the numerical scheme and establish strong convergence rates of order $O(k^{1/2} + k^{-1/2}h^2)$, where $k$ and $h$ denote the time and spatial step sizes, respectively. Notably, the presence of stochastic forcing leads to an inverse dependence on $k$ in the error estimates, distinguishing the convergence behavior from that of the deterministic case. Numerical experiments are presented to validate the theoretical results and demonstrate the effectiveness and accuracy of the proposed methods.
- [24] arXiv:2503.05674 (replaced) [pdf, html, other]
-
Title: Multiple solutions to the static forward free-boundary Grad-Shafranov problem on MAST-USubjects: Plasma Physics (physics.plasm-ph); Numerical Analysis (math.NA)
The Grad-Shafranov (GS) equation is a nonlinear elliptic partial differential equation that governs the ideal magnetohydrodynamic equilibrium of a tokamak plasma. Previous studies have demonstrated the existence of multiple solutions to the GS equation when solved in idealistic geometries with simplified plasma current density profiles and boundary conditions. Until now, the question of whether multiple equilibria might exist in real-world tokamak geometries with more complex current density profiles and integral free-boundary conditions (commonly used in production-level equilibrium codes) has remained unanswered. In this work, we discover multiple solutions to the static forward free-boundary GS problem in the MAST-U tokamak geometry using the validated evolutive equilibrium solver FreeGSNKE and the deflated continuation algorithm. By varying the plasma current, current density profile coefficients, or coil currents in the GS equation, we identify and characterise distinct equilibrium solutions, including both deeply and more shallowly confined plasma states. We suggest that the existence of even more equilibria is likely prohibited by the restrictive nature of the integral free-boundary condition, which globally couples poloidal fluxes on the computational boundary with those on the interior. We conclude by discussing the implications of these findings for wider equilibrium modelling and emphasise the need to explore whether multiple solutions are present in other equilibrium codes and tokamaks, as well as their potential impact on downstream simulations that rely on GS equilibria.