We present an algorithm for efficient evaluation of Boys functions F0,…,Fkmax tailored to modern computing architectures, in particular graphical processing units, where maximum throughput is high and data movement is costly. The method combines rational minimax approximations with upward and downward recurrence relations. The non-negative real axis is partitioned into three regions, [0, ∞⟩ = A ∪ B ∪ C, where regions A and B are treated using rational minimax approximations and region C by an asymptotic approximation. This formulation avoids lookup tables and irregular memory access, making it well-suited for hardware with high maximum throughput and low latency. The rational minimax coefficients are generated using the rational Remez algorithm. For a target maximum absolute error of ɛtol = 5 × 10-14, the corresponding approximation regions and coefficients for Boys functions F0, …, F32 are provided in AppendixD.

The Mittag-Leffler function holds significant importance in fractional calculus due to its extensive applications in addressing challenges across science, engineering, biology, hydrology, and earth sciences. Notably, the closed-form solution of a time-fractional model naturally emerges as the Mittag-Leffler function (MLF), necessitating precise and efficient computations. Consequently, numerical approximations are essential for accurately calculating the Mittag-Leffler function. In this study, we develop a straightforward yet precise real pole rational approximation for the Mittag-Leffler function. We demonstrate first-order convergence and L-acceptability, which aid in mitigating unwanted oscillations. Additionally, we create an effective and precise first-order generalized exponential time differencing scheme to solve the time-fractional reaction–diffusion equations. We obtain and prove the convergence result using Grönwall-type inequality. Several numerical experiments are conducted to confirm the efficiency and accuracy of the proposed numerical scheme compared with exact solutions. The computational efficiency of the proposed method is compared with another existing first-order numerical technique. Furthermore, our proposed scheme is crucial for developing higher-order predictor–corrector schemes for solving time-fractional models.

This paper presents a systematic methodology for time-domain modeling of three-phase power transformers aimed at electromagnetic transient analysis in shipboard and embedded electrical systems. Accurate modeling of transformers in such environments is critical, as naval power systems are subject to strict electromagnetic compatibility (EMC) requirements and are particularly susceptible to fast transients caused by switching operations, fault events, and nonlinear loads operating in confined and isolated grids. The proposed approach combines the Vector Fitting (VF) algorithm with Clarke modal decomposition to obtain stable, passive, and causal rational approximations of the frequency-dependent admittance matrix over a wide frequency range. The admittance matrix is first identified from frequency-domain measurements or simulations, capturing the transformer’s terminal behavior across multiple frequency sub-bands. Clarke’s transformation is then applied to decouple the three-phase system into independent modal components—namely the zero-sequence and positive-sequence modes, reducing the original multi-phase problem to a set of independent single-phase systems. This modal decoupling significantly improves computational efficiency without sacrificing accuracy, as each mode can be fitted and simulated independently.

Abstract We define and analyze preconditioners for the Riesz operator $$-(- \Delta )^{\frac{\alpha }{2}}$$ , $$\alpha \in (1,2]$$ commonly used in fractional models, such as anomalous diffusion. For $$\alpha$$ close to 2 there are various effective preconditioners at disposal with linear computational cost. Seminal results on treatment of the case $$\alpha$$ near 1 that still maintains linear computational complexity has been obtained approximating the Riesz operator as a fractional power of a discretized Laplacian, using Gauss-Jacobi formula. In this work, we extend this rational preconditioning approach by leveraging additional quadrature rules with exponential convergence. More precisely, we investigate both sinc and Gauss-Laguerre quadratures and show that, after an opportune choice of the involved parameters, both allow us to construct preconditioners based on a sum of a few shifted Laplacian inverses, and achieve high computational efficiency, ensuring numerical optimality. Several numerical results show that the sinc-based preconditioner is more versatile than the Gauss-Laguerre preconditioner, and that both outperform the Gauss-Jacobi one.

Abstract A classical problem in data-driven model order reduction (MOR) of linear time-invariant (LTI) systems is the preservation of structural properties of the underlying large-scale dynamics. When dealing with MOR based on transfer function measurements, one relevant problem is how to force the reduced-order model (ROM) to inherit the asymptotic stability of the reference system, i.e., to enforce the poles of the ROM transfer function to have strictly negative real part. In this work, we tackle the more general problem of placing such poles in arbitrary linear matrix inequality (LMI) regions of the complex plane, which include a rich class of convex sets symmetric with respect to the real axis. LTI systems with poles constrained to this kind of regions are called $$\mathcal {D}$$ D -stable. Combining well-established results from control theory and recent developments in asymptotically stable rational approximation algorithms, we show that the problem can be solved efficiently via standard convex optimization routines. Several numerical testbenches of engineering interest confirm the effectiveness of the proposed methodology in practical applications.

Topology optimization is a powerful computational tool for determining optimal material layouts within a design domain, aiming to minimize structural compliance under a prescribed volume constraint. Conventional density-based approaches, such as the solid isotropic material with penalization (SIMP) and the rational approximation of material properties (RAMP), are widely employed to interpolate material properties and drive the solution toward discrete "black-and-white" topologies. Although most topology optimization frameworks rely on the finite element method (FEM), alternative discretization strategies, such as the finite-volume theory (FVT), have demonstrated promising capabilities, particularly in addressing numerical issues like checkerboard patterns in the absence of filtering techniques. The proposed formulation employs compatibility and continuity conditions in the surface-averaged sense at the subvolume faces and locally satisfies the equilibrium equations at the subvolume level, providing a physically consistent basis for structural analysis. This work presents an extension of a density-based topology optimization algorithm grounded in the zeroth-order finite-volume theory for three-dimensional linear elastic structures. The numerical formulation adopts structured meshes and an artificial interpolation scheme to enhance solution discreteness and stability. In addition, an adaptation of the optimality criteria (OC) method is adopted, aimed at suppressing gray-scales, employing a linear relationship between material stiffness and density (Voigt model of micromechanics) instead of a penalization of intermediate densities approach. The formulation's effectiveness is demonstrated through numerical examples, establishing it as a promising alternative for 3D topology optimization.

Nuclear data consist of nuclide-specific files that tabulate a wide range of quantities, including reaction cross sections, scattering and absorption probabilities, and decay coefficients. Here we reinterpret this body of data through a graph-theoretic lens, representing it algebraically as an adjacency matrix. We then analyze how the structural and spectral properties of this matrix influence numerical solutions of the nuclear inventory equation, with particular emphasis on the Chebyshev rational approximation method and the backward differentiation formula methods. Activation-decay simulations of nuclear inventories typically proceed in two stages: an activation phase under irradiation, followed by a decay phase after the source is removed. While both phases solve the same stiff system of linear ordinary differential equations, the decay phase has a distinctive property: the transmutation graph becomes acyclic. We exploit this by applying a topological ordering of isotopes, which transforms the decay matrix into strictly triangular form. This removes the need for expensive LU factorizations during the decay phase, replacing it with a single forward or back-substitution step. The result is a gain in computational speed with increased accuracy.

<bold>Background</bold>The reactor burnup calculation is crucial for the safe operation and fuel management of nuclear power plants. In recent years, the Chebyshev Rational Approximation Method (CRAM) has become a primary approach for solving burnup equations. When solving the burnup equations using the CRAM, the Sparse Gaussian Elimination (SGE) method is usually used for complex matrix calculation, but the improvement of computational efficiency is limited.<bold>Purpose</bold>This study aims to develop a Gauss-Seidel (GS)-based acceleration method for solving the CRAM burnup equations to enhance the computational efficiency of the burnup equation solver.<bold>Methods</bold>Firstly, based on the self-developed burnup calculation code AMAC, an acceleration method for solving CRAM burnup equations was developed on the basis of the GS method. Then, three burnup databases (containing 71, 221, and 1 487 nuclides) were utilized to analyze the computational accuracy and efficiency of a light-water reactor benchmark, and the results of the SGE and GS methods were compared for evaluating computational accuracy, thereby demonstrating the precision of the GS method. Subsequently, detailed analyses of the computational results of the Partial Fraction Decomposition (PFD) and Incomplete Partial Fractions (IPF) formulations based on the GS method were conducted with the Transmutation Trajectory Analysis (TTA) as the reference solution. Finally, a comparative analysis of the computational efficiency of the SGE and GS methods was performed.<bold>Results</bold>The computational results show that the numerical precision of GS method is comparable for solving IPF and PFD burnup equations under different scale burnup databases. For the calculation of short-lived nuclides, the calculation accuracy of IPF is better than that of PFD. In terms of efficiency, the GS method significantly surpasses SGE method, achieving up to 80.17% acceleration across the three databases.<bold>Conclusions</bold>Results of this study recommend the adoption of the GS-accelerated IPF formalism for practical burnup calculations to effectively balance computational accuracy and efficiency.

Abstract In the quantum infinite square well with hard walls at x = 0 and L the energy eigenvalues are proportional to ϵ n 2 , where ϵ n = n is a dimensionless reduced wave number. We examine how the values of ϵ n are modified by the addition of one or more Dirac delta barriers, scaled so that the probability T for a particle to transmit through the barrier is not dependent on the particle’s energy. If the barriers are placed at rational locations x = iL / D , where D &gt; 1 is an integer common denominator and i is any integer such that 0 &lt; i &lt; D , then ϵ n follows a pattern that repeats every D values. Once the first D values of ϵ n are found, all other values follow from this pattern. Analytical solutions for ϵ n are derived for D = 2 and 3. Numerical solutions for D = 7 are shown to illustrate phenomena such as band structure, defects, and avoided crossings. When the barriers are placed at irrational locations there is no repeating pattern, but a finite sequence of ϵ values can be well-approximated by using results with the barriers placed at rational approximations of the irrational locations. Because this system admits analytical solutions, or numerical solutions obtained with minimal effort, it could serve as the basis for interesting exercises or projects in undergraduate or beginning graduate courses in quantum mechanics.

A bstract We revisit the solution to the Schwinger-Dyson equations in the simple case of the 0-dimensional $$ \frac{1}{2}{m}^2{\phi}^2+\frac{\lambda }{4}{\phi}^4 $$ 1 2 m 2 ϕ 2 + λ 4 ϕ 4 theory with m 2 &gt; 0 and λ ≥ 0. We argue that the truncated Schwinger-Dyson equations are solved by rational approximants to all n-point functions ⟨ ϕ 2 k ⟩, and provide strikingly simple recursive relations for them. These rational approximants are constructed without any reference to ordinary perturbative expansions. They turn out to be Padé approximants for ⟨ ϕ 2 ⟩ and for half of the truncations in the case of ⟨ ϕ 4 ⟩, but they are not Padé approximants for higher n-point functions. This difference is related to the fact that ⟨ ϕ 2 ⟩ and ⟨ ϕ 4 ⟩ are Stieltjes functions, while higher n-point functions are not. We prove that as the size of the truncation tends to infinity, these rational approximants converge to the full non-perturbative n-point functions for all positive values of the coupling λ . Thus, in the example studied in this work, these new rational approximants are much easier to derive than the usual Padé approximants, and when different, they are better suited to approximate the full non-perturbative n-point functions.

Simultaneous multiplicative rational approximation to a real and a p-adic numbers

The maximum principle is a widely used qualitative property of linear (and not only) elliptic boundary value problems. A natural goal for developing numerical methods is for the approximate solution to have a similar property. In this case, we say that a discrete maximum principle holds. In many cases, such a requirement is critical to ensuring the reliability of computational models. Here, we consider multidimensional linear elliptic problems with diffusion and reaction terms. Such problems have been studied and analyzed for many decades. Since relatively recently, scientists have faced conceptually new challenges when considering anomalous (fractional) diffusion. In the present paper, we concentrate on the case of spectral fractional diffusion. Discretization was carried out using the finite difference method and the finite element method with a lumped mass matrix. In large-scale multidimensional problems, the computational complexity of dense matrix operations is critical. To overcome this problem, BURA (best uniform rational approximation) methods were applied to find the efficient numerical solutions of emerging dense linear systems. Thus, along with the need to satisfy the discrete maximum principle associated with the mesh method applied for discretization of the differential operator, the issue of the monotonicity of BURA numerical solution arises. The presented results are three-fold and include the following: (i) maximum principles for fractional diffusion–reaction problems; (ii) sufficient conditions for discrete maximum principles; and (iii) sufficient conditions for monotonicity of the investigated BURA- or BURA-like approximation methods. A novel, systematic theoretical analysis is developed for sub-diffusion with a fractional power α∈(1/2,1) and a constant reaction coefficient. The theoretical findings are further supported by numerical examples.

Accurate evaluation of the standard normal cumulative distribution function is fundamental in many areas of mathematics, statistics, and applied computation, yet no closed-form expression in elementary functions exists. We present a simple analytic approximation based on a logistic function with a cubic argument, designed to preserve symmetry, monotonicity, and analytic invertibility. The parameters of the approximation are obtained through numerical optimization over a wide domain, targeting both maximum absolute error and root-mean-square error. The resulting function achieves uniformly low approximation error and significantly reduces error relative to the classical logistic approximation, while remaining competitive with commonly used high-accuracy numerical methods. Unlike rational or high-degree polynomial approximations, the proposed form admits an explicit inverse, making it convenient for applications requiring analytic quantile evaluation or inverse transform sampling. Numerical error analysis and illustrative examples demonstrate that the approximation provides a practical balance between accuracy, simplicity, and analytic tractability.

ABSTRACT We perform a comprehensive cosmographic analysis of the late-time Universe using the latest Dark Energy Spectroscopic Instrument (DESI) Data Release 2 (DR2) baryon acoustic oscillation (BAO) measurements, comparing Taylor, Padé, and Chebyshev expansions as model-independent reconstructions of the background expansion. We consider Padé approximants of order (2,1) and (2,2), a Chebyshev expansion, and a third-order Taylor series. Due to its limited radius of convergence, the Taylor expansion is constrained using only the low-redshift DESI sub-set ($z&amp;lt; 1$), while the rational Padé forms and the Chebyshev expansion are applied over the full DESI DR2 redshift range. Cosmographic parameters are inferred through a Bayesian Markov chain Monte Carlo (MCMC) analysis, and the resulting best-fitting reconstructions of $H(z)$, $d_L(z)$, and BAO distance indicators are compared with the predictions of the Lambda cold dark matter ($\Lambda$CDM) model. All methods are consistent with $\Lambda$CDM at low redshift, but the Chebyshev expansion exhibits noticeable deviations at higher redshifts, while the Padé$_{(2,1)}$ and Padé$_{(2,2)}$ reconstructions remain closely aligned with $\Lambda$CDM across the DESI DR2 range. A model-selection analysis based on Akaike Information Criterion and Bayesian Information Criterion shows a clear statistical preference for the Taylor expansion over low-z$\Lambda$CDM, and a strong preference for Padé cosmography over $\Lambda$CDM when the full DESI DR2 data set is used. These results demonstrate the constraining power of DESI DR2 for cosmographic studies and highlight the utility of rational approximants, especially Padé forms, in extending cosmography reliably to higher redshifts beyond the domain of traditional Taylor series.

Rational function approximation (RFA) techniques, including vector fitting (VF), are widely employed for modeling the frequency-dependent behavior of S-parameters. However, when applied to multi-port devices and wideband data, the performance of VF is often limited by reduced accuracy resulting from the processing of large S-parameter datasets. To address this issue, this paper proposes a gated rational complex neural network (GRCNN) for modeling multi-port S-parameters in passive microwave systems. Based on a rational complex neural network (RCNN) framework derived from rational function approximation, the GRCNN introduces a gating mechanism to enhance modeling adaptability. This mechanism adaptively selects the effective number of complex-conjugate and real pole–residue pairs. Model parameters are optimized end-to-end via backpropagation within a GPU-accelerated training framework. By incorporating causality constraints and a soft passivity penalty into the loss function, the GRCNN promotes causal behavior and reduces passivity violations on a sampled frequency grid, although strict global passivity over the continuous frequency domain is not guaranteed. Moreover, the models generated by the GRCNN can be directly converted into equivalent lumped circuits for macromodeling purposes. Comprehensive comparisons are conducted among VF, RCNN, weighted RCNN (WRCNN), GRCNN, and a gated complex-valued neural network (GCVNN). Experimental results on four test cases, ranging from 2 to 138 ports with wideband responses extending into the sub-THz range, demonstrate that the GRCNN achieves lower fitting error while maintaining causality and passivity.

We compute the constant of approximation for an arbitrary rational point on an arbitrary smooth cubic hypersurface X over a number field k, provided that there is a k-rational line somewhere on X.In the process, we verify the Coba conjecture for X.

Improving EMT simulations using frequency-shifted rational approximations

• A quadratically convergent rational integrator is developed for replicator dynamics on the probability simplex. • The method unconditionally preserves positivity, mass conservation, invariant faces, and both corner and internal equilibria. • A discrete analogue of the quotient rule governing component ratios evolution is established. • An embedded first-order auxiliary scheme enables adaptive time-stepping via efficient local error estimation. • Extensive numerical experiments confirm second-order accuracy and show superior performance compared to standard solvers in preserving qualitative dynamics. In this work, we introduce a quadratically convergent and dynamically consistent integrator specifically designed for the replicator dynamics. The proposed scheme combines a two-stage rational approximation with a normalization step to ensure confinement to the probability simplex and unconditional preservation of non-negativity, invariant sets and equilibria. A rigorous convergence analysis is provided to establish the scheme’s second-order accuracy, and an embedded auxiliary method is devised for adaptive time-stepping based on local error estimation. Furthermore, a discrete analogue of the quotient rule, which governs the evolution of component ratios, is shown to hold. Numerical experiments validate the theoretical results, illustrating the method’s ability to reproduce complex dynamics and to outperform well-established solvers in particularly challenging scenarios.

A new methodology for ultra-fast and accurate statistical EMT analysis in electric power-systems

Abstract Various computational problems as, e.g., equations with fractional diffusion operators, evaluation of high-dimensional integrals, the Møller–Plesset approach in quantum chemistry, etc., are easily solved by using approximations by rational functions or by exponential sums. In the case of Cauchy–Stieltjes or, respectively, Lebesgue–Stieltjes functions we provide a uniform proof of upper bounds of the convergence rates of their best approximations by rational functions or exponential sums. It turns out that the convergence rate by rational approximation is better than for exponential sums. We extend the analysis also to the approximation on infinite intervals and to the best approximation of the relative error. Instead of looking for the best approximation one can use the computationally cheaper quadrature method, in particular the sinc quadrature. The corresponding sharp error estimates are determined. The theoretical results are supported by numerical results.