Sinc Numeric Methods for Fox-H, Aleph (ℵ), and Saxena-I Functions

Abstract

The purpose of this study is to offer a systematic, unified approach to the Mellin-Barnes integrals and associated special functions as Fox H, Aleph ℵ, and Saxena I function, encompassing the fundamental features and important conclusions under natural minimal assumptions on the functions in question. The approach’s pillars are the concept of a Mellin-Barnes integral and the Mellin representation of the given function. A Sinc quadrature is used in conjunction with a Sinc approximation of the function to achieve the numerical approximation of the Mellin-Barnes integral. The method converges exponentially and can handle endpoint singularities. We give numerical representations of the Aleph ℵ and Saxena I functions for the first time.

1. Introduction

In the past 40 years, the field of fractional calculus has undergone extraordinary development. The analytic and numeric approaches in fractional calculus created tremendous progress, especially the analytic side generated diverse directions which increased the improvement tremendously [1]. A large number of methods and approaches were developed, generating a consistent framework for analysis and symbolic computations [1,2]. However, the numeric developments are far behind the analytic achievements especially the numeric representation of special functions like Fox H, Aleph $(\aleph$), and Saxena I functions [3]. Such kinds of functions exist nowadays utilized in the analysis of fractional calculus. The special functions also found its way to applications in physics, engineering, and computer science [4]. It turned out during the years that linear transforms like Laplace-, Fourier-, and Mellin transforms play a vital role to generate special functions like Fox-H, ℵ, and Saxena’s I function [1]. For the generation of such function, it is always essential to use the inverse of linear transforms which analytically exists but finally are difficult to compute numerically. The issue with unknown functions is that they could have previously unknown singularities. Naturally, this affects both the choice of the numerical method and the convergence at these singularities. The convergence of the employed numerical technique itself may also be a concern. This was the case if the approximation was calculated using an inverse Laplace transform, as mentioned in [5,6]. In a recent paper, we demonstrated by using Sinc methods that it becomes pretty efficient when Mittag-Leffler functions, a subset of Fox H functions, are the target in connection with an inverse Laplace transform [7]. Mittag-Leffler functions are frequently used in representing solutions of fractional differential or integral equations [8,9]. However, these functions are only a subset of the analytic functions needed to represent the large assortment of possible solutions to fractional equations.

Generalizations of Fox H functions are ℵ functions which also include the class of Saxena I functions [10]. These exceptional functions, which have been investigated analytically but are difficult to obtain numerically, are still a painstaking foundation for fractional calculus today. We aim to offer a numerical technique that solves most numerical difficulties like convergence and the occurrence of singularities by applying Sinc methods to the computation of these functions. Sinc methods initially introduced by Frank Stenger are a powerful numerical tool that allows representing nearly any calculus operation in an efficient and exponentially converging way [11]. For example, Sinc methods allow for computing of definite or indefinite integrals, convolution integrals, linear integral transforms and their inverse, to solve fractional differential and integral equations, and many other practical computations [12,13]. One essential characteristic of Sinc methods is the use of a small number of computing aids; i.e., small programs, a small number of discretization points, less memory, etc., in connection with a high precision output of numerical results [14]. We shall apply these approaches to the numerical computation of Mellin-Barnes integrals used in the presentation of special functions.

Next, we will introduce the definition of Fox, ℵ, and Saxena I functions. In Section 2, we shall introduce the approximation methods needed for this work. The application of these Sinc methods is demonstrated in Section 3. Section 4 summarizes the results and addresses open problems with the current approach.

1.1. The Fox H Function

The Fox H function was introduced by Charles Fox in connection with dual integral equations in 1965 [15]. As he stated at that time “These H functions contain Bessel functions as special cases and my aim is to show that, with the help of a suitable terminology, it is possible to write down a solution by inspection”. Today we know that Fox H functions are a remarkably broad set of functions and include the elementary as well as special functions. The application of these functions is versatile and permits to derive solutions just by “inspection” as Fox noted. A collection of such applications are comprised in the book by Mathai et al. [4] which extends the classical text by Mathai and Saxena [16]. Both monographs concentrate on the part of getting solutions “by inspection”. However, since then there exists a tremendous need and pressure to boil down the solutions to numbers; i.e., to represent the symbolic Fox H solutions as numerical estimations at least or as accurate numbers. Our aim here is to use the analytic and symbolic ideas of Fox and his successors to evaluate such kinds of solutions numerically. To this end let us introduce some notations for these functions.

A Fox function $H_{p,q}^{m,n}\left(z\right)$ is defined via a Mellin-Barnes type integral using integers m, n, p, q such that $0\le m\le q$, $0\le n\le p$, for $a_i,b_j\in \mathbb{C}$ with $\mathbb{C}$, the set of complex numbers, and for ${\alpha }_i,{\beta }_j\in {\mathbb{R}}^{+}=(0,\infty )$$(i=1,2,\dots ,p;j=1,2,\dots ,q)$ in the form

$$H_{p,q}^{m,n}\left(z\left|\begin{array}{c}{\left(a_i,{\alpha }_i\right)}_{1,p}\\ {\left(b_j,{\beta }_j\right)}_{1,q}\end{array}\right.\right)=\frac{1}{2\pi i}{\int }_{\mathcal{C}}{\mathcal{H}}_{p,q}^{m,n}\left(s\right)z^{-s}ds$$

(1)

with

$${\mathcal{H}}_{p,q}^{m,n}\left(\left.\begin{array}{c}{\left(a_i,{\alpha }_i\right)}_{1,p}\\ {\left(b_j,{\beta }_j\right)}_{1,q}\end{array}\right|s\right)=\frac{{\prod }_{j=1}^m\Gamma \left(b_j+{\beta }_{j}s\right){\prod }_{j=1}^n\Gamma \left(1-a_j-{\alpha }_{j}s\right)}{{\prod }_{j=m+1}^q\Gamma \left(1-b_j-{\beta }_{j}s\right){\prod }_{j=n+1}^p\Gamma \left(a_j+{\alpha }_{j}s\right)}.$$

(2)

Here

$$z^{-s}=exp[-s\{log\left|z\right|+iargz\left\}\right],z\ne 0,i=\sqrt{-1},$$

(3)

where $log\left|z\right|$ represents the natural logarithm of $\left|z\right|$ and $argz$ is not necessarily the principal value. An empty product in (2), if it occurs, is taken to be one, and the poles

$$b_{j,l}=\frac{-b_j-l}{{\beta }_j}(j=1,\dots ,m;l=0,1,2,\dots )$$

(4)

of the gamma function $\Gamma \left(b_j+{\beta }_{j}s\right)$ and the poles

$$a_{i,k}=\frac{1-a_i+k}{{\alpha }_i}(i=1,\dots ,n;k=0,1,2,\dots )$$

(5)

of the gamma function $\Gamma \left(1-a_i-{\alpha }_{i}s\right)$ do not coincide:

$${\alpha }_i\left(b_j+l\right)\ne {\beta }_j\left(a_i-k-1\right)(i=1,\dots ,n;j=1,\dots ,m;l,k=0,1,2,\dots ).$$

(6)

The contour $\mathcal{C}$ in (1) is the infinite contour which separates all the poles $b_{j,l}$ in (4) to the left and all the poles $a_{i,k}$ in (5) to the right of $\mathcal{C}$. In fact there are many ways to define $\mathcal{C}$ in the complex plane. However, we will concentrate on the cases where $\mathcal{C}$ is parallel to the imaginary axis in the complex plane. For this reason let $s=\gamma +i\sigma$, where $\gamma$ and $\sigma$ are real; then the contour $\mathcal{C}$ along which the integral of (1) is taken is the straight line whose equation is $\gamma ={\gamma }_0$, where ${\gamma }_0$ is a constant. This line is parallel to the imaginary axis in the complex s plane and separates the poles.

For numeric integration we take $\mathcal{C}$ as a contour starting at the point $\gamma -i\infty$ and terminating at the point $\gamma +i\infty$, where $\gamma \in \mathbb{R}=(-\infty ,\infty )$. To simplify the integration we use the substitutions $s=\gamma +i\sigma$ and $ds=id\sigma$ which delivers

$$\begin{array}{ccc}{H}_{p,q}^{m,n}\left(z\left|\begin{array}{c}{\left(a_i,{\alpha }_i\right)}_{1,p}\\ {\left(b_j,{\beta }_j\right)}_{1,q}\end{array}\right.\right)\hfill & =\hfill & \frac{1}{2\pi i}{\int }_{\mathcal{C}}{\mathcal{H}}_{p,q}^{m,n}\left(s\right)z^{-s}ds=\frac{1}{2\pi i}{\int }_{\gamma -i\infty }^{\gamma +i\infty }{\mathcal{H}}_{p,q}^{m,n}\left(s\right)z^{-s}ds\hfill \\ \hfill & =\hfill & \frac{1}{2\pi }{\int }_{-\infty }^{\infty }{\mathcal{H}}_{p,q}^{m,n}(\gamma +i\sigma )z^{-\gamma -i\sigma }d\sigma \hfill \end{array}$$

(7)

allowing a direct integration to represent the Fox H function at points $z\in \mathbb{C}$. Note that the rightmost integral in (7) is a highly oscillating integral. Such kinds of integrals need special care if treated by standard quadrature methods. We will deal with this problem by using the Sinc quadrature discussed in Section 2. The absolute convergence of the integral can be guaranteed under certain conditions on the parameters of the Fox H function; for details see [3].

1.2. The ℵ Function

The ℵ function was introduced by the authors during the examination of fractional differential equations particularly the drift less Fokker-Planck equation [17,18]. The function is a generalization of Fox H function and allows to handle different initial conditions. Today, the ℵ function is well established and in use in several applications [19,20,21,22,23].

An ℵ function ${\aleph }_{p,q}^{m,n}\left(z\right)$ is defined via a Mellin-Barnes type integral using integers m, n, $p_k$, $q_k$ such that $0\le m\le q_k$, $0\le n\le p_k$, for $a_i,b_j,a_{i,k},b_{i,k}\in \mathbb{C}$ with $\mathbb{C}$, the set of complex numbers, and for ${\alpha }_i,{\beta }_j,{\alpha }_{i,k},{\beta }_{i,k}\in {\mathbb{R}}^{+}=(0,\infty )$$\left(i=1,2,\dots ,p_k;j=1,2,\dots ,q_k\right)$, ${\tau }_kϵ\mathbb{R}$ for k = 1,…,r. The integration path $\mathcal{C}$ extends from $\gamma -\mathrm{i}\infty$ to $\gamma +\mathrm{i}\infty$, and is such that the poles of the gamma functions in the numerator $\Gamma \left(1-a_j-{\alpha }_{j}s\right)$, $j=1,\dots ,n$ do not coincide with the poles of the gamma functions $\Gamma \left(b_j+{\beta }_{j}s\right)$, $j=1,\dots ,m$. The parameters $p_k$ and $q_k$ are non-negative integers satisfying $0\le n\le p_k$, $0\le m\le q_k$. All the poles of the integrand (8) are often assumed to be simple, and the empty product is interpreted as unity. The ℵ function is defined as follows

$${\aleph }_{p_k,q_k,{\tau }_k;r}^{m,n}\left(z\left|\begin{array}{c}{\left(a_i,{\alpha }_i\right)}_{1,p}\\ {\left(b_j,{\beta }_j\right)}_{1,q}\end{array}\right.\right)=\frac{1}{2\pi i}{\int }_{\mathcal{C}}{\mathcal{A}}_{p_k,q_k,{\tau }_k;r}^{m,n}\left(s\right)z^{-s}ds$$

(8)

with the Mellin representation of the kernel ${\mathcal{A}}_{p_k,q_k,{\tau }_k;r}^{m,n}\left(s\right)$

$${\mathcal{A}}_{p_k,q_k,{\tau }_k;r}^{m,n}\left(\left.\begin{array}{c}{\left(a_i,{\alpha }_i\right)}_{1,p}\\ {\left(b_j,{\beta }_j\right)}_{1,q}\end{array}\right|s\right)=\frac{{\prod }_{j=1}^m\Gamma \left(b_j+{\beta }_{j}s\right){\prod }_{j=1}^n\Gamma \left(1-a_j-{\alpha }_{j}s\right)}{{\sum }_{k=1}^r{\tau }_k{\prod }_{j=m+1}^{q_k}\Gamma \left(1-b_{j,k}-{\beta }_{j,k}s\right){\prod }_{j=n+1}^{p_k}\Gamma \left(a_{j,k}+{\alpha }_{j,k}s\right)}.$$

(9)

Here

$$z^{-s}=exp[-s\{log\left|z\right|+iargz\left\}\right],z\ne 0,i=\sqrt{-1}.$$

(10)

Note that the Fox H function follows from the ℵ function in case when $r=1$ and ${\tau }_k=1$. If ${\tau }_k=1$ for $k=1,\dots ,r$ the Aleph function reduces to a Saxena I function [10]. According to the definition of $\mathcal{C}$ we are numerically dealing with a Bromwich integral in the form

$$\begin{array}{ccc}\hfill {\aleph }_{p_k,q_k,{\tau }_k;r}^{m,n}\left(z\left|\begin{array}{c}{\left(a_i,{\alpha }_i\right)}_{1,p}\\ {\left(b_j,{\beta }_j\right)}_{1,q}\end{array}\right.\right)& =& \frac{1}{2\pi i}{\int }_{\mathcal{C}}{\mathcal{A}}_{p_k,q_k,{\tau }_k;r}^{m,n}\left(s\right)z^{-s}ds\hfill \\ & =& \frac{1}{2\pi i}{\int }_{\gamma -i\infty }^{\gamma +i\infty }{\mathcal{A}}_{p_k,q_k,{\tau }_k;r}^{m,n}\left(s\right)z^{-s}ds\hfill \\ & =& \frac{1}{2\pi }{\int }_{-\infty }^{\infty }{\mathcal{A}}_{p_k,q_k,{\tau }_k;r}^{m,n}(\gamma +i\sigma )z^{-\gamma -i\sigma }d\sigma .\hfill \end{array}$$

(11)

The Aleph function in the present form is the result of solving integral and differential equations using linear transform techniques and is now considered the most generalized special function of a function representation [22].

1.3. The Saxena $\mathcal{I}$ Function

A Saxena I function $I_{p_k,q_k,1;r}^{m,n}\left(z\right)$ is similarly defined as an ℵ function via a Mellin-Barnes type integral using integers m, n, $p_k$, $q_k$ such that $0\le m\le q_k$, $0\le n\le p_k$, for $a_i,b_j,a_{i,k},b_{j,k}\in \mathbb{C}$ with $\mathbb{C}$, the set of complex numbers, and for ${\alpha }_i,{\beta }_j,{\alpha }_{i,k},{\beta }_{j,k}\in {\mathbb{R}}^{+}=(0,\infty )$$\left(i=1,2,\dots ,p_k;j=1,2,\dots ,q_k\right)$, and $k=1,\dots ,r$ [10], in the form

$$I_{p_k,q_k,1;r}^{m,n}\left(z\left|\begin{array}{c}{\left(a_i,{\alpha }_i\right)}_{1,p}\\ {\left(b_j,{\beta }_j\right)}_{1,q}\end{array}\right.\right)=\frac{1}{2\pi i}{\int }_{\mathcal{C}}{\mathcal{I}}_{p_k,q_k,1;r}^{m,n}\left(s\right)z^{-s}ds,$$

(12)

with the Mellin representation of the kernel ${\mathcal{I}}_{p_k,q_k,1;r}^{m,n}\left(s\right)$

$${\mathcal{I}}_{p_k,q_k,1;r}^{m,n}\left(\left.\begin{array}{c}{\left(a_i,{\alpha }_i\right)}_{1,p}\\ {\left(b_j,{\beta }_j\right)}_{1,q}\end{array}\right|s\right)=\frac{{\prod }_{j=1}^m\Gamma \left(b_j+{\beta }_{j}s\right){\prod }_{j=1}^n\Gamma \left(1-a_j-{\alpha }_{j}s\right)}{{\sum }_{k=1}^r{\prod }_{j=m+1}^{q_k}\Gamma \left(1-b_{j,k}-{\beta }_{j,k}s\right){\prod }_{j=n+1}^{p_k}\Gamma \left(a_{j,k}+{\alpha }_{j,k}s\right)}.$$

(13)

Here

$$z^{-s}=exp[-s\{log\left|z\right|+iargz\left\}\right],z\ne 0,i=\sqrt{-1}.$$

(14)

The conditions for the contour $\mathcal{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 $\mathcal{D}$ be a simply connected domain in the complex plane and $z\in \mathbb{C}$ having a boundary $\partial \mathcal{D}$. Let a and b denote two distinct points of $\partial \mathcal{D}$ and ϕ denote a conformal map of $\mathcal{D}$ onto ${\mathcal{D}}_d$, where ${\mathcal{D}}_d=\{z\in \mathbb{C}:|\mathcal{I}\mathit{m}\left(z\right)|<d\}$, such that $\varphi \left(a\right)=-\infty$ and $\varphi \left(b\right)=\infty$. Let $\psi ={\varphi }^{-1}$ denote the inverse conformal map, and let Γ be an arc defined by $\Gamma =\{z\in \mathbb{C}:z=\psi (x),x\in \mathbb{R}\}$. Given ϕ, ψ, and a positive number h, let us set $z_k=\psi \left(kh\right)$, $k\in \mathbb{Z}$ to be the Sinc points, let us also define $\rho \left(z\right)=e^{\varphi \left(z\right)}$.

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\in (0,\pi )$, and let the domains $\mathcal{D}$ and ${\mathcal{D}}_d$ be given as in Definition 1. If $d^{\prime }$ is a number such that $d^{\prime }>d$, and if the function ϕ provides a conformal map of ${\mathcal{D}}^{\prime }$ onto ${\mathcal{D}}_{d^{\prime }}$, then $\mathcal{D}\subset {\mathcal{D}}^{\prime }$. Let α and β denote positive numbers, and let ${\mathit{L}}_{\alpha ,\beta }\left(\mathcal{D}\right)$ denote the family of functions $u\in \mathrm{Hol}\left(\mathcal{D}\right)$, for which there exists a positive constant $c_1$ such that, for all $z\in \mathcal{D}$. Let ${\mathit{L}}_{\alpha ,\beta }\left(\mathcal{D}\right)$ be the set of all analytic functions, for which there exists a constant $c_1$, such that

$$\left|u\left(z\right)\right|\le c_1\frac{{\left|\rho \left(z\right)\right|}^{\alpha }}{{(1+|\rho \left(z\right)\left|\right)}^{\alpha +\beta }}.$$

(15)

Now let the positive numbers α and β belong to $(0,1]$, and let ${\mathit{M}}_{\alpha ,\beta }\left(\mathcal{D}\right)$ denote the family of all functions $g\in \mathrm{Hol}\left(\mathcal{D}\right)$, such that $g\left(a\right)$ and $g\left(b\right)$ are finite numbers, where $g\left(a\right)={lim}_{z\to a}g\left(z\right)$ and $g\left(b\right)={lim}_{z\to b}g\left(z\right)$, and such that $u\in {\mathit{L}}_{\alpha ,\beta }\left(\mathcal{D}\right)$ where

$$u\left(z\right)=g\left(z\right)-\frac{g\left(a\right)+\rho \left(z\right)g\left(b\right)}{1+\rho \left(z\right)}.$$

(16)

These definitions directly allow the formulation of an algorithm for a Sinc approximation. Let $\mathbb{Z}$ denote the set of all integers. Select positive integers N and $M=[\beta N/\alpha ]$ so that $m=M+N+1$. The step length is determined by $h={(\pi d/\left(\beta N\right))}^{1/2}$ where $\alpha$, $\beta$, and d are real parameters. In addition assume there is a conformal map $\varphi$ and its inverse $\psi$ such that we can define Sinc points $z_j=\psi \left(jh\right)$, $j\in \mathbb{Z}$ [25]. The following relations define the basis of a Sinc approximation:

$$\mathrm{Sin}\mathrm{c}\left(z\right)=\frac{sin\left(\pi z\right)}{\pi z}.$$

(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 \left(z\right)=sin\left(\pi \right(z/h-k\left)\right)/\left(\pi \right(z/h-k\left)\right)=\mathrm{Sin}\mathrm{c}(z/h-k).$$

(18)

The approximation of a function $f\left(z\right)$ results to the representation

$$C_{h,M,N}\left[f\right]\left(z\right)=\sum _{k=-M}^{N}f\left(z_k\right)\mathrm{Sin}\mathrm{c}\left(\frac{\varphi \left(z\right)}{h}-k\right)=\sum _{k=-M}^{N}f\left(z_k\right)S(k,h)\circ \varphi \left(z\right),$$

(19)

using the set of orthogonal functions

$$B_S={\{S(k,h)\circ \varphi \left(z\right)\}}_{k=-M}^N,$$

(20)

where $\varphi \left(z\right)$ is the conformal map. This type of approximation allows to represent a function $f\left(z\right)$ on an arc $\Gamma$ 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\in {\mathit{L}}_{\alpha ,\beta }\left(\mathcal{D}\right)$ for $\alpha >0$ and $\beta >0$, take M = [β N/α], where [x] denotes the greatest integer in x, and then set $m=M+N+1$. If $u\in {\mathit{M}}_{\alpha ,\beta }\left(\mathcal{D}\right)$, and if $h={(\pi d/\left(\beta N\right))}^{1/2}$ then there exist positive constants $K_1$ and $k_1$ independent of N, such that

$${ϵ}_N=∥f\left(z\right)-C_{h,M,N}\left[f\right]\left(z\right)∥=K_1{N}^{1/2}exp\left(-k_1{N}^{1/2}\right).$$

(21)

with $w_k$ the base function (see (22)).

The proof of this theorem is given in [13]. Note the choice $h={(\pi d/\left(\beta N\right))}^{1/2}$ is close to optimal for an approximation in the space ${\mathit{M}}_{\alpha ,\beta }\left(\mathcal{D}\right)$ 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 $z_k=\psi \left(kh\right)={\varphi }^{-1}\left(kh\right)$ 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 $w_j\left(z\right)$ by

$$w_j=\left\{\begin{array}{cc}\frac{1}{1+\rho \left(z\right)}-\sum _{k=-M+1}^N\frac{1}{1+e^{kh}}S(k,h)\circ \varphi \left(z\right)\hfill & j=-M\hfill \\ S(k,h)\circ \varphi \left(z\right)\hfill & j=-M+1,\dots ,N-1\hfill \\ \frac{\rho \left(z\right)}{1+\rho \left(z\right)}-\sum _{k=-M}^{N-1}\frac{e^{kh}}{1+e^{kh}}S(k,h)\circ \varphi \left(z\right)\hfill & j=N\hfill \end{array}\right.,$$

(22)

where $\rho \left(z\right)=exp\left(\varphi \right(z\left)\right)$. 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 ${\mathit{V}}_m\left(S\right)$ of basis functions

$${\mathit{V}}_m\left(S\right)=\left(w_{-M},\dots ,w_N\right)$$

(23)

with $w_j$ defined as in (22). For a given vector ${\mathit{V}}_m\left(u\right)=\left(u_{-M},\dots ,u_N\right){}^T$ we now introduce the dot product as an approximation of the function $u\left(z\right)$ by

$$u\left(z\right)\approx V_m\left(S\right).{\mathit{V}}_m\left(u\right)=\sum _{k=-M}^N{u}_k{w}_k.$$

(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 $\mathbb{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\left(f\right)={\int }_a^{b}f\left(t\right)dt,$$

(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\left(f\right)$,

$$\begin{array}{ccc}\hfill J\left(f\right)& =& {\int }_a^{b}f\left(t\right)dt\approx \hfill \end{array}$$

$$\begin{array}{ccc}\hfill J_{h,M,N}\left[f\right]\left(x\right)& =& {\int }_a^b\sum _{k=-M}^{N}f\left(t_k\right)S(k,h)\circ \varphi \left(t\right)dt\hfill \\ \hfill \end{array}$$

(26)

$$\begin{array}{ccc}& =& \sum _{k=-M}^{N}f\left(t_k\right){\int }_a^{b}S(k,h)\circ \varphi \left(t\right)dt.\hfill \end{array}$$

(27)

Scaling the variable $\xi =t/h$ and collocating the expression with respect to $\xi$, we end up with the representation

$$J\left(f\right)\approx h\sum _{k=-M}^{N}f\left(t_k\right)\frac{1}{{\varphi }^{\prime }\left(t_k\right)},\mathrm{with}{t}_k=\psi \left(kh\right).$$

(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 $t_k=\psi \left(kh\right)$ with $k=-M,\dots ,N$ and $h=\pi /\sqrt{N}$. Note since the conformal map $\varphi$ depends on the structure of the arc $\Gamma$; i.e., finite, semi-infinite or infinite, the approximation is defined for domains $[a,b]$, $(0,\infty )$, or $(-\infty ,\infty )$, 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 $\mathbb{R}$, let h denote a fixed positive number, and let the Sinc points be defined on $(a,b)$ by $z_k={\varphi }^{-1}\left(kh\right)$, $k\in \mathbb{Z}$, where ${\varphi }^{-1}=\psi$ 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 ${\mathit{V}}_m\left(f\right)=\left(f\left(z_{-M}\right),\dots ,f\left(z_N\right)\right)$, and a vector ${\mathit{V}}_m(1/{\varphi }^{\prime })=\left(1/{\varphi }^{\prime }\left(z_{-M}\right),\dots ,1/{\varphi }^{\prime }\left(z_N\right)\right)$, then the definite integral is approximated by

$$J\left(f\right)\approx h{\mathit{V}}_m\left(f\right).{\mathit{V}}_m(1/{\varphi }^{\prime })=J_{h,M,N}\left[f\right]\left(z\right),$$

(29)

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

$${ϵ}_N=∥J\left(f\right)-J_{h,M,N}\left[f\right]\left(z\right)∥\sim K_2{N}^{1/2}exp\left(-k_2{N}^{1/2}\right),$$

(30)

where $K_2$ and $k_2$ 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 $\Gamma$ 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 $\gamma$ parallel to the imaginary axis in $s\in \mathbb{C}$. In formulas this means

$$\begin{array}{ccc}{H}_{p,q}^{m,n}\left(z_k|\begin{array}{c}{\left(a_i,{\alpha }_i\right)}_{1,p}\\ {\left(b_j,{\beta }_j\right)}_{1,q}\end{array}\right)\hfill & =\hfill & \frac{1}{2\pi }{\int }_{-\infty }^{\infty }{\mathcal{H}}_{p,q}^{m,n}(\gamma +i\sigma )z_k{}^{-\gamma -i\sigma }d\sigma \hfill \\ \hfill & \approx \hfill & \frac{h}{2\pi }{\mathit{V}}_m\left({\mathcal{H}}_{p,q}^{m,n}\left(\gamma +i{\sigma }_l\right)z_k{}^{-\gamma -i{\sigma }_l}\right).{\mathit{V}}_m\left(1/{\varphi }^{\prime }\left({\sigma }_l\right)\right),\hfill \end{array}$$

(31)

where ${\sigma }_l=\psi \left(lh\right),l=-M,\dots ,N,\mathrm{and}{z}_k=\psi \left(kh\right),k=-M,\dots ,N$, generating a vector ${\mathit{V}}_m\left(H_{p,q}^{m,n}\right)$ which approximates the Fox H function pointwise at $z_k=\psi \left(kh\right)$ on an arc $\Gamma$. $\gamma$ is selected according to the pole structure of the $\Gamma$ functions of the numerator in the Mellin representation. Using the basis functions $w_k$ (22) on the arc $\Gamma$, we shall approximate the Fox H function by

$$H_{p,q}^{m,n}\left(z\right)\approx V_m\left(S\right).{\mathit{V}}_m\left(H_{p,q}^{m,n}\right)=\sum _{k=-M}^N{H}_{p,q}^{m,n}\left(z_k\right)w_k\left(z\right).$$

(32)

This approach guarantees that the error ${ϵ}_N$ 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 $\Gamma$.

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;r}^{m,n}\left(z\right)\approx V_m\left(S\right).{\mathit{V}}_m\left({\aleph }_{p_k,q_k,{\tau }_k;r}^{m,n}\right)=\sum _{k=-M}^N{\aleph }_{p_k,q_k,{\tau }_k;r}^{m,n}\left(z_k\right)w_k\left(z\right),$$

(33)

and the Saxena I function

$$I_{p_k,q_k,1;r}^{m,n}\left(z\right)\approx V_m\left(S\right).{\mathit{V}}_m\left(I_{p_k,q_k,1;r}^{m,n}\right)=\sum _{k=-M}^N{I}_{p_k,q_k,1;r}^{m,n}\left(z_k\right)w_k\left(z\right).$$

(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^{-x}=\frac{1}{2\pi i}{\int }_{\gamma -i\infty }^{\gamma +i\infty }\Gamma \left(s\right)x^{-s}ds,\gamma >0,$$

(35)

and

$$\Gamma \left(s\right)={\int }_0^{\infty }{e}^{-x}{x}^{s-1}dx,$$

(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 $\gamma =\mathrm{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^{-x}=H_{0,2}^{2,0}\left(\frac{x}{2}\left|\begin{array}{cc}& \\ \left(0,\frac{1}{2}\right)& \left(\frac{1}{2},\frac{1}{2}\right)\end{array}\right.\right),with{\mathcal{H}}_{0,2}^{2,0}\left(s\right)=\frac{\Gamma \left(\frac{s}{2}+\frac{1}{2}\right)\Gamma \left(\frac{s}{2}\right)}{2\sqrt{\pi }}.$$

(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{array}{ccc}\hfill e^{-z}& =& \frac{1}{2\pi i}{\int }_{\mathcal{C}}{\mathcal{H}}_{0,2}^{2,0}\left(s\right)z^{-s}ds=\frac{1}{2\pi i}{\int }_{\gamma -i\infty }^{\gamma +i\infty }{\mathcal{H}}_{0,2}^{2,0}\left(s\right)z^{-s}ds\hfill \\ & =& \frac{1}{2\pi }{\int }_{-\infty }^{\infty }{\mathcal{H}}_{0,2}^{2,0}(\gamma -i\sigma )z^{-\gamma +i\sigma }d\sigma ,\hfill \end{array}$$

(38)

where the path $\mathcal{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 $\sigma \in (-\infty ,\infty )$. The set of values ${\left\{z_k,I_k\right\}}_{k=-N}^N$ is used in a Sinc approximation delivering an analytic representation of the function $f\left(z\right)=e^{-z}$ 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_2$ and $k_2$ 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).

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_{\alpha }\left(z\right)=H_{1,2}^{1,1}\left(-z\left|\begin{array}{cc}(0,1)& \\ (0,1)& (0,\alpha )\end{array}\right.\right),with{\mathcal{H}}_{1,2}^{1,1}\left(s\right)=\frac{\Gamma (1-s)\Gamma \left(s\right)}{\Gamma (1-s\alpha )},$$

(39)

$$E_{\alpha ,\beta }\left(z\right)=H_{1,2}^{1,1}\left(-z\left|\begin{array}{cc}(0,1)& \\ (0,1)& (1-\beta ,\alpha )\end{array}\right.\right),with{\mathcal{H}}_{1,2}^{1,1}\left(s\right)=\frac{\Gamma (1-s)\Gamma \left(s\right)}{\Gamma (\beta -s\alpha )},$$

(40)

and

$$E_{\alpha ,\beta }^{\mu }\left(z\right)=H_{1,2}^{1,1}\left(-z\left|\begin{array}{cc}(1-\mu ,1)& \\ (0,1)& (1-\beta ,\alpha )\end{array}\right.\right),with{\mathcal{H}}_{1,2}^{1,1}\left(s\right)=\frac{\Gamma (\mu -s)\Gamma \left(s\right)}{\Gamma (\beta -s\alpha )}.$$

(41)

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

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

(42)

where the path $\mathcal{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 $\sigma \in (-\infty ,\infty )$. The set of values ${\left\{z_k,I_k\right\}}_{k=-N}^N$ 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.

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_{\rho }^{\nu }\left(z\right)$ was defined by Krätzel [31] as the kernel of an integral transform as follows:

$$Z_{\rho }^{\nu }\left(z\right)={\int }_0^{\infty }{y}^{\nu -1}{e}^{-y^{\rho }-x/y}dy.$$

(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

$$\begin{array}{ccc}\hfill f\left(z\right)& =& c{z}^{\alpha -1}{Z}_{\rho }^{\nu }\left(z\right)=\frac{1}{\Gamma \left(\alpha \right)\Gamma \left(\frac{\alpha +\nu }{\rho }\right)}{H}_{0,2}^{2,0}\left(z\left|\begin{array}{cc}& \\ (\alpha -1,1)& \left(\frac{\alpha +\nu -1}{\rho },\frac{1}{\rho }\right)\end{array}\right.\right),with\hfill \\ & & {\mathcal{H}}_{0,2}^{2,0}\left(s\right)=\Gamma (s+\alpha -1)\Gamma \left(\frac{s}{\rho }+\frac{\alpha +\nu -1}{\rho }\right)\hfill \end{array}$$

(44)

here $x\ge 0$, $\alpha >0$, $\rho >0$, and $\nu >0$. The corresponding line integral representation is computed by the Mellin-Barnes integral via

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

(45)

where the path $\mathcal{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 $\sigma \in (-\infty ,\infty )$. The set of values ${\left\{z_k,I_k\right\}}_{k=-N}^N$ is used in a Sinc approximation. An example using specific parameters $\alpha =2$ and $\rho =\nu =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 $z_k$ 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.

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

For the special choice $\rho =1$, the Krätzel density is reduced to the modified Bessel function $K_{\nu }\left(x\right)$ with the specific representation

$$f\left(z\right)=c{z}^{\alpha -1}{Z}_{\rho }^{\nu }\left(z\right)=\frac{2}{\Gamma \left(\alpha \right)\Gamma (\alpha +\nu )}{x}^{\alpha +\nu /2-1}{K}_{\nu }\left(2\sqrt{x}\right),withx\ge 0,\alpha >0,\nu >0.$$

(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.

Example 4.

Abel Equation

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

$$u\left(x\right)=f\left(x\right)+\frac{\lambda }{\Gamma \left(\alpha \right)}{\int }_0^x{(x-\xi )}^{\alpha -1}u\left(\xi \right)d\xi ,$$

(47)

where λ and $0<\alpha <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

$$u\left(x\right)=f\left(x\right)+\lambda {\mathcal{D}}_{0,x}^{-\alpha }u\left(x\right).$$

(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 $f\left(x\right)=N_0$ and $\lambda =-c^{\alpha }$. 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\left(x\right)$ 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

$$u\left(x\right)=f\left(x\right)+\lambda {\int }_0^x{(x-\xi )}^{\alpha -1}{E}_{\alpha ,\alpha }\left(\lambda {(x-\xi )}^{\alpha }\right)f\left(\xi \right)d\xi .$$

(49)

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

$$u\left(x\right)=f\left(x\right)+\lambda {\int }_0^x{\eta }^{\alpha -1}{H}_{1,2}^{1,1}\left[-\lambda {\eta }^{\alpha }\left|\begin{array}{cc}(0,1)& \\ (0,1)& (1-\beta ,\alpha )\end{array}\right.\right]{}_{\alpha ,\beta =\alpha }\left(\lambda {\eta }^{\alpha }\right)f(x-\eta )d\eta ,$$

(50)

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

$$u\left(x\right)=f\left(x\right)+\lambda {\int }_0^x{H}_{1,2}^{1,1}\left[-\lambda {\eta }^{\alpha }\left|\begin{array}{cc}(\alpha -1,1)& \\ (\alpha -1,1)& (1-\beta +\alpha (\alpha -1),\alpha )\end{array}\right.\right]{}_{\alpha ,\beta =\alpha }f(x-\eta )d\eta .$$

(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 $\alpha =\beta$ (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 $\alpha =\beta$ is changed. The location of the maximum shifts to the right if the values for $\alpha =\beta$ become smaller.

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).

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 $\Gamma$ functions in the denominator in Mellin space. We will present different simple and advanced rational expressions of $\Gamma$ 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,40,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}^{m,n}\left(s\right)=\frac{\Gamma (2s+1)}{2\Gamma \left(s\right)+1},$$

(52)

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

$$\begin{array}{ccc}I\hfill & =\hfill & \frac{1}{2\pi i}{\int }_{\mathcal{C}}{\mathcal{A}}_{p_k,q_k,{\tau }_k;r}^{m,n}\left(s\right)z^{-s}ds=\frac{1}{2\pi i}{\int }_{\gamma -i\infty }^{\gamma +i\infty }{\mathcal{A}}_{p_k,q_k,{\tau }_k;r}^{m,n}\left(s\right)z^{-s}ds\hfill \\ \hfill & =\hfill & \frac{1}{2\pi }{\int }_{-\infty }^{\infty }{\mathcal{A}}_{p_k,q_k,{\tau }_k;r}^{m,n}(\gamma -i\sigma )z^{-\gamma +i\sigma }d\sigma .\hfill \\ \hfill \end{array}$$

(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}^{-8}$ which is acceptable for applications.

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_2$ norm of the relative local error. The error ${ϵ}_N$ 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.

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}}_{\alpha }\left(s\right)=\frac{\Gamma (2s+1)}{\alpha \Gamma \left(s\right)+1},{\mathcal{A}}_{\beta }\left(s\right)=\frac{\Gamma (\beta s+1)}{\Gamma \left(s\right)/2+1},$$

(54)

and

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

(55)

The Mellin-Barnes integral yields the equivalent function

$$\begin{array}{ccc}\hfill I& =& \frac{1}{2\pi i}{\int }_{\mathcal{C}}{\mathcal{A}}_{p_k,q_k,{\tau }_k;r}^{m,n}\left(s\right)z^{-s}ds=\frac{1}{2\pi i}{\int }_{\gamma -i\infty }^{\gamma +i\infty }{\mathcal{A}}_{p_k,q_k,{\tau }_k;r}^{m,n}\left(s\right)z^{-s}ds\hfill \\ & =& \frac{1}{2\pi }{\int }_{-\infty }^{\infty }{\mathcal{A}}_{p_k,q_k,{\tau }_k;r}^{m,n}(\gamma -i\sigma )z^{-\gamma +i\sigma }d\sigma ,\hfill \end{array}$$

(56)

where α, β, δ, and η are all derived from ${\mathbb{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.

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_{\alpha ,\beta }^{\delta ,\gamma }\left(s\right)=\frac{\Gamma (1-s)\Gamma \left(s\right)}{\delta \Gamma (\beta -s\alpha )+\gamma \Gamma (1-s)}=\mathcal{A}\left(s\right).$$

(57)

$$\begin{array}{ccc}I\hfill & =\hfill & \frac{1}{2\pi i}{\int }_{\mathcal{C}}{\mathcal{A}}_{p_k,q_k,{\tau }_k;r}^{m,n}\left(s\right)z^{-s}ds=\frac{1}{2\pi i}{\int }_{\gamma -i\infty }^{\gamma +i\infty }{\mathcal{A}}_{p_k,q_k,{\tau }_k;r}^{m,n}\left(s\right)z^{-s}ds\hfill \\ \hfill & =\hfill & \frac{1}{2\pi }{\int }_{-\infty }^{\infty }{\mathcal{A}}_{p_k,q_k,{\tau }_k;r}^{m,n}(\gamma -i\sigma )z^{-\gamma +i\sigma }d\sigma ,\hfill \\ \hfill \end{array}$$

(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 $\gamma =0$ and $\delta =1$, we reproduce the two parameter Mittag-Leffler function $E_{\alpha ,\beta }(-x)$. To examine the influence of the two weight parameters δ, $\gamma \in {\mathbb{R}}^{+}$, 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 $\gamma =1$ (top row) the asymptotic slope is nearly the same for $\delta =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.

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\left(x\right)$. If we view (51) as a convolution integral, we may alter the restrictions $\alpha =\beta$ to $\alpha \ne \beta$ 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\left(x\right)=f\left(x\right)+\lambda {\int }_0^x{\aleph }_{1,2,{\tau }_k;2}^{1,1}\left[-\lambda {\eta }^{\alpha }\left|\begin{array}{cc}(\alpha -1,1)& \\ (\alpha -1,1)& (1-\beta +\alpha (\alpha -1),\alpha )\end{array}\right.\right]f(x-\eta )d\eta ;$$

(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\left(x\right)$ once more, [7,13]. Consider first the case when $\alpha <\beta$ with $\alpha =1/3$ and $\beta =1/2$. As previously stated, the weights γ and δ allow for the adaptation of the function $u\left(x\right)$ 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.

The second situation considered is given by $\alpha >\beta$ with $\alpha =1/2$ and $\beta =1/3$. We also modified the weights δ and γ for this set of parameters, resulting in Figure 16. A comparison of Figure 15 and Figure 16 illustrates that changing the relationship between α and β modifies the original function $f\left(x\right)=x^{2}exp(-3x)$. If $\alpha <\beta$ 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.

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:

$${\widehat{f}}_{\alpha ,\beta }\left(s\right)=ϵ\frac{\Gamma \left(s\right)\Gamma (ϵ-ϵs)}{\Gamma (\eta -\eta s)\Gamma (1-\eta +\eta s)},withϵ={\alpha }^{-1}\mathrm{and}\eta =\frac{\alpha -\beta }{2\alpha },$$

(60)

where $0<\alpha <2$ and $\beta =max(\alpha ,\alpha -2)$. In Mellin space, however, relation (60) is nothing more than a Fox H representation. As a result, we may write

$${\widehat{\mathcal{H}}}_{2,2}^{1,1}\left(s\right)=ϵ{\mathcal{H}}_{2,2}^{1,1}\left(s\left|\begin{array}{cc}(1-ϵ,ϵ)& (1-\eta ,\eta )\\ (0,1)& (1-\eta ,\eta )\end{array}\right.\right)={\widehat{f}}_{\alpha ,\beta }\left(s\right),$$

(61)

which will deliver using the Mellin-Barnes integral the probability density $f_{\alpha ,\beta }$

$$\begin{array}{ccc}\hfill f_{\alpha ,\beta }\left(z\right)& =& \frac{1}{2\pi i}{\int }_{\mathcal{C}}{\widehat{\mathcal{H}}}_{2,2}^{1,1}\left(s\right)z^{-s}ds=\frac{1}{2\pi i}{\int }_{\gamma -i\infty }^{\gamma +i\infty }{\widehat{\mathcal{H}}}_{2,2}^{1,1}\left(s\right)z^{-s}ds\hfill \\ & =& \frac{1}{2\pi }{\int }_{-\infty }^{\infty }{\widehat{\mathcal{H}}}_{2,2}^{1,1}(\gamma -i\sigma )z^{-\gamma +i\sigma }d\sigma ,\hfill \end{array}$$

(62)

with $0<\gamma <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

$${\widehat{\mathcal{A}}}_{p_k,q_k,{\tau }_k;2}^{m,n}\left(s\right)=ϵ{\mathcal{A}}_{p_k,q_k,{\tau }_k;2}^{m,n}\left(s\left|\begin{array}{cc}(1-ϵ,ϵ)& (1-\eta ,\eta )\\ (0,1)& (1-\eta ,\eta )\end{array}\right.\right)={\widehat{a}}_{\alpha ,\beta }\left(s\right),$$

(63)

or

$${\widehat{\mathcal{I}}}_{p_k,q_k,1;2}^{m,n}\left(s\right)=ϵ{\mathcal{I}}_{p_k,q_k,1;2}^{m,n}\left(s\left|\begin{array}{cc}(1-ϵ,ϵ)& (1-\eta ,\eta )\\ (0,1)& (1-\eta ,\eta )\end{array}\right.\right)={\widehat{i}}_{\alpha ,\beta }\left(s\right),$$

(64)

as integrand in (62) and get $a_{\alpha ,\beta }\left(z\right)$ or $i_{\alpha ,\beta }\left(z\right)$, 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 $(\alpha ,\beta )=(1,0)$, the Lévy case using $(\alpha ,\beta )=(1/2,-1/2)$, and one of the Whittaker cases using $(\alpha ,\beta )=(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}\left(z\right)=1/\left(z^2+1\right)$, the Lévy distribution is $f_{1/2,-1/2}\left(z\right)=\frac{x^{-3/2}{e}^{-\frac{1}{4x}}}{2\sqrt{\pi }}$, and the Whittaker distribution is $f_{2/3,-2/3}\left(z\right)=\sqrt{3/\pi }exp\left(-\frac{4}{2\left(27z^2\right)}\right)z^{-1}{W}_{1/2,1/6}\left(\frac{4}{27z^2}\right)$, where $W_{\alpha ,\beta }\left(z\right)$ 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 ${\tau }_1\to 0$ and ${\tau }_2\to 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 ${ϵ}_N$ 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.

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 $\mathcal{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.