*1.3. The Saxena* I *Function*

A Saxena *I* function *I m*,*n pk*,*qk*,1;*r*(*z*) is similarly defined as an ℵ function via a Mellin-Barnes type integral using integers *m*, *n*, *pk*, *qk* such that 0 ≤ *m* ≤ *qk*, 0 ≤ *n* ≤ *pk*, for *ai*, *bj*, *ai*,*k*, *bj*,*<sup>k</sup>* <sup>∈</sup> <sup>C</sup> with <sup>C</sup>, the set of complex numbers, and for *<sup>α</sup>i*, *<sup>β</sup>j*, *<sup>α</sup>i*,*k*, *<sup>β</sup>j*,*<sup>k</sup>* <sup>∈</sup> <sup>R</sup><sup>+</sup> = (0, <sup>∞</sup>) (*i* = 1, 2, . . . , *pk*; *j* = 1, 2, . . . , *qk*), and *k* = 1, . . . ,*r* [10], in the form

$$I\_{p\_k, q\_k, 1; r}^{m, n} \left( z \middle| \begin{array}{c} (a\_i, a\_i)\_{1, p} \\ (b\_j, \beta\_j)\_{1, q} \end{array} \right) = \frac{1}{2 \pi i} \int\_{\mathcal{C}} \mathcal{L}\_{p\_k, q\_k, 1; r}^{m, n} (s) \ z^{-s} ds,\tag{12}$$

with the Mellin representation of the kernel <sup>I</sup>*m*,*<sup>n</sup> pk*,*qk*,1;*r*(*s*)

$$\mathcal{L}\_{p\_k, q\_k, 1, r}^{m, n} \left( \begin{array}{c} (a\_i, a\_i)\_{1, p} \\ (b\_j, \beta\_j)\_{1, q} \end{array} \bigg| s \right) = \frac{\prod\_{j=1}^{m} \Gamma \big( b\_j + \beta\_j s \big) \prod\_{j=1}^{n} \Gamma \big( 1 - a\_j - a\_j s \big)}{\sum\_{k=1}^{r} \prod\_{j=m+1}^{q\_k} \Gamma \big( 1 - b\_{j,k} - \beta\_{j,k} s \big) \prod\_{j=n+1}^{p\_k} \Gamma \big( a\_{j,k} + a\_{j,k} s \big)}}. \tag{13}$$

Here

$$z^{-s} = \exp\left[-s\{\log|z| + i\arg z\}\right], \; z \neq 0, \; i = \sqrt{-1}. \tag{14}$$

The conditions for the contour C are the same as for the ℵ function.

#### **2. Approximations**

The approaches for approximating using Sinc functions are discussed in this section. First, the basic concepts of Sinc approximations are introduced describing the terms and notion. The second part deals with approximations of definite integrals. Based on these definitions we introduce the approximation of Mellin-Barnes integrals in the next step. We use the properties of Sinc functions allowing a stable and accurate approximation based on Sinc points [24]. For a detailed representation we refer the reader to [11,13].

#### *2.1. Sinc Basis*

To start with we first introduce some definitions and theorems allowing us to specify the space of functions, domains, and arcs for a Sinc approximation.

#### **Definition 1.** *Domain and Conditions.*

*Let* D *be a simply connected domain in the complex plane and z* ∈ C *having a boundary ∂*D*. Let a and b denote two distinct points of ∂*D *and φ denote a conformal map of* D *onto* D*d, where* <sup>D</sup>*<sup>d</sup>* <sup>=</sup> {*<sup>z</sup>* <sup>∈</sup> <sup>C</sup> : |I*m*(*z*)<sup>|</sup> <sup>&</sup>lt; *<sup>d</sup>*}*, such that <sup>φ</sup>*(*a*) = <sup>−</sup><sup>∞</sup> *and <sup>φ</sup>*(*b*) = <sup>∞</sup>*. Let <sup>ψ</sup>* <sup>=</sup> *<sup>φ</sup>*−<sup>1</sup> *denote the inverse conformal map, and let* Γ *be an arc defined by* Γ = {*z* ∈ C : *z* = *ψ*(*x*), *x* ∈ R}*. Given φ, ψ, and a positive number h, let us set zk* = *ψ*(*kh*)*, k* ∈ Z *to be the Sinc points, let us also define ρ*(*z*) = *eφ*(*z*)*.*

Note the Sinc points are an optimal choice of approximation points in the sense of Lebesgue measures for Sinc approximations [24].

#### **Definition 2.** *Function Space.*

*Let d* ∈ (0, *π*)*, and let the domains* D *and* D*<sup>d</sup> be given as in Definition 1. If d is a number such that d* > *d, and if the function φ provides a conformal map of* D *onto* D*d , then* D⊂D *. Let α and β denote positive numbers, and let Lα*,*β*(D) *denote the family of functions u* ∈ *Hol* (D)*, for which there exists a positive constant c*<sup>1</sup> *such that, for all z* ∈ D*. Let Lα*,*β*(D) *be the set of all analytic functions, for which there exists a constant c*1*, such that*

$$|u(z)| \le c\_1 \frac{|\rho(z)|^a}{(1 + |\rho(z)|)^{a + \beta}}.\tag{15}$$

*Now let the positive numbers α and β belong to* (0, 1]*, and let Mα*,*β*(D) *denote the family of all functions g* ∈ *Hol* (D)*, such that g*(*a*) *and g*(*b*) *are finite numbers, where g*(*a*) = lim*z*→*<sup>a</sup> g*(*z*) *and g*(*b*) = lim*z*→*<sup>b</sup> <sup>g</sup>*(*z*)*, and such that u* ∈ *<sup>L</sup>α*,*β*(D) *where*

$$
\mu(z) = g(z) - \frac{g(a) + \rho(z)g(b)}{1 + \rho(z)}.\tag{16}
$$

These definitions directly allow the formulation of an algorithm for a Sinc approximation. Let Z denote the set of all integers. Select positive integers *N* and *M* = [*βN*/*α*] so that *m* = *M* + *N* + 1. The step length is determined by *h* = (*πd*/(*βN*))1/2 where *α*, *β*, and *d* are real parameters. In addition assume there is a conformal map *φ* and its inverse *ψ* such that we can define Sinc points *zj* = *ψ*(*jh*), *j* ∈ Z [25]. The following relations define the basis of a Sinc approximation:

$$\text{Sinc}\left(z\right) = \frac{\sin\left(\pi z\right)}{\pi z}.\tag{17}$$

The shifted Sinc is derived from relation (17) by translating the argument by integer steps of length *h* and applying the conformal map to the independent variable

$$S(k,h)\circ(z) = \sin(\pi(z/h-k))/(\pi(z/h-k)) = \text{Sinc}\left(z/h-k\right).\tag{18}$$

The approximation of a function *f*(*z*) results to the representation

$$\mathbb{C}\_{h,M,N}[f](z) = \sum\_{k=-M}^{N} f(z\_k) \text{Sinc}\left(\frac{\phi(z)}{h} - k\right) = \sum\_{k=-M}^{N} f(z\_k) \text{S}(k,h) \circ \phi(z),\tag{19}$$

using the set of orthogonal functions

$$B\_S = \{ S(k, h) \circ \phi(z) \}\_{k=-M'}^{N} \tag{20}$$

where *φ*(*z*) is the conformal map. This type of approximation allows to represent a function *f*(*z*) on an arc Γ with an exponential decaying accuracy [11]. As proved in [11,13] the approximation works effectively for analytic functions. The approximations (19) allow us to formulate the following theorem for Sinc approximations.

#### **Theorem 1.** *Sinc Approximation [25].*

*Let u* ∈ *Lα*,*β*(D) *for α* > 0 *and β* > 0*, take M=[β N/α], where [x] denotes the greatest integer in x, and then set <sup>m</sup>* <sup>=</sup> *<sup>M</sup>* <sup>+</sup> *<sup>N</sup>* <sup>+</sup> <sup>1</sup>*. If <sup>u</sup>* <sup>∈</sup> *<sup>M</sup>α*,*β*(D)*, and if <sup>h</sup>* = (*πd*/(*βN*))1/2 *then there exist positive constants K*<sup>1</sup> *and k*<sup>1</sup> *independent of N, such that*

$$\epsilon\_N = \left\| f(z) - \mathbb{C}\_{h,M,N}[f](z) \right\| = K\_1 N^{1/2} \exp\left( -k\_1 N^{1/2} \right). \tag{21}$$

*with wk the base function (see (22)).*

The proof of this theorem is given in [13]. Note the choice *h* = (*πd*/(*βN*))1/2 is close to optimal for an approximation in the space *Mα*,*β*(D) in the sense that the error bound in Theorem 1 cannot be appreciably improved regardless of the basis [25]. It is also optimal in the sense of the Lebesgue measure achieving an optimal value less than Chebyshev approximations [24].

Here *zk* = *ψ*(*kh*) = *φ*−1(*kh*) are the discrete points based on Sinc points *kh*. Note that the discrete shifting allows us to cover the approximation interval (*a*, *b*) in a dense way while the conformal map is used to map the interval of approximation from an infinite range of values to a finite one. Using the Sinc basis we are able to represent the basis functions as a piecewise defined function *wj*(*z*) by

$$w\_{j} = \begin{cases} \frac{1}{1+\rho(z)} - \sum\_{k=-M+1}^{N} \frac{1}{1+e^{kh}} S(k,h) \diamond \phi(z) & j=-M\\ S(k,h) \diamond \phi(z) & j=-M+1, \dots, N-1\\ \frac{\rho(z)}{1+\rho(z)} - \sum\_{k=-M}^{N-1} \frac{e^{kh}}{1+e^{kh}} S(k,h) \diamond \phi(z) & j=N \end{cases},\tag{22}$$

where *ρ*(*z*) = exp(*φ*(*z*)). This form of the Sinc basis is chosen as to satisfy the interpolation at the boundaries. The basis functions defined in (22) suffice for purposes of uniform−norm approximation over (*a*, *b*).

This notation allows us to define a row vector *Vm*(*S*) of basis functions

$$\mathbf{V}\_{m}(\mathbf{S}) = (w\_{-M\prime}, \dots, w\_{N}) \tag{23}$$

with *wj* defined as in (22). For a given vector *Vm*(*u*) = (*u*−*M*,..., *uN*)*<sup>T</sup>* we now introduce the dot product as an approximation of the function *u*(*z*) by

$$
\mu(z) \approx V\_{\mathfrak{m}}(S). \mathbf{V}\_{\mathfrak{m}}(\mu) = \sum\_{k=-M}^{N} u\_k w\_k. \tag{24}
$$

Based on this notation, we will introduce in the next subsection approximations of definite integrals [13].

#### *2.2. Definite Integral Approximation*

In this section, we pose the query of how to approximate definite integrals on a domain over R. The approximation will use our basis system introduced in Section 2.1. It turns out that for the basis systems (20), we can get an approximation converging exponentially. Specifically, we are interested in a quadrature formula for definite integrals of the type

$$J(f) = \int\_{a}^{b} f(t) \, dt,\tag{25}$$

where *a* and *b* can be finite or infinite. If the function *f* is approximated by the approximation given in Section 2.1, we write for *J*(*f*),

$$J(f) \quad = \int\_{a}^{b} f(t) \, dt \approx$$

$$J\_{h,M,N}[f](\mathbf{x}) \quad = \int\_{a}^{b} \sum\_{k=-M}^{N} f(t\_k) S(k,h) \circ \phi(t) \, dt \tag{26}$$

$$\xi = \sum\_{k=-M}^{N} f(t\_k) \int\_a^b S(k, h) \diamond \phi(t) dt. \tag{27}$$

Scaling the variable *ξ* = *t*/*h* and collocating the expression with respect to *ξ*, we end up with the representation

$$J(f) \approx h \sum\_{k=-M}^{N} f(t\_k) \frac{1}{\phi'(t\_k)}, \text{ with } t\_k = \psi(kh). \tag{28}$$

Equation (28) represents a quadrature formula applicable to different domains. Compared with a Gaussian quadrature formula (28) delivers the weights as well as the function values if the discretization is known by *tk* = *ψ*(*kh*) with *k* = −*M*, ... , *N* and *h* = *π* √*<sup>N</sup>* . Note since the conformal map *φ* depends on the structure of the arc Γ; i.e., finite, semiinfinite or infinite, the approximation is defined for domains [*a*, *b*], (0, ∞), or (−∞, ∞), respectively [11]. The following Theorem summarizes these results.

#### **Theorem 2.** *Definite Integrals.*

*If φ denotes a one*−*to*−*one transformation of the interval* (*a*, *b*) *onto the real line* R*, let h denote a fixed positive number, and let the Sinc points be defined on* (*a*, *<sup>b</sup>*) *by zk* <sup>=</sup> *<sup>φ</sup>*−1(*kh*)*, <sup>k</sup>* <sup>∈</sup> <sup>Z</sup>*, where φ*−<sup>1</sup> = *ψ denotes the inverse function of the conformal map φ. Let M and N be positive integers, set m* = *M* + *N* + 1*, and for a given function f defined on* (*a*, *b*)*, define the vector Vm*(*f*) = (*f*(*z*−*M*),..., *f*(*zN*))*, and a vector Vm*(1/*φ* ) = (1/*φ* (*z*−*M*), . . . , 1/*φ* (*zN*))*, then the definite integral is approximated by*

$$J(f) \approx h \Psi\_m(f). \Psi\_m(1/\phi') = f\_{h,M,N}[f](z),\tag{29}$$

*and the error of this approximation was estimated in [11] as*

$$\epsilon\_N = \left\| \left[ f(f) - f\_{\hbar, M, N}[f](z) \right] \right\| \sim K\_2 N^{1/2} \exp\left( -k\_2 N^{1/2} \right),\tag{30}$$

*where K*<sup>2</sup> *and k*<sup>2</sup> *are constants independent of N.*

#### *2.3. Sinc Approximation of Mellin-Barnes Integrals*

The Sinc approximation of the Fox *H* functions needs two discretization steps. Foremost the Sinc discretization of the arc Γ on which the function should be represented finite or semi-infinite and second the Sinc quadrature of the Mellin-Barnes integral at these Sinc points on the arc. The second discretization is a Sinc quadrature on an infinite interval corresponding to the line at *γ* parallel to the imaginary axis in *s* ∈ C. In formulas this means

$$\begin{split} \left( \mathcal{H}\_{p,\eta}^{m,\mathfrak{n}} \begin{pmatrix} \begin{matrix} a\_{i},a\_{i} \\ b\_{\boldsymbol{\upbeta}},\boldsymbol{\upbeta}\_{\boldsymbol{\upbeta}} \end{matrix} \end{pmatrix}\_{1,\boldsymbol{\uprho}} \right)\_{1,\boldsymbol{\uprho}} &= \begin{array}{c} \frac{1}{2\pi} \int\_{-\infty}^{\infty} \mathcal{H}\_{p,\boldsymbol{\uprho}}^{m,\mathfrak{n}} (\boldsymbol{\upgamma} + i\boldsymbol{\uprho}) \, z\_{k}^{-\gamma-i\sigma} d\sigma \\ & \approx \begin{array}{c} \frac{\mathrm{h}}{2\pi} \mathbf{V}\_{\mathfrak{m}} \left( \mathcal{H}\_{p,\boldsymbol{\uprho}}^{m,\mathfrak{n}} (\boldsymbol{\upgamma} + i\boldsymbol{\uprho}) z\_{k}^{-\gamma-i\sigma} \right) . \mathbf{V}\_{\mathfrak{m}} (1/\boldsymbol{\uprho}'(\boldsymbol{\upsigma})) , \end{array} \right) \end{split} \tag{31}$$

where *σ<sup>l</sup>* = *ψ*(*lh*), *l* = −*M*, ... , *N*, and *zk* = *ψ*(*kh*), *k* = −*M*, ... , *N*, generating a vector *Vm Hm*,*<sup>n</sup> <sup>p</sup>*,*<sup>q</sup>* which approximates the Fox *H* function pointwise at *zk* = *ψ*(*kh*) on an arc Γ. *γ* is selected according to the pole structure of the Γ functions of the numerator in the Mellin representation. Using the basis functions *wk* (22) on the arc Γ, we shall approximate the Fox H function by

$$H\_{p, \boldsymbol{\eta}}^{m, n}(z) \approx V\_m(\boldsymbol{\mathcal{S}}) \mathbf{V}\_m\left(H\_{p, \boldsymbol{\mathcal{q}}}^{m, n}\right) = \sum\_{k=-M}^{N} H\_{p, \boldsymbol{\eta}}^{m, n}(z\_k) w\_k(z). \tag{32}$$

This approach guarantees that the error  *<sup>N</sup>* will decay as given by (21). Note that a Sinc point based interpolation is able to deal with singularities at the end points of the arc Γ.

Since the ℵ function and the Saxena *I* functions are similarly defined by Mellin-Barnes integrals as the Fox *H* function, the procedure for the approximation follows the same line delivering for the ℵ function

$$
\aleph\_{p\_k, q\_k, \tau\_k \mathcal{I}}^{m, n}(z) \approx V\_m(S). \\
\mathbf{V}\_m \left( \aleph\_{p\_k, q\_k, \tau\_k \mathcal{I}}^{m, n} \right) = \sum\_{k=-M}^{N} \aleph\_{p\_k, q\_k, \tau\_k \mathcal{I}}^{m, n}(z\_k) w\_k(z), \tag{33}
$$

and the Saxena *I* function

$$I\_{p\_k, q\_k, 1; r}^{\mathfrak{m}, \mathfrak{n}}(z) \approx V\_{\mathfrak{m}}(S). \mathcal{V}\_{\mathfrak{m}}\left(I\_{p\_k, q\_k, 1; r}^{\mathfrak{m}, \mathfrak{n}}\right) = \sum\_{k=-M}^{N} I\_{p\_k, q\_k, 1; r}^{\mathfrak{m}, \mathfrak{n}}(z\_k) w\_k(z). \tag{34}$$

All three approximations will satisfy the *a priory* error formula (21) and converge exponentially to their exact value.

The following section shall demonstrate by a few examples the numerical representation of the three generalized functions.

#### **3. Numerical Examples**

This section collects some examples demonstrating the efficient and accurate numerical evaluation of Sinc-based methods applied to Mellin-Barnes integrals, Fox H-, Saxena I-, and ℵ functions. We selected the examples concerning their analytic representation applied to specific physical or engineering problems. Allowing us to compare our results with exact expressions and thus delivering an *a priori* estimation of the numerical errors. For every comparison in which we calculated local errors, the analytical solution was employed. The following graphics display the analytical answer as a solid line. Over the solid line, the equivalent numerical approximation is displayed as a dashed curve. We may depend on the previously mentioned exponential convergence of the Sinc quadrature and Sinc approximation if no local errors are calculated in the graphs. The analytic Sinc approximations, however, are typically not represented because of their large symbolic representation since we are utilizing a computer algebra system.

#### *3.1. Fox H Functions*

In this subsection, we present some examples which are taken from the literature where the symbolic representation was given [1,2,4]. In most of these cases a numerical evaluation of the functions is not offered or discussed. For the first time we will show the numerical evaluation and the *a priori* error estimation using Sinc methods.

#### **Example 1.** *Exponential function*

*It is well-known that the exponential function is connected with the Mellin transform via the Euler* Γ *function. The relations between the two functions are as follows*

$$e^{-\mathbf{x}} = \frac{1}{2\pi i} \int\_{\gamma - i\infty}^{\gamma + i\infty} \Gamma(s) \mathbf{x}^{-s} ds,\ \gamma > 0,\tag{35}$$

*and*

$$
\Gamma(s) = \int\_0^\infty e^{-\mathbf{x}} \mathbf{x}^{s-1} d\mathbf{x}\_\prime \tag{36}
$$

*where the later of these two expressions is the Mellin transform of the exponential function and the first its inverse given by a Bromwich integral over the vertical line γ* = *const , of the complex s-plane. Whereas the last formula (36) is due to L. Euler who communicated it in two letters of 13 October 1729, and 8 January 1730, to Goldbach. The first formula (35) in a somewhat different form goes back to Pincherle [9,26] as Mellin reported in his paper [27]; but it was Mellin who first realized its great importance.*

*Let us first numerically examine the representation of the exponential function using the Fox H representation. The exponential function is defined in terms of a Fox H function by*

$$e^{-\mathbf{x}} = H\_{0,2}^{2,0} \left( \begin{matrix} \mathbf{x} \\ \mathbf{2} \end{matrix} \; \middle| \; \begin{matrix} \mathbf{0}\_{\prime} \frac{1}{2} \end{matrix} \; \begin{matrix} \mathbf{1}\_{\prime} \end{matrix} \; \middle| \; \begin{matrix} \frac{1}{2}, \frac{1}{2} \end{matrix} \right) , \text{ with } \; \mathcal{H}\_{0,2}^{2,0}(\mathbf{s}) = \frac{\Gamma\left(\frac{\mathbf{s}}{2} + \frac{1}{2}\right) \Gamma\left(\frac{\mathbf{s}}{2}\right)}{2\sqrt{\pi}}. \tag{37}$$

*Using the duplication formula for the* Γ *function in the right most term of (37), we end up with the Mellin representation in (35). The line integral is then calculated by the Mellin-Barnes integral via*

$$\begin{split} e^{-z} &=& \frac{1}{2\pi i} \int\_{\mathcal{C}} \mathcal{H}\_{0,2}^{2,0}(s) \, z^{-s} ds = \frac{1}{2\pi i} \int\_{\gamma-i\infty}^{\gamma+i\infty} \mathcal{H}\_{0,2}^{2,0}(s) \, z^{-s} ds \\ &=& \frac{1}{2\pi i} \int\_{-\infty}^{\infty} \mathcal{H}\_{0,2}^{2,0}(\gamma-i\sigma) \, z^{-\gamma+i\sigma} d\sigma,\end{split} \tag{38}$$

*where the path* C *is chosen according to the pole structure of the* Γ *functions so that γ is a fixed real number. For different z values selected as Sinc points, we compute the numerical values of the integral in the domain for <sup>σ</sup>* <sup>∈</sup> (−∞, <sup>∞</sup>)*. The set of values* {*zk*, *Ik*}*<sup>N</sup> <sup>k</sup>*=−*<sup>N</sup> is used in a Sinc approximation delivering an analytic representation of the function f*(*z*) = *e*−*<sup>z</sup> in terms of the Sinc basis. Results of such an approximation are illustrated in the following Figure 1. Figure 1 shows the exact function (solid line), the Sinc approximation (dashed), and the discrete set of Sinc points (dots). The right panel shows the local error of the approximation. In Figure 2, we depict the error decay of the Sinc approximation as a function of the used number of Sinc points N. A least square fit was used to get the two parameters K*<sup>2</sup> *and k*<sup>2</sup> *which determine the error formula. The example demonstrates that the approximation of the Fox H function is accurate and the precision of the approximation can be a priori estimated using the error decay Formula (30).*

**Figure 1.** Numerical representation of the Fox H function *H*2,0 0,2 (*x*) = *<sup>e</sup>*−*<sup>x</sup>* as a Sinc approximation based on (35) using *N* = 24 Sinc points for *x*. The *L*<sup>2</sup> norm over the local error delivers a value of 1.784 10−6.

**Figure 2.** Error decay  *<sup>N</sup>* <sup>=</sup> *<sup>u</sup>* <sup>−</sup> *<sup>u</sup>*ex <sup>∼</sup> *<sup>K</sup>*2*N*1/2 exp −*k*2*N*1/2 with *K*<sup>2</sup> = 2.623 and *k*<sup>2</sup> = 3.017. The two parameters *K*<sup>2</sup> and *k*<sup>2</sup> are the result of a least square fit.

#### **Example 2.** *Mittag-Leffler functions*

*Over the years it turned out that Mittag-Leffler (ML) functions are of central importance in fractional calculus [5,8,14,28]. This kind of functions are used in many applications and in theoretical work [2]. However, these functions need a special approach when it comes to their numerical representation. Quite recently, we presented an approach to tackle the numerical representation by an inverse Laplace transform [7]. It was shown that a Sinc based approach has some advantages over the standard Talbot or Weideman approach used by Garrappa et al. [6,29,30]. A second approach to numerically represent Mittag-Leffler functions is now available with high accuracy and precision using Fox H functions in connection with Sinc quadrature. We will restrict our discussions to the main three types of Mittag-Leffler functions which are classified as single-, two-, and three-parameter ML functions. The three-parameter ML function is also known as Prabhakar's function. There are in fact higher parameter ML functions in use which are defined via Fox H functions [28]. We note that these functions are also numerically available with high accuracy. However, due to lack of space, we will not present them here.*

*The three ML functions we have in mind are defined by the following relations:*

$$E\_a(z) = H\_{1,2}^{1,1} \left( -z \middle| \begin{array}{c} (0,1) \\ (0,1) & (0,a) \end{array} \right), \text{ with } \mathcal{H}\_{1,2}^{1,1}(s) = \frac{\Gamma(1-s)\Gamma(s)}{\Gamma(1-sa)},\tag{39}$$

$$E\_{\mathfrak{a},\mathfrak{\mathfrak{\mathfrak{F}}}}(z) = H\_{1,2}^{1,1} \left( -z \Big| \begin{array}{c} (0,1) \\ (0,1) & (1-\beta,\mathfrak{a}) \end{array} \right), \text{ with } \mathcal{H}\_{1,2}^{1,1}(s) = \frac{\Gamma(1-s)\Gamma(s)}{\Gamma(\beta-s\mathfrak{a})},\tag{40}$$

*and*

$$E\_{a,\beta}^{\mu}(z) = H\_{1,2}^{1,1}\left(-z \middle| \begin{array}{c} (1-\mu, 1) \\ (0,1) \end{array} \right. \\ \left.(1-\beta, a) \right. \right), \text{ with } \mathcal{H}\_{1,2}^{1,1}(s) = \frac{\Gamma(\mu-s)\Gamma(s)}{\Gamma(\beta-sa)}.\tag{41}$$

*The line integral is computed by the Mellin-Barnes integral via*

$$\begin{split} I &=& \frac{1}{2\pi i} \int\_{\mathcal{C}} \mathcal{H}^{1,1}\_{1,2}(s) \, z^{-s} ds = \frac{1}{2\pi i} \int\_{\gamma - i\infty}^{\gamma + i\infty} \mathcal{H}^{1,1}\_{1,2}(s) \, z^{-s} ds \\ &=& \frac{1}{2\pi} \int\_{-\infty}^{\infty} \mathcal{H}^{1,1}\_{1,2}(\gamma - i\sigma) \, z^{-\gamma + i\sigma} d\sigma, \end{split} \tag{42}$$

*where the path* C *is selected according to the pole structure of the* Γ *functions so that γ is a fixed real number. For different z values selected as Sinc points, we compute the numerical values of the integral in the domain for <sup>σ</sup>* <sup>∈</sup> (−∞, <sup>∞</sup>)*. The set of values* {*zk*, *Ik*}*<sup>N</sup> <sup>k</sup>*=−*<sup>N</sup> is used in a Sinc approximation.*

*In Figure 3 we graph the solution for the one-, two-, and three-parameter ML function for a specific selection of parameters (left column of Figure 3). The dashed lines correspond to the Sinc approximation while the solid line represents the Mathematica approximation for the first two ML functions and a higher order Sinc approximation with N* = 128 *for the Prabhakar ML function. Obviously there is no visible difference for all three types of ML functions. This becomes obvious if we examine the right column of Figure 3. Here the absolute value of the local error is shown for the two first ML functions and a relative error for Prabhakar's function. The magnitude in all cases is less than* 10−10*. Figure 4 collects the error decay as a function of the number of Sinc points N used in the approximations. It is obvious that for all three types of ML functions nearly the same functional error decay results. This indicates that the upper bound estimation of Equation (30) is satisfied and allows an a priori estimation of the expected error.*

**Figure 3.** Function representation for *<sup>E</sup>α*(*x*), *<sup>E</sup>α*,*β*(*x*), and *<sup>E</sup><sup>μ</sup> <sup>α</sup>*,*β*(*x*) from top to bottom, left column with *N* = 56 Sinc points. The parameters *α*, *β*, and *μ* are indicated in the axis labels, respectively. The local errors are shown in the right column accordingly. For Prabhakar's function we have to use a relative error because the exact representation is not available. The reference solution *u*(*x*) was computed with *N* = 128 Sinc points.

**Figure 4.** Error decay  *<sup>N</sup>* <sup>=</sup> *<sup>u</sup>* <sup>−</sup> *<sup>u</sup>*ex <sup>∼</sup> *<sup>K</sup>*2*N*1/2 exp −*k*2*N*1/2 for the three ML functions *Eα*, *Eα*,*β*, and *E<sup>μ</sup> <sup>α</sup>*,*<sup>β</sup>* (dots, solid), (diamond, dot dashed), and (square, dashed), respectively. The least square fit for each function delivered the values (*K*2, *k*2)=(0.421, 3.288), (*K*2, *k*2)=(0.347, 3.292), and (*K*2, *k*2)=(0.126, 3.123).

*The least square parameters of the three error decays also indicates that all three functions belong to the same function space introduced in Definition 2.*

#### **Example 3.** *Krätzel Function*

*The Krätzel integral Z<sup>ν</sup> <sup>ρ</sup>* (*z*) *was defined by Krätzel [31] as the kernel of an integral transform as follows:*

$$Z^{\nu}\_{\rho}(z) = \int\_0^{\infty} y^{\nu - 1} e^{-y^{\rho} - x/y} dy. \tag{43}$$

*Today we know that the Krätzel function occurs in many fields of applications, such as the study of astrophysical thermonuclear reactions, reaction rate probability integrals in the theory of nuclear reaction rates, in applied analysis, inverse Gaussian distribution, generalized families of distributions in statistical distribution theory, and in statistical mechanics as well as the general pathway model are all shown to be connected to the integral (43) [4]. The generalized Krätzel function was examined by Kilbas and Kumar [32]. Solar radiation data were examined recently in connection with probability distribution functions by Princy [33]. Due to the versatile applications of the Krätzel function, we demonstrate the numerical representation as a density function.*

*The normalized probability density can be represented as Fox H function by*

$$f(z) = \ c z^{a-1} Z\_{\rho}^{\upsilon}(z) = \frac{1}{\Gamma(a) \Gamma\left(\frac{a+\upsilon}{\rho}\right)} H\_{0,2}^{2,0} \left( z \middle| \begin{array}{c} (a-1,1) \ \end{array} \left( \frac{a+\upsilon-1}{\rho}, \frac{1}{\rho} \right) \right), \text{with} \\\mathcal{H}\_{0,2}^{2,0}(s) = \Gamma(s+a-1) \Gamma\left(\frac{s}{\rho} + \frac{a+\upsilon-1}{\rho}\right) \end{array} \tag{44}$$

*here x* ≥ 0*, α* > 0*, ρ* > 0*, and ν* > 0*. The corresponding line integral representation is computed by the Mellin-Barnes integral via*

$$\begin{split} I &=& \frac{1}{2\pi i} \int\_{\mathcal{C}} \mathcal{H}\_{0,2}^{2,0}(s) \, z^{-s} ds = \frac{1}{2\pi i} \int\_{\gamma - i\infty}^{\gamma + i\infty} \mathcal{H}\_{0,2}^{2,0}(s) \, z^{-s} ds \\ &=& \frac{1}{2\pi} \int\_{-\infty}^{\infty} \mathcal{H}\_{0,2}^{2,0}(\gamma - i\sigma) \, z^{-\gamma + i\sigma} d\sigma, \end{split} \tag{45}$$

*where the path* C *is selected according to the pole structure of the* Γ *functions, so that γ is a fixed real number. For different discrete z values using Sinc points, we compute the numerical values of the integral in the real domain <sup>σ</sup>* <sup>∈</sup> (−∞, <sup>∞</sup>)*. The set of values* {*zk*, *Ik*}*<sup>N</sup> <sup>k</sup>*=−*<sup>N</sup> is used in a Sinc approximation. An example using specific parameters α* = 2 *and ρ* = *ν* = 3 *is shown in Figure 5. The left panel of Figure 5 represents the probability density based on N* = 128 *and N* = 56 *(dashed) Sinc points. The right panel of Figure 5 shows the relative local error of these two versions of the approximation. It is obvious that the local error is small and nearly homogenous on the domain* [0, 10]*. In Figure 6, we present the error decay of computations with different number of Sinc points zk used in the approximation. The solid line represents the least square approximation of the two parameter error formula (30). Note that for small numbers of Sinc points we already reach an acceptable small error. If we double the number of Sinc points, we gain more than one decade of accuracy. Figure 6 clearly demonstrates that we have an exponential decay of the error.*

**Figure 5.** The Krätzel function used as a probability density function for *ρ* = 3, *ν* = 2, and *α* = 2 (left panel). The relative local error of the distribution *f*(*x*) = *xα*−1*Z<sup>ν</sup> <sup>ρ</sup>* (*x*) (right panel). Number of Sinc points *N* = 56.

**Figure 6.** Error decay  *<sup>N</sup>* <sup>=</sup> *<sup>u</sup>* <sup>−</sup> *<sup>u</sup>*ex <sup>∼</sup> *<sup>K</sup>*2*N*1/2 exp −*k*2*N*1/2 for the Krätzel function with *ρ* = *ν* = 3, *α* = 2. The least square fit delivered the values (*K*2, *k*2)=(5.713, 2.807).

*In Figure 7 we demonstrate the variation of the Krätzel function if we change one of the parameters α or ν, respectively.*

**Figure 7.** The Krätzel function for *ρ* = *ν* = 3 with *α* variation (left graph) and *ν* variation (right graph).

*For the special choice ρ* = 1*, the Krätzel density is reduced to the modified Bessel function Kν*(*x*) *with the specific representation*

$$f(z) = cz^{a-1}Z\_{\rho}^{\upsilon}(z) = \frac{2}{\Gamma(a)\Gamma(a+\upsilon)}x^{a+\upsilon/2-1}K\_{\upsilon}(2\sqrt{x}), \text{with } \ge 0, a>0, \upsilon>0. \tag{46}$$

*This relation can be used to check the accuracy of our numerical computations shown in Figure 8. For a few number of approximation points N* = 56*, we are able to reach a low level of local errors.*

**Figure 8.** The Krätzel function for *ρ* = 1, *ν* = 1/3 and *α* = 2 (left graph). Right panel local error between the approximation and the Bessel representation. Number of Sinc points *N* = 56.

x

**Example 4.** *Abel Equation*

*In the current example, we will examine the following Abel integral equation*

$$u(\mathbf{x}) = f(\mathbf{x}) + \frac{\lambda}{\Gamma(\alpha)} \int\_0^\mathbf{x} (\mathbf{x} - \boldsymbol{\xi})^{\alpha - 1} u(\boldsymbol{\xi}) d\boldsymbol{\xi},\tag{47}$$

*where λ and* 0 < *α* < 1 *are real parameters [8,34]. The integral can also be identified as a Riemann-Liouville fractional integral so that (47) can be written as*

$$
\mu(\mathbf{x}) = f(\mathbf{x}) + \lambda \mathcal{D}\_{0,\mathbf{x}}^{-\mathbf{a}} \mu(\mathbf{x}).\tag{48}
$$

*This in turn opens the connection to the fractional kinetic equation for astrophysical systems examined by Haubold and Mathai [35] using a constant function <sup>f</sup>*(*x*) = *<sup>N</sup>*<sup>0</sup> *and <sup>λ</sup>* <sup>=</sup> <sup>−</sup>*cα. However, such kind of equation was already examined in 1991 by Glöckle and Nonnemacher in connection with viscoelastic materials [36]. We will numerically demonstrate that the use of Fox H functions in connection with a Sinc convolution integral allows a general solution as long as the Laplace transform of f*(*x*) *exists. In their work on Mittag-Leffler functions, Gorenflow et al. [8] state in Theorem 4.2 the problem based on [37] in which the solution of (47) using convolution integrals are examined. The solution of (47) using two parameter Mittag-Leffler functions [35] reads*

$$
\mu(\mathbf{x}) = f(\mathbf{x}) + \lambda \int\_0^\mathbf{x} (\mathbf{x} - \boldsymbol{\xi})^{a-1} E\_{\mathbf{a},a} (\lambda (\mathbf{x} - \boldsymbol{\xi})^a) f(\boldsymbol{\xi}) d\boldsymbol{\xi}.\tag{49}
$$

*This result goes back to the pioneering work of Hille and Tamarkin in 1930 [37]. Introducing new variables in (49) by η* = *x* − *ξ and using the Fox H representation of the ML function allows us to rewrite the convolution integral as*

$$u(\mathbf{x}) = f(\mathbf{x}) + \lambda \int\_0^\mathbf{x} \eta^{a-1} H\_{1,2}^{1,1} \begin{bmatrix} -\lambda \eta^a \end{bmatrix} \begin{pmatrix} 0, 1 \\ 0, 1 \end{pmatrix}\_{\begin{pmatrix} 1-\beta, a \end{pmatrix}} \begin{bmatrix} \lambda \eta^a \end{bmatrix} f(\mathbf{x} - \eta) d\eta,\tag{50}$$

*which can be rewritten by using properties of the Fox H function [4] as*

$$u(\mathbf{x}) = f(\mathbf{x}) + \lambda \int\_0^\mathbf{x} H\_{1,2}^{1,1} \left[ -\lambda \eta^a \Big| \begin{array}{c} (a-1,1) \\ (a-1,1) \end{array} \begin{array}{c} (a-1,1) \\ (1-\beta+a(a-1),a) \end{array} \right]\_{a,\beta=a} f(\mathbf{x}-\eta)d\eta. \tag{51}$$

*This representation of the solution is suitable for Sinc convolution computations because we already know how to get the discrete Sinc representation of the Fox H function. This data at hand, we only need the Stenger Laplace representation of f to perform the convolution. For details of the numerical approach and implementation see [7,13].*

*Applying the numerical Sinc methods, we can generate solution approximations for different fractional orders α* = *β (see Figure 9). The results in Figure 9 demonstrate that the solution starts at x* = 0 *like the function f (dashed line). A single peak is observed with varying maximum value and location when α* = *β is changed. The location of the maximum shifts to the right if the values for α* = *β become smaller.*

**Figure 9.** Solutions of Equation (47) represented by (51) for varying *α* = *β* (see inset) using an inhomogeneity *<sup>f</sup>*(*x*) = *<sup>x</sup>*3/2 exp(−*x*) (dashed line). The parameter value for *<sup>λ</sup>* <sup>=</sup> 1. Number of Sinc points *N* = 48.

*In astrophysical applications, the Maxwell-Boltzmann distribution is used as a standard model for idealized gases. Abel's Equation (47) can be solved by using a Maxwell-Boltzmann distribution as inhomogeneity. Convolution with the Fox H kernel delivers a shifted distribution as shown in Figure 10 (top graph). The double logarithmic graph (bottom) of Figure 10 for the same computation indicates a similarity observed for the solar neutrino spectrum [38] with a sharp decline on the right side of the function and a linear relation on the left end (in scaled figures).*

**Figure 10.** *Cont*.

**Figure 10.** Solutions of Equation (47) represented by (51) for varying *α* = *β* (see inset) using an inhomogeneity of Maxwell-Boltzmann type *f*(*x*) = *x*<sup>2</sup> exp <sup>−</sup><sup>3</sup> *<sup>x</sup>*<sup>2</sup> 2 (dashed line). The parameter value for *λ* = 1. Number of Sinc points *N* = 32.

The examples demonstrate that our approach to represent numerically Fox *H* functions is highly efficient.

## *3.2.* ℵ *Functions*

This subsection is discussing numerical representations of different ℵ functions as "simple" functions and advanced applications, as introduced in Equation (8). One characteristic of the ℵ function is the weighted sum of Γ functions in the denominator in Mellin space. We will present different simple and advanced rational expressions of Γ functions to represent functions for which their analytic expressions are unknown. In general, we do not know the specific names of the resulting functions, but we are able to generate analytic representations based on Sinc approximations. The Sinc representation in turn allows generalizations of the Fox *H* and Saxena *I* functions. In turn, the Sinc representation allows a numeric computation of function values. Since we do not know what kind of analytic function is represented by the fractions in Mellin space, we shall use a reference representation generated by numerous Sinc points (*N* = 128) to estimate a relative error. This approach to estimate the error of the approximation can be used because we know from (21) and (30) that the Sinc approximation converges exponentially. Due to this fact, the error given in the following examples is always a relative error with respect to this large Sinc point approximation.

Recently some applications of ℵ functions to engineering and biomedical applications were discussed in the literature [39–41] while applications are discussed in [42,43]. However, the authors stop at the analytical representation of their results and do not generate a numerical verification. We will extend their approaches by going one step further to show that a numerical representation of ℵ functions is possible. Thus, the discussed models can be numerically verified. Due to lack of space, we will restrict ourselves to the function representation only.

#### **Example 5.** *Characteristic simple ratio*

*This example uses a rational expression of* Γ *functions in such a way that it cannot be reduced to a Fox H or a Saxena I function. The ratio keeps the characteristic of the* ℵ *function, with a sum of weighted* Γ *functions in the denominator. We examine the following ratio numerically*

$$\mathcal{A}\_{p\_k, q\_k, \tau\_k; r}^{\text{m,n}}(s) = \frac{\Gamma(2s+1)}{2\Gamma(s)+1},\tag{52}$$

*by evaluating the following Mellin-Barnes integral on a line parallel to the imaginary axis*

$$\begin{array}{rcl} I &=& \frac{1}{2\pi i} \int\_{\mathcal{C}} \mathcal{A}\_{p\_k \mathfrak{f}\_k, \mathfrak{r}\_k \mathfrak{r}}^{\mathfrak{m}, \mathfrak{n}}(\mathbf{s}) \ z^{-s} ds = \frac{1}{2\pi i} \int\_{\gamma - i\infty}^{\gamma + i\infty} \mathcal{A}\_{p\_k \mathfrak{f}\_k, \mathfrak{r}\_k \mathfrak{r}}^{\mathfrak{m}, \mathfrak{n}}(\mathbf{s}) \ z^{-s} ds \\\\ &=& \frac{1}{2\pi} \int\_{-\infty}^{\infty} \mathcal{A}\_{p\_k \mathfrak{f}\_k, \mathfrak{r}\_k \mathfrak{r}}^{\mathfrak{m}, \mathfrak{n}}(\gamma - i\sigma) \ z^{-\gamma + i\sigma} d\sigma. \end{array} \tag{53}$$

*A reference approximation using N* = 128 *Sinc points is plotted in the same graph as solid line. Obviously there is no visual difference between the low and high approximation. The used Sinc points m* = 49 *for approximation are shown as dots in the graph. The right panel of Figure 11 shows the local relative error of the approximation. The reference in this graph is the Sinc approximation with the large number of Sinc points. The right panel reflects a mean error of approximately* 10−<sup>8</sup> *which is acceptable for applications.*

**Figure 11.** A true ℵ function using a Mellin representation Γ(1 + 2*s*)/(1 + 2Γ(*s*)) approximated with *N* = 24 Sinc points (left panel). The right panel shows the local relative error where *u*(*x*) is given by a Sinc approximation with *N* = 128.

*First, we verify the accurate representation of the approximation shown in Figure 11. The left panel shows the function approximation using N* = 24*. In Figure 12 we examine the error based on the L*<sup>2</sup> *norm of the relative local error. The error <sup>N</sup> shows an exponential decay if the number of approximation points N is increased. The upper bound of this decay follows the relation (30) which is indicated by the solid line.*

x

**Figure 12.** Error decay  *<sup>N</sup>* <sup>=</sup> *<sup>u</sup>* <sup>−</sup> *<sup>u</sup>*ex <sup>∼</sup> *<sup>K</sup>*2*N*1/2 exp −*k*2*N*1/2 . The least square fit delivered the values (*K*2, *k*2)=(0.007, 2.391). The • and the indicate the error estimation based on the relative *L*<sup>2</sup> norm error and the *L*<sup>2</sup> norm, respectively. In the *L*<sup>2</sup> norms, the exact function is replaced by the Sinc approximation using *N* = 128 Sinc points. It is obvious that the two error estimations deliver nearly the same numerical values.

*We were also curious to see how the* ℵ *function changes when some Mellin representation elements are changed. As a result, we included four factors at various points in (53). Because the* *current form in (53) is so straightforward, we modified it to a single parameter representation as follows:*

$$\mathcal{A}\_{\mathfrak{a}}(\mathbf{s}) = \frac{\Gamma(2s+1)}{\mathfrak{a}\Gamma(s)+1}, \mathcal{A}\_{\mathfrak{f}}(\mathbf{s}) = \frac{\Gamma(\mathfrak{f}s+1)}{\Gamma(s)/2+1},\tag{54}$$

*and*

$$\mathcal{A}\_{\delta}(s) = \frac{\Gamma(s+1)}{\Gamma(s)/2 + \delta'}, \mathcal{A}\_{\eta}(s) = \frac{\Gamma(s+1)}{\Gamma(\eta s)/2 + 1}. \tag{55}$$

*The Mellin-Barnes integral yields the equivalent function*

$$I\_{\parallel} = \frac{1}{2\pi i} \int\_{\mathcal{C}} \mathcal{A}^{m,n}\_{p\_k, q\_k, \eta\_k; r}(s) \, z^{-s} ds = \frac{1}{2\pi i} \int\_{\gamma - i\infty}^{\gamma + i\infty} \mathcal{A}^{m,n}\_{p\_k, q\_k, \eta\_k; r}(s) \, z^{-s} ds$$

$$= \frac{1}{2\pi} \int\_{-\infty}^{\infty} \mathcal{A}^{m,n}\_{p\_k, q\_k, \eta\_k; r}(\gamma - i\sigma) \, z^{-\gamma + i\sigma} d\sigma,\tag{56}$$

*where α, β, δ, and η are all derived from* R+*. Variation of the parameters α, β, δ, and η result in distinct variations of the function according to the definition of the* ℵ *function in (54) and (55). Figure 13 depicts the variation results. We can see that the function behavior varies continuously within a defined region of parameters.*

**Figure 13.** Variation of parameters (indicated in the panels) for A*α*(*s*), A*β*(*s*), A*δ*(*s*), and A*η*(*s*).

**Example 6.** *Two parameter weighted Mittag-Leffler (wML) function*

*It is simple to generalize some types of functions when a tool like the* ℵ *function is available. In the next example, we will demonstrate this. A Fox H function defines the two-parameter Mittag-Leffler (ML) function, as shown in (40). We may introduce the following representation by following the same line of representation as an* ℵ *function.*

$$E\_{a, \beta}^{\delta, \gamma}(s) = \frac{\Gamma(1 - s)\Gamma(s)}{\delta \Gamma(\beta - s a) + \gamma \Gamma(1 - s)} = \mathcal{A}(s). \tag{57}$$

$$\begin{array}{rcl} I &=& \frac{1}{2\pi i} \int\_{\mathcal{C}} \mathcal{A}\_{p\_k q\_k, \mathfrak{r}\_k \mathcal{F}}^{m, n}(s) \ z^{-s} ds = \frac{1}{2\pi i} \int\_{\gamma - i\infty}^{\gamma + i\infty} \mathcal{A}\_{p\_k q\_k, \mathfrak{r}\_k \mathcal{F}}^{m, n}(s) \ z^{-s} ds \\\\ &=& \frac{1}{2\pi} \int\_{-\infty}^{\infty} \mathcal{A}\_{p\_k q\_k, \mathfrak{r}\_k \mathcal{F}}^{m, n}(\gamma - i\sigma) \ z^{-\gamma + i\sigma} d\sigma, \end{array} \tag{58}$$

*The Mellin representation of the extended Mittag-Leffler function is similar to the two parameter Mittag-Leffler function represented as Fox H function with three* Γ *terms. The numerator is exactly the same as the Equation (40). The denominator term is weighted by the coefficients δ*

*and γ. The term containing the weight of δ is the same as (40). The difference exists only in the γ weighted term. If we select γ* = 0 *and δ* = 1*, we reproduce the two parameter Mittag-Leffler function <sup>E</sup>α*,*β*(−*x*)*. To examine the influence of the two weight parameters <sup>δ</sup>, <sup>γ</sup>* <sup>∈</sup> <sup>R</sup>+*, we vary them and represent the results in Figure 14. The Figure indicates that a variation of either of δ or γ in a fractional or integer way will approach the two parameter Mittag-Leffler function from above or below, respectively. In Figure 14 we represent in each row the variation of the parameter δ and γ. The left column of the Figure show changes approaching values larger than one, while the right column shows an approach to smaller values than one. Larger values than one separate the graph of the wML function from the two parameter ML function. Smaller values than one move the graphs of the wML function towards the unweighted ML function, or transpass the ML function in case of δ variation. In addition, we observe that the asymptotic values for larger x values approach a different asymptotic behavior of the ML function. While for γ* = 1 *(top row) the asymptotic slope is nearly the same for δ* = 1 *(bottom row) we observe a fanning out of the asymptotic behavior. In general, we can state that the asymptotic behavior of the wML function assumes a different slope than the ML function, which is typically larger than the slope of the ML function.*

**Figure 14.** Variation of parameters (indicated in the panels) for the weighted Mittag-Leffler function *Eδ*,*<sup>γ</sup> <sup>α</sup>*,*β*(−*x*) = *<sup>E</sup>*˜ *<sup>α</sup>*,*β*(−*x*) using *α* = 1/3 and *β* = 1/2. The top row uses *γ* = 1and the bottom row *δ* = 1. The Mittag-Leffler function *Eα*.*β*(−*x*) is shown as reference (dashed line). The number of Sinc points in all approximations is *N* = 72.

This example demonstrates that well-known functions can be modified straightforwardly if we introduce weights in an ℵ representation in Mellin space. This bears the potential that we can adapt specific functions which do not fit to practical data properly within a well-defined function environment. Such a tool is of utmost practical importance if the data set is of experimental origin.

#### **Example 7.** *Solution of Abel's Integral Equation (47) Generalized*

*The present example deals with Abel's equation's solution (51). We assume that (51) is a function (u(x)) produced by a convolution based on a known function f*(*x*)*. If we view (51) as a convolution integral, we may alter the restrictions α* = *β to α* = *β and/or substitute the ML function with a wML function based on* ℵ *functions. These modifications will produce a convolution integral representation of the type:*

$$u(\mathbf{x}) = f(\mathbf{x}) + \lambda \int\_0^\mathbf{x} \mathbb{N}^{1,1}\_{1,2,\eta,2} \left[ -\lambda \eta^a \Big|\begin{matrix} (a-1,1) \\ (a-1,1) \end{matrix} \begin{matrix} (a-1,1) \\ (1-\beta+a(a-1),a) \end{matrix} \right] f(\mathbf{x}-\eta)d\eta;\tag{59}$$

*i.e., in the convolution, we replace the Fox H function with an* ℵ *function. This is a simple assignment in terms of basic arithmetic. The related integral equation, on the other hand, will convert to an unknown equation. We will employ Sinc convolution techniques to produce the function u*(*x*) *once more, [7,13]. Consider first the case when α* < *β with α* = 1/3 *and β* = 1/2*. As previously stated, the weights γ and δ allow for the adaptation of the function u*(*x*) *to a desired structure, which is more adaptable to a specific instance than the Fox H example. This trait is depicted in Figure 15. The weights allow the curves in the double logarithmic representations to be spread out by approximately four decades in amplitude.*

**Figure 15.** Variation of parameters (indicated in the panels) for the convolution integral (56) using a weighted Mittag-Leffler function *Eδ*,*<sup>γ</sup> <sup>α</sup>*,*β*(−*x*) = *<sup>E</sup>*˜ *<sup>α</sup>*,*β*(−*x*) based on ℵ functions. The parameters are *α* = 1/3 and *β* = 1/2. The top row uses *γ* = 2 and the bottom row *δ* = 2. The number of Sinc points in all approximations is *<sup>N</sup>* <sup>=</sup> 48. The inhomogeneity used is *<sup>f</sup>*(*x*) = *<sup>x</sup>*<sup>2</sup> exp(−3*x*) (dashed line).

*The second situation considered is given by α* > *β with α* = 1/2 *and β* = 1/3*. We also modified the weights δ and γ for this set of parameters, resulting in Figure 16. A comparison of Figures 15 and 16 illustrates that changing the relationship between α and β modifies the original function <sup>f</sup>*(*x*) = *<sup>x</sup>*<sup>2</sup> exp(−3*x*)*. If <sup>α</sup>* <sup>&</sup>lt; *<sup>β</sup> pronounced peaks and minima appear, while just a prominent maximum and shallow minima appear in the opposite case. The position of maxima is marked by dots in both Figures.*

**Figure 16.** *Cont*.

**Figure 16.** Variation of parameters (indicated in the panels) for the convolution integral (56) using a weighted Mittag-Leffler function *Eδ*,*<sup>γ</sup> <sup>α</sup>*,*β*(−*x*) = *<sup>E</sup>*˜ *<sup>α</sup>*,*β*(−*x*) based on ℵ functions. The parameters are *α* = 1/2 and *β* = 1/3. The top row uses *γ* = 2 and the bottom row *δ* = 2. The number of Sinc points in all approximations is *<sup>N</sup>* <sup>=</sup> 48. The inhomogeneity used is *<sup>f</sup>*(*x*) = *<sup>x</sup>*<sup>2</sup> exp(−3*x*) (dashed line).

The examples revealed that we can get numerical values for unknown functions with a few discretization points. The parameter variation allows for a continual modification in the functions, resulting in the desired representation.

#### *3.3. Saxena I Functions*

Rather of examining the Saxena *I* function, we compare the three special functions in this section. As a result, we utilize the case of stable distributions, which Schneider initially looked at in conjunction with Fox *H* functions [44]. Later, Mainardi and Pagnini used Mellin-Barnes integrals to investigate the Schneider technique [45]. This representation of Fox *H* functions in Mellin space will be extended to ℵ- and Saxena *I* functions. We already know that this change in the Mellin model of "stable distributions" leads to various distributions. However, it is fascinating how close some of these functions are to one another. We can use the restrictions on the weights to go back to the Fox *H* function for the ℵ function. The purpose of this section is to compare these classes of functions graphically by graphing their numerical representations.

#### **Example 8.** *Stable Distributions*

*Schneider established stable one-sided Lévy distributions using a Fourier-Stieltjes transform in conjunction with the Mellin transform [44]. As a result, the density representation in Mellin space is as follows:*

$$f\_{\mathfrak{a},\mathfrak{f}}(\mathfrak{s}) = \mathfrak{e} \frac{\Gamma(\mathfrak{s})\Gamma(\mathfrak{e} - \mathfrak{e}\mathfrak{s})}{\Gamma(\mathfrak{\eta} - \mathfrak{\eta}\mathfrak{s})\Gamma(1 - \mathfrak{\eta} + \mathfrak{\eta}\mathfrak{s})}, \\
\text{with } \mathfrak{e} = \mathfrak{a}^{-1} \\
\text{and } \mathfrak{\eta} = \frac{\mathfrak{a} - \mathfrak{\beta}}{2\mathfrak{a}}, \tag{60}$$

*where* 0 < *α* < 2 *and β* = max(*α*, *α* − 2)*. In Mellin space, however, relation (60) is nothing more than a Fox H representation. As a result, we may write*

$$
\hat{\mathcal{H}}\_{2,2}^{1,1}(s) = \epsilon \mathcal{H}\_{2,2}^{1,1} \left( s \middle| \begin{array}{c} (1 - \epsilon, \epsilon) & (1 - \eta, \eta) \\ (0, 1) & (1 - \eta, \eta) \end{array} \right) = \hat{f}\_{a, \emptyset}(s), \tag{61}
$$

*which will deliver using the Mellin-Barnes integral the probability density fα*,*<sup>β</sup>*

$$f\_{a, \delta}(z) = \frac{1}{2\pi i} \int\_{\mathcal{C}} \mathcal{H}^{1, 1}\_{2, 2}(s) \, z^{-s} ds = \frac{1}{2\pi i} \int\_{\gamma - i\infty}^{\gamma + i\infty} \mathcal{H}^{1, 1}\_{2, 2}(s) \, z^{-s} ds$$

$$= \frac{1}{2\pi} \int\_{-\infty}^{\infty} \mathcal{H}^{1, 1}\_{2, 2}(\gamma - i\sigma) \, z^{-\gamma + i\sigma} d\sigma,\tag{62}$$

*with* 0 < *γ* < 1*. Now we may ask what happens if we replace H with either the* ℵ *or the Saxena I function in (62); i.e., we set*

$$\mathcal{A}\_{p\_k,q\_k,\tau\_k;2}^{\mathfrak{m},\mathfrak{n}}(s) = \mathfrak{e} \mathcal{A}\_{p\_k,q\_k,\tau\_k;2}^{\mathfrak{m},\mathfrak{n}}\left(\mathfrak{s}\begin{array}{c} (1-\mathfrak{e},\mathfrak{e}) & (1-\eta,\eta) \\ (0,1) & (1-\eta,\eta) \end{array}\right) = \mathfrak{d}\_{\mathfrak{a},\mathfrak{f}}(s),\tag{63}$$

*or*

$$\hat{\mathcal{L}}\_{p\_k,\rho\_k,1,2}^{\mathfrak{m},\mathfrak{n}}(\mathfrak{s}) = \mathfrak{c} \mathcal{L}\_{p\_k,\rho\_k,1,2}^{\mathfrak{m},\mathfrak{n}}\left(\mathfrak{s} \begin{array}{c} (1-\mathfrak{c},\mathfrak{c}) & (1-\eta,\eta) \\ (0,1) & (1-\eta,\eta) \end{array}\right) = \hat{i}\_{\mathfrak{a},\mathfrak{f}}(\mathfrak{s}),\tag{64}$$

*as integrand in (62) and get aα*,*β*(*z*) *or iα*,*β*(*z*)*, respectively. The results are quite interesting and have some similarities with the original definition of the one-sided stable distribution function. The results are collected for a view examples in Figure 17. We depict the Cauchy-Lorentz case with* (*α*, *β*)=(1, 0)*, the Lévy case using* (*α*, *β*)=(1/2, −1/2)*, and one of the Whittaker cases using* (*α*, *β*)=(2/3, −2/3) *shown in each row, respectively. The Fox column in Figure 17 shows the one-sided stable distribution. The Aleph- and the Saxena columns represent the graphs based on the respective function representation. As a reference, we graph in each row the analytic representation of the function as a dashed line. For Cauchy-Lorentz we have f*1,0(*z*) = 1 *z*<sup>2</sup> + 1 *, the Lévy distribution is <sup>f</sup>*1/2,−1/2(*z*) = *<sup>x</sup>*−3/2*<sup>e</sup>* <sup>−</sup> <sup>1</sup> 4*x* 2 <sup>√</sup>*<sup>π</sup> , and the Whittaker distribution is <sup>f</sup>*2/3,−2/3(*z*) = <sup>√</sup>3/*<sup>π</sup>* exp <sup>−</sup> <sup>4</sup> 2(27*z*2) *<sup>z</sup>*−1*W*1/2,1/6 <sup>4</sup> 27*z*<sup>2</sup> *, where Wα*,*β*(*z*) *is Whittaker's W function. The "distributions" based on Aleph or Saxena are graphed as solid lines. We also included an inset to each graph that shows the approximation's local absolute error. The Aleph and Saxena-based approximations clearly resemble the original distribution in various ways. If the weights are experiencing some limiting process like τ*<sup>1</sup> → 0 *and τ*<sup>2</sup> → 1*, the original distribution can be restored using the Aleph approximation. We don't have the freedom to adjust the numerical approximation in the Saxena approximation; the numerical approximation is set. The local absolute and relative error <sup>N</sup> is a minor number, as seen in the insets. Furthermore, we may infer from Figure 17 that the asymptotic behavior of the Aleph and Saxena approximations take the values of the Fox H function.*

The example shows that, in addition to the Fox and Aleph approximations, the Mellin-Barnes representation technique also works for the Saxena *I* function.

**Figure 17.** Three one sided stable distributions and their Aleph and Saxena variants. The number of Sinc points in all approximations is *N* = 48.

#### **4. Conclusions**

We showed that a Sinc approximation of the Mellin-Barnes integral produces acceptable numerical results for three classes of special functions. We were able to construct numerical representations of the Aleph- and Saxena *I* functions for the first time. The above example demonstrates how simple it is to deal numerically with new and unknown functions. It was observed that the Sinc approximation works effectively with an exponential convergence rate. The calculation benefits greatly from the very small number of discretization points. Furthermore, the approximation approach can handle singularities at the approximation interval's endpoints. This is only foreseeable by an asymptotic analysis for a certain choice of parameters in the Mellin representation of the functions. In this work, we looked at standard integration pathways parallel to the imaginary axis. However, the Mellin-Barnes contour C is not limited to this particular example. In future studies, we are going to investigate alternative contours that may be useful in numerically representing special functions.

**Author Contributions:** Conceptualization, G.B. and N.S.; methodology, G.B.; software, G.B.; validation, G.B. and N.S.; formal analysis, G.B.; investigation, G.B.; writing—original draft preparation, G.B.; writing—review and editing, G.B. and N.S.; visualization, G.B.; supervision, G.B.; project administration, G.B.; All authors have read and agreed to the published version of the manuscript.

**Funding:** This research received no external funding.

**Acknowledgments:** The authors acknowledge the inspiring discussions and valuable ideas by Frank Stenger.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**

