On the analysis of fractional calculus operators with bivariate Mittag Leffler function in the kernel

Abstract

Bivariate Mittag-Leffler (ML) functions are a substantial generalization of the univariate ML functions, which are widely recognized for their significance in fractional calculus. In the present paper, our initial focus is to investigate the fractional calculus properties of the integral and derivative operators with kernels including the Bivariate ML functions. Further, certain fractional Cauchy-type problems including these operators are considered. Also the numerical approximations of the Caputo type derivative operator are investigated. The theoretical results are justified by applications on examples. Furthermore, the theory of applying the same operators with respect to arbitrary monotonic functions is analyzed in this research.

1 Introduction

Fractional calculus is of great importance in applied mathematics and mathematical analysis, with substantial connections to ML functions. The function is defined by the following series when the real part of $\alpha$ is strictly positive, [1]

$$\begin{aligned} E_{\alpha }(t)=\sum _{s=0}^\infty \frac{t^s}{\Gamma (\alpha s+1)}, \end{aligned}$$

where $\Gamma$ is the gamma function. For $\alpha =1$ the ML function reduces to the exponential function and ML functions were used to define the solutions to fractional differential equations. It is well known that fractional differential equations are among the strongest tools of mathematical modeling and are successfully employed to model complex physical and biological phenomena. Such as impulsive neutral Hilfer fractional evolation equations [2, 3], fractional order Zika virus model [4] and fractional Lagevin equation [5].

The bivariate ML function is a generalization of the classical ML function applied to two variables and has been employed to solve fundamental fractional differential equations involving two independent fractional orders. Recent studies include but not limited with the coupled-Laplacian fractional differential equations with nonlinear boundary conditions [6], singular multi-order fractional differential equations of Lane–Emden type [7]. In last quarter century there has been a growing interest in studying different variants of bivariate ML functions and fractional calculus operators with these functions in the kernel (see [8,9,10,11,12,13,14,15,16,17,18,19,20,21,22,23,24,25,26]).

This work focuses on the bivariate ML functions $E_{\alpha ,\beta ,\gamma }^{\delta }(x,y)$, introduced in [13] as a double series:

$$\begin{aligned} E_{\alpha ,\beta ,\gamma }^{\delta }(x,y)=\sum _{s=0}^\infty \sum _{r=0}^\infty \frac{(\delta )_{s+r}x^sy^r}{\Gamma (\alpha s+\beta r+\gamma )r!s!}, \end{aligned}$$

(1)

where $\alpha , \beta ,\gamma ,\delta$ are complex parameters with $Re(\alpha )>0$ and $Re(\beta )>0$, and the notation $(a)_s$ is used for the Pochhammer symbol $\frac{\Gamma (a+s)}{\Gamma (a)}$. The given series converges absolutely and locally uniformly for $Re(\alpha )>0$ and $Re(\beta )>0$. The authers of [13] gave the proof by using the convergence conditions studied in [27] for the generalised Lauricella series in two variables. If all parameters are 1 in Eq. (1) we recover the double exponential function which is the natural analogue of the fact that the ML function reduces to exponential function.

For $x>a$ the following fractional integral operator was also defined in [13], with the kernel containing (1)

$$\begin{aligned} \big ({}_{{a}}{I}_{\alpha ,\beta ,\gamma }^{\delta ,w_1,w_2}\psi \big )(x)= \int _a^x (x-t)^{\gamma -1} E_{\alpha ,\beta , \gamma }^{\delta }(w_1(x-t)^\alpha ,w_2(x-t)^\beta )\psi (t)dt, \end{aligned}$$

(2)

where $\psi (x) \in L^1(a,b)$, being the space of absolutely integrable functions. As shown in Theorem 8 of [13], the integral operator is bounded on the space $L^1(a,b)$. The extra restriction $Re(\gamma )>0$ is to avoid non-integrable singularity at the end point $t=x$. In the case $\gamma =0$, the fractional integral operator ${}_{{a}}I_{\alpha ,\beta ,\gamma }^{\delta ,w_1,w_2}$ reduces to the Riemann - Liouville fractional integral operator which is written as follows:

$$\begin{aligned} {}^{{RL}}_{{a}}I_{x}^{\mu }\psi (x)=\frac{1}{\Gamma (\mu )}\int _a^x (x-\epsilon )^{\mu -1}\psi (\epsilon )d\epsilon , \end{aligned}$$

where $x>a,Re(\mu )>0.$

Additionally the fractional derivatives are given by,

$$\begin{aligned} {}^{{RL}}_{{a}}D_{x}^{\mu }\psi (x) =\frac{d^n}{dx^n}\bigg ({}^{{RL}}_{{a}}I_{x}^{n-\mu }\psi (x)\bigg ), \quad n=\lfloor Re(\mu )\rfloor +1, \quad Re(\mu )\ge 0. \end{aligned}$$

Distinguishing itself from the Riemann–Liouville fractional derivative, the Caputo fractional derivative is defined based on the Riemann–Liouville fractional integral,

$$\begin{aligned} {}^{{C}}_{{a}}D_{x}^{\mu }\psi (x) ={}^{{RL}}_{{a}}I_{x}^{n-\mu }\frac{d^n}{dx^n}\psi (x), \quad n=\lfloor Re(\mu )\rfloor +1, \quad Re(\mu )\ge 0. \end{aligned}$$

In [13], the function in (1) and the fractional integral operator (2) were introduced and their properties were studied, including the crucial semigroup property of the fractional integral operator for any summable function $\psi \in L^1(a,b)$:

$$\begin{aligned} \begin{aligned} \big ({}_{{a}}{I}_{\alpha ,\beta ,\gamma _1}^{\delta _1,w_1,w_2}{}_{{a}}{I}_{\alpha , \beta ,\gamma _2}^{\delta _2,w_1,w_2}\psi \big )(x) =\big ({}_{{a}}{I}_{\alpha ,\beta ,\gamma _1+\gamma _2}^{\delta _1+\delta _2,w_1,w_2}\psi \big )(x). \end{aligned} \end{aligned}$$

(3)

By using the semigroup property the corresponding left inverse operator to the fractional integral operator (2) was defined in [13] as follows,

$$\begin{aligned} \big ({}^{{RL}}_{{a}}{D}_{\alpha ,\beta ,\gamma }^{\delta ,w_1,w_2}\psi \big )(x) ={}^{{RL}}_{{a}}D_{x}^{\gamma +\sigma }{}_{{a}}{I}_{\alpha ,\beta , \gamma }^{-\delta ,w_1,w_2}\psi (x)=\frac{d^n}{dx^n}\big ({}_{{a}}{I}_{\alpha ,\beta , n-\gamma }^{-\delta ,w_1,w_2}\psi (x)\big ), \end{aligned}$$

(4)

and the corresponding Caputo fractional derivative was defined by,

$$\begin{aligned} \big ({}^{{C}}_{{a}}{D}_{\alpha , \beta ,\gamma }^{\delta ,w_1,w_2}\psi \big )(x)={}_{{a}}{I}_{\alpha , \beta ,n-\gamma }^{-\delta ,w_1,w_2}\big (\frac{d^n}{dx^n}\psi (x)\big ), \quad \quad n=\lfloor Re(\gamma )\rfloor +1. \end{aligned}$$

(5)

In the present study, we undertake additional analysis of the integral operator (2) and derivative operators (4), (5), naturally extending the results [13]. In Sect. 2 we give Laplace transforms of the operators and the product rule for the operator (2). In Sect. 3 we investigate the eigenvalue problem for the Caputo version derivative and the solution of nonhomogenuous Cauchy problem involving Riemann–Liouville derivative and integral operator (2). Section 4 is devoted to the numerical approximation of the Caputo type derivative operator through the use of linear and quadratic univariate Lagrange interpolation. In Sect. 5, we extend and generalize the theory by using fractional calculus operators with respect to functions. In Sect. 6, we conclude the paper.

2 Fractional-calculus of the operators with bivariate ML kernel

In this section, we remember some properties of the fractional calculus operators from [13] and we give some new results for these operators. The following theoretical results are crucial for the derivation of the theoretical findings in the sebsequent sections.

For $\alpha , \beta ,\gamma , \delta , w_1, w_2 \in \mathbb {C}, Re(\beta )>0, Re(\alpha )>0, Re(\gamma )>0$ and for any function $\psi \in L^1( a,b)$, the fractional integral operator (2) can be written as [13, Theorem 7]

$$\begin{aligned} \big ({}_{{a}}{I}_{\alpha ,\beta ,\gamma }^{\delta ,w_1,w_2}\psi \big )(x) =\sum _{s=0}^\infty \sum _{r=0}^\infty \frac{(\delta )_{s+r}w_1^sw_2^r}{r!s!} \big ({}^{{RL}}_{{a}}I_{x}^{\alpha s+\beta r+\gamma }\psi \big )(x), \end{aligned}$$

(6)

where the series on the right hand side is locally uniformly convergent. To demonstrate the effectiveness of (6) we give the following resullts which show the compositon of the fractional operators ${}^{{RL}}_{{a}}I_{x}^{\sigma }$ and ${}^{{RL}}_{{a}}D_{x}^{\sigma }$ with the operator (2). The relation [13, Corollary 4]

$$\begin{aligned} {}^{{RL}}_{{a}}I_{x}^{\sigma }\bigg ({}_{{a}}{I}_{\alpha ,\beta , \gamma }^{\delta ,w_1,w_2}\psi (x)\bigg )={}_{{a}}{I}_{\alpha ,\beta , \gamma +\sigma }^{\delta ,w_1,w_2}\psi (x)={}_{{a}}{I}_{\alpha ,\beta , \gamma }^{\delta ,w_1,w_2}\bigg ({}_{a}I_{x}^{\sigma }\psi (x)\bigg ), \end{aligned}$$

holds true for any summable function $\psi \in L^1[a,b]$ and any $\alpha , \beta ,\gamma , \delta , w_1, w_2 \in \mathbb {C}$, $Re(\beta )>0, Re(\alpha )>0, Re(\gamma )>0.$ Also the relation

$$\begin{aligned} {}^{{RL}}_{{a}}D_{x}^{\sigma }\bigg ({}_{{a}}{I}_{\alpha ,\beta , \gamma }^{\delta ,w_1,w_2}\psi (x)\bigg )={}_{{a}}{I}_{\alpha ,\beta ,\gamma -\sigma }^{\delta ,w_1,w_2}\psi (x), \end{aligned}$$

holds true for any continuous function $\psi \in C[a,b]$ and any $\alpha , \beta ,\gamma , \delta , w_1, w_2 \in \mathbb {C}$, $Re(\beta )>0, Re(\alpha )>0, Re(\gamma )>0.$

In order to highlight the usefullness of the series formula in (6), we will now illustrate a straightforward approach for computing the outcome of applying the integral operator (2) to an elementary power function. This method offers a simplified alternative for using the original definition (2).

Proposition 1

The integral operator ${}_{{a}}{I}_{\alpha ,\beta ,\gamma }^{\delta ,w_1,w_2}$ on a power function is given as:

$$\begin{aligned} {}_{{a}}{I}_{\alpha ,\beta ,\gamma }^{\delta ,w_1,w_2}\big [(x-a)^\sigma \big ] =\Gamma (\sigma +1)(x-a)^{\sigma +\gamma }E_{\alpha ,\beta ,\gamma +\sigma }^{\delta }[w_1(x-a)^\alpha ,w_2(x-a)^\beta ]. \end{aligned}$$

(7)

Proof

According to Riemann–Liouville the fractional integral of a power function can be presented as

$$\begin{aligned} {}^{{RL}}_{{a}}I_{x}^{\alpha }(x-a)^\delta =\frac{\Gamma (\delta +1)}{\Gamma (\delta +\alpha +1)}(x-a)^{\delta +\alpha }, \quad Re(\delta )>-1, Re(\alpha )>0. \end{aligned}$$

By combining this established result with the series formula (6), we derive the following relation for $Re(\delta )>-1$

$$\begin{aligned} \begin{aligned}&{}_{{a}}{I}_{\alpha ,\beta ,\gamma }^{\delta ,w_1,w_2}\big [(x-a)^\sigma \big ]{} =\sum _{s=0}^\infty \sum _{r=0}^\infty \frac{(\delta )_{s+r}w_1^sw_2^r}{r!s!}\big ({}^{{RL}}_{{a}}I_{x}^{\alpha s+\beta r+\gamma }(x-a)^\sigma \big )\\&=\sum _{s=0}^\infty \sum _{r=0}^\infty \frac{(\delta )_{s+r}w_1^sw_2^r}{r!s!}\frac{\Gamma (\sigma +1)}{\Gamma (\sigma +\alpha s+\beta r+\gamma )}(x-a)^{\sigma +\alpha s+\beta r+\gamma } \\&=\Gamma (\sigma +1)(x-a)^{\sigma +\gamma } \sum _{s=0}^\infty \sum _{r=0}^\infty \frac{(\delta )_{s+r}w_1^sw_2^r}{\Gamma (\sigma +\alpha s+\beta r+\gamma )r!s!}(x-a)^{\alpha s+\beta r}. \end{aligned} \end{aligned}$$

$\square$

Example 1

In this example we calculate the integral operator $_{0}I_{\alpha ,\beta ,\gamma }^{\delta ,w_{1},w_{2}}\left[ x^{\sigma }\right]$ given in (7) for $w_{1}=w_{2}=1,$ over the interval $\Omega =\left[ 0,1\right]$ at the mesh points

$$\begin{aligned} \Omega ^{\zeta }=\left\{ x:x=x_{i}=i\zeta ,\text { }\zeta =\frac{1}{N_{1}},\text { } i=0,1,2,\ldots ,N_{1}\right\} . \end{aligned}$$

(8)

For $N_{1}=64$ and $\alpha =0.3,\beta =0.9$ and $\delta =\sigma =-0.7$ the integral $_{0}I_{0.3,0.9,\gamma }^{-0.7,1,1}\left[ x^{-0.7}\right]$ is presented in Fig. 1 for various values of $\gamma$ as 0.1, 0.5, 0.9. Further, the integral $_{0}I_{0.3,0.9,0.6}^{-0.7,1,1}\left[ x^{\sigma }\right]$ for $\sigma =0.5,0.9,1.9$ is illustrated in Fig. 2.

Fig. 1

The integral $_{0}I_{0.3,0.9,\gamma }^{-0.7,1,1}\left[ x^{-0.7}\right]$ for various values of $\gamma$

Fig. 2

The integral $_{0}I_{0.3,0.9,0.6}^{-0.7,1,1}\left[ x^{\sigma }\right]$ for various values of $\sigma$

Proposition 2

The integral operator ${}_{{a}}{I}_{\alpha ,\beta ,\gamma }^{\delta ,w_1,w_2}$ on an exponential function is given as:

$$\begin{aligned} {}_{{-\infty }}{I}_{\alpha ,\beta ,\gamma }^{\delta ,w_1,w_2}[e^{\sigma x}]=\sigma ^{-\gamma }e^{\sigma x}\bigg [1-\frac{w_1}{\sigma ^\alpha }-\frac{w_2}{\sigma ^\beta }\bigg ]^{-\delta }. \end{aligned}$$

(9)

Proof

The fractional integral of an exponential function according to Riemann–Liouville can be presented as

$$\begin{aligned} {}^{{RL}}_{{-\infty }}I_{y}^{\beta }e^{\sigma y}=\sigma ^{-\beta }e^{\sigma y}, \quad Re(\sigma )>0, Re(\beta )>0. \end{aligned}$$

Using this result with the series formula (6), we establish the following for $Re(\delta )>0$

$$\begin{aligned} \begin{aligned} {}_{{-\infty }}{I}_{\alpha ,\beta ,\gamma }^{\delta ,w_1,w_2}[e^{\sigma x}]{}&= \sum _{s=0}^\infty \sum _{r=0}^\infty \frac{(\delta )_{s+r}w_1^sw_2^r}{r!s!}\big ({}^{{RL}}_{{a}}I_{x}^{\alpha s+\beta r+\gamma }e^{\sigma x}\big ) \\&=\sigma ^{-\gamma }e^{\sigma x}\sum _{s=0}^\infty \frac{(\delta )_s}{s!}\bigg (\frac{w_1}{\sigma ^\alpha }\bigg )^s\sum _{r=0}^\infty \frac{(\delta +s)_r}{r!}\bigg (\frac{w_2}{\sigma ^\beta }\bigg )^{ r} \\&=\sigma ^{-\gamma }e^{\sigma x}\bigg [1-\frac{w_2}{\sigma ^\beta }\bigg ]^{-\delta }\sum _{s=0}^\infty \frac{(\delta )_s}{s!}\bigg (\frac{w_1}{\sigma ^\alpha }\bigg )^s\bigg (\frac{\sigma ^\beta }{\sigma ^\beta -w_2}\bigg )^s \\&=\sigma ^{-\gamma }e^{\sigma x}\bigg [1-\frac{w_2}{\sigma ^\beta }\bigg ]^{-\delta }\bigg [1-\bigg (\frac{w_1}{\sigma ^\alpha }\bigg )\bigg (\frac{\sigma ^\beta }{\sigma ^\beta -w_2}\bigg )\bigg ]^{-\delta }\\&=\sigma ^{-\gamma }e^{\sigma x}\bigg [1-\frac{w_1}{\sigma ^\alpha }-\frac{w_2}{\sigma ^\beta }\bigg ]^{-\delta }. \end{aligned} \end{aligned}$$

$\square$

Subsequently, our attention is directed towards the Laplace transform, with specific reference to [13], which briefly discussed the Laplace transform of the bivariate ML function (1). Because the gamma function on the denominator involves both s and r together, the property of formula (1) inables computation of the double Laplace transform with respect to x and y. Instead of taking the double Laplace transform they considered the univariate version of (1) and calculated the Laplace transform with respect to one variable t.

It was obtained in [13] that

$$\begin{aligned} \mathbb {L}\bigg [t^{\gamma -1}E_{\alpha ,\beta ,\gamma }^{\delta }(w_1t,w_2t) \bigg ](p)=\frac{1}{p^\gamma }\bigg (1-\frac{w_1}{p^\alpha }-\frac{w_2}{p^\beta }\bigg )^{-\delta } , \end{aligned}$$

(10)

where $Re(\beta )>0, Re(\alpha )>0, Re(\gamma )>0,Re(p)>0$ and $\mathbb {L}[f(x)](p)$ is the univariate Laplace transform operator defined by

$$\begin{aligned} \mathbb {L}[f(x)](p)=\int _0^\infty e^{-px}f(x)dx, \quad Re(p)>0. \end{aligned}$$

Now, we obtain the Laplace transforms of the fractional integral (2), Riemann–Liouville type and Caputo type fractional derivatives corresponding to the bivariate ML functions.

Theorem 3

The Laplace transform of integral operator (2) can be expressed by

$$\begin{aligned} \mathbb {L}\bigg [{}_{{0}}{I}_{\alpha ,\beta ,\gamma }^{\delta ,w_1,w_2}\psi (x)\bigg ](p) =\frac{\mathbb {L}\big [\psi (x)\big ](p)}{p^\gamma }\bigg (1-\frac{w_1}{p^\alpha }-\frac{w_2}{p^\beta }\bigg )^{-\delta }, \end{aligned}$$

(11)

where $Re(\beta )>0, Re(\alpha )>0, Re(\gamma )>0$ and p is the Laplace transformed parameter satisfying $Re(p)>0.$

Proof

By using series formula (6) we get,

$$\begin{aligned} \begin{aligned} \mathbb {L}\bigg [{}_{{0}}{I}_{\alpha ,\beta , \gamma }^{\delta ,w_1,w_2}\psi (x)\bigg ](p){}&=\mathbb {L}\bigg [\sum _{s=0}^\infty \sum _{r=0}^\infty \frac{(\delta )_{s+r}w_1^sw_2^r}{r!s!}\big ({}^{{RL}}_{{0}}I_{x}^{\alpha s+\beta r+\gamma }\psi (x)\big )\bigg ](p)\\ {}&=\sum _{s=0}^\infty \sum _{r=0}^\infty \frac{(\delta )_{s+r} w_1^sw_2^r}{r!s!}\mathbb {L}\bigg [\big ({}^{{RL}}_{{0}}I_{x}^{\alpha s+\beta r+\gamma }\psi (x)\big )\bigg ](p). \end{aligned} \end{aligned}$$

The series manipulation is same as in the proof of Proposition 2 and by using the fact [26] that Laplace transforms of Riemann–Liouville integrals,

$$\begin{aligned} \mathbb {L}\bigg [{}_{y}I_{0^+}^{\beta }f(x)\bigg ](p)=\frac{1}{p^\beta }\mathbb {L}[f(x)](p), \quad Re(\beta )>0, Re(p)>0, \end{aligned}$$

we obtain (11). $\square$

Theorem 4

The Laplace transform of Riemann–Liouville type fractional derivative operator ${}^{{RL}}_{{0}}{D}_{\alpha ,\beta ,\gamma }^{\delta ,w_1,w_2}$ can be expressed as below

$$\begin{aligned} \begin{aligned} \mathbb {L} {}&\big [{}^{{RL}}_{{0}}{D}_{\alpha ,\beta ,\gamma }^{\delta ,w_1,w_2}\psi (x)\big ](p) =p^\gamma \bigg (1-\frac{w_1}{s^\alpha }-\frac{w_2}{s^\beta }\bigg )^\delta \mathbb {L}[\psi (x)](p) \\ {}&-\sum _{s=0}^{\infty }\sum _{r=0}^{\infty }\frac{(-\delta )_{s+r}w_1^sw_2^r}{s!r!}\sum _{m=0}^n p^{m-1}\frac{d^{n-m}}{dx^{n-m}}\bigg ({}^{{RL}}_{{0}}I_{x}^{\alpha s+\beta r-\gamma +n}\psi \bigg )(0), \end{aligned} \end{aligned}$$

where $Re(\beta )>0, Re(\alpha )>0, Re(\gamma )>0$ and p is the Laplace transform parameter satisfying $Re(p)>0$ and $n-1 \ne Re(\gamma )< n$, $n \in \mathbb {N}.$

Proof

By using the series representation for the operator ${}^{{RL}}_{{0}}{D}_{\alpha ,\beta ,\gamma }^{\delta ,w_1,w_2}\psi (x)$, we get

$$\begin{aligned} \begin{aligned} \mathbb {L}\big [{}^{{RL}}_{{0}}{D}_{\alpha ,\beta ,\gamma }^{\delta ,w_1,w_2}\psi (x)\big ](p) {}&=\mathbb {L}\bigg [\frac{d^n}{dx^n}\big ({}_{{0}}{I}_{\alpha ,\beta ,\gamma }^{-\delta ,w_1,w_2}\psi (x)\big )\bigg ](p) \\&=\mathbb {L}\bigg [\frac{d^n}{dx^n}\sum _{s=0}^\infty \sum _{r=0}^\infty \frac{(-\delta )_{s+r}w_1^sw_2^r}{r!s!}\big ({}^{{RL}}_{{0}}I_{x}^{\alpha s+\beta r-\gamma +n}\psi (x)\big ) \bigg ](p). \end{aligned} \end{aligned}$$

Based on the locally uniformly convergent property of the series and applying the known formula for the Laplace transform of the $n-$th derivative of the function, we obtain

$$\begin{aligned} \begin{aligned} {}&\mathbb {L} \big [{}^{{RL}}_{{0}}{D}_{\alpha ,\beta ,\gamma }^{\delta ,w_1,w_2}\psi (x)\big ](p)\\&=p^n\sum _{s=0}^\infty \sum _{r=0}^\infty \frac{(-\delta )_{s+r}w_1^sw_2^r}{r!s!}\mathbb {L} \bigg [\big ({}^{{RL}}_{{0}}I_{x}^{\alpha s+\beta r-\gamma +n}\psi (x)\big )\bigg ](p) \\&\quad -\sum _{s=0}^{\infty }\sum _{r=0}^{\infty }\frac{(-\delta )_{s+r}w_1^sw_2^r}{s!r!}\sum _{m=0}^n p^{m-1}\frac{d^{n-m}}{dx^{n-m}}\bigg ({}^{{RL}}_{{0}}I_{x}^{\alpha s+\beta r-\gamma +n}\psi (x) \bigg )(0)\\&=p^n\sum _{s=0}^\infty \sum _{r=0}^\infty \frac{(-\delta )_{s+r}w_1^sw_2^r}{r!s!}p^{-(\alpha s+\beta r -\gamma +n)}\mathbb {L}[\psi (x)](p) \\&\quad -\sum _{s=0}^{\infty }\sum _{r=0}^{\infty }\frac{(-\delta )_{s+r}w_1^sw_2^r}{s!r!}\sum _{m=0}^n p^{m-1}\frac{d^{n-m}}{dx^{n-m}}\bigg ({}^{{RL}}_{{0}}I_{x}^{\alpha s+\beta r-\gamma +n}\psi (x)\bigg )(0) \\&=p^\gamma \bigg (1-\frac{w_1}{s^\alpha }-\frac{w_2}{s^\beta }\bigg )^\delta \mathbb {L}[\psi (x)](p)\\&\quad -\sum _{s=0}^{\infty }\sum _{r=0}^{\infty }\frac{(-\delta )_{s+r}w_1^sw_2^r}{s!r!}\sum _{m=0}^n p^{m-1}\frac{d^{n-m}}{dx^{n-m}}\bigg ({}^{{RL}}_{{0}}I_{x}^{\alpha s+\beta r-\gamma +n}\psi (x)\bigg )(0). \end{aligned} \end{aligned}$$

$\square$

Theorem 5

The Laplace transform of the Caputo type fractional derivative ${}^{{C}}_{{0}}{D}_{\alpha ,\beta ,\gamma }^{\delta ,w_1,w_2}$ can be given as follows

$$\begin{aligned} \begin{aligned} \mathbb {L}\big [{}^{{C}}_{{0}}{D}_{\alpha , \beta ,\gamma }^{\delta ,w_1,w_2}\psi (x)\big ](p)=p^{\gamma -n} \bigg (1-\frac{w_1}{s^\alpha }-\frac{w_2}{s^\beta }\bigg )^\delta \\ \bigg [p^n\mathbb {L}[\psi (x)](p)-\sum _{m=1}^n p^{n-m}\frac{d^{m-1}}{dx^{m-1}}\psi (0)\bigg ], \end{aligned} \end{aligned}$$

where the conditions on the parameters are same as in Theorem 4.

Proof

By using Theroem 3 and applying formula for the Laplace transform of the $n-$th derivative of the function, we obtain

$$\begin{aligned}&\mathbb {L}\big [{}^{{C}}_{{0}}{D}_{\alpha ,\beta ,\gamma }^{\delta ,w_1,w_2}\psi (x)\big ](p) = \mathbb {L}\bigg [{}_{{0}}{I}_{\alpha ,\beta ,\gamma }^{-\delta ,w_1,w_2}\big (\frac{d^n}{dx^n}\psi (x)\big )\bigg ](p)\\&=p^{-n+\gamma }\bigg (1-\frac{w_1}{p^\alpha }-\frac{w_2}{p^\beta }\bigg )^{\delta }\mathbb {L}\bigg [\psi (x)\bigg ] (p)\\&=p^{\gamma -n}\bigg (1-\frac{w_1}{s^\alpha }-\frac{w_2}{s^\beta }\bigg )^\delta \bigg [p^n\mathbb {L}[\psi (x)](p)-\sum _{m=1}^n p^{n-m}\frac{d^{m-1}}{dx^{m-1}}\psi (0)\bigg ], \end{aligned}$$

which gives the desired result. $\square$

In the literature the product rule is extensively studied and generalized to the fractional scenarios [28]. For the integral operator (2) using the series formula (6) we give a version of the product rule.

Theorem 6

Let f and g be complex functions such that f,g and f(x)g(x) are all in the form $x^\eta \Delta (x)$ with $Re(\eta )>0$ and $\Delta$ holomorphic on a complex domain $U \subset \mathbb {C}$. Under the conditions $Re(\beta )>0, Re(\alpha )>0, Re(\gamma )>0$, the integral operator (2) satisfies the following product rule

$$\begin{aligned} \begin{aligned}&{}_{{a}}{I}_{\alpha ,\beta ,\gamma }^{\delta ,w_1,w_2} \big (f(x)g(x)\big ){}\\&=\sum _{n=0}^\infty \frac{1}{n!}\frac{d^{n}g(x)}{dx^n}\bigg [\sum _{s=0}^\infty \sum _{r=0}^\infty \frac{\Gamma (\delta +s+r)\Gamma (1-\alpha s- \beta r-\gamma )w_1^sw_2^r}{\Gamma (\delta )\Gamma (1-\alpha s- \beta r-\gamma -n)s!r!} \\&\quad \times {}_{x}I_{a^+}^{\alpha s+\beta r+\gamma +n}f(x)\bigg ] . \end{aligned} \end{aligned}$$

(12)

Proof

By using series formula (6) with the known results [29, 30] on the product rule for Riemann–Liouville fractional differintegrals we get the following

$$\begin{aligned} \begin{aligned} {}&{}_{{a}}{I}_{\alpha ,\beta ,\gamma }^{\delta ,w_1,w_2}\big (f(x)g(x)\big ) =\sum _{s=0}^\infty \sum _{r=0}^\infty \frac{(\delta )_{s+r}w_1^sw_2^r}{r!s!}\big [{}^{{RL}}_{{a}}I_{x}^{\alpha s+\beta r+\gamma }\big (f(x)g(x)\big ] \\&\quad =\sum _{s=0}^\infty \sum _{r=0}^\infty \frac{(\delta )_{s+r}w_1^sw_2^r}{r!s!}\sum _{n=0}^\infty \left( {\begin{array}{c}-\alpha s- \beta r-\gamma \\ n\end{array}}\right) {}^{{RL}}_{{a}}I_{x}^{\alpha s+\beta r+\gamma +n}f(x)\frac{d^{n}g(x)}{dx^n} \\&\quad =\sum _{n=0}^\infty \frac{d^{n}g(x)}{dx^n} \bigg [\sum _{s=0}^\infty \sum _{r=0}^\infty \frac{(\delta )_{s+r}w_1^sw_2^r}{r!s!}\left( {\begin{array}{c}-\alpha s- \beta r-\gamma \\ n\end{array}}\right) {}^{{RL}}_{{a}}I_{x}^{\alpha s+\beta r+\gamma +n}f(x) \bigg ], \end{aligned} \end{aligned}$$

given (12). $\square$

Example 2

To verify the results of Theorem 6, we take the functions $f(x)=e^{cx}$ and $g(x)=x$ and the constant of differintegration $a=i\infty$. Consequently, the outer series in (12) has only two non-trivial terms. Thus;

$$\begin{aligned}&\sum _{n=0}^1\frac{1}{n!}\frac{d^{n}x}{dx^n}\bigg [\sum _{s=0}^\infty \sum _{r=0}^\infty \frac{\Gamma (\delta +s+r)\Gamma (1-\alpha s- \beta r-\gamma )w_1^sw_2^r}{\Gamma (\delta )\Gamma (1-\alpha s- \beta r-\gamma -n)s!r!}{}_{x}I_{a^+}^{\alpha s+\beta r+\gamma +n}e^{cx}\bigg ] \\&\quad =\sum _{n=0}^1\frac{1}{n!}\frac{d^{n}x}{dx^n}\bigg [\sum _{s=0}^\infty \sum _{r=0}^\infty \frac{\Gamma (\delta +s+r)\Gamma (1-\alpha s- \beta r-\gamma )w_1^sw_2^r}{\Gamma (\delta )\Gamma (1-\alpha s- \beta r-\gamma -n)s!r!}c^{\alpha s+ \beta r+\gamma +1}e^{cx}\bigg ] \\&\quad =x\bigg [\sum _{s=0}^\infty \sum _{r=0}^\infty \frac{\Gamma (\delta +s+r)\Gamma (1-\alpha s- \beta r-\gamma )(c^\alpha w_1)^s(c^\beta w_2)^r}{\Gamma (\delta )\Gamma (1-\alpha s- \beta r-\gamma )s!r!}c^\gamma e^{cx}\bigg ] \\&\qquad +\bigg [\sum _{s=0}^\infty \sum _{r=0}^\infty \frac{\Gamma (\delta +s+r)\Gamma (1-\alpha s- \beta r-\gamma )(c^\alpha w_1)^s(c^\beta w_2)^r}{\Gamma (\delta )\Gamma (-\alpha s- \beta r-\gamma )s!r!}c^{\gamma +1} e^{cx}\bigg ] \\&\quad =x\bigg [\sum _{s=0}^\infty \sum _{r=0}^\infty \frac{\Gamma (\delta +s+r)(c^\alpha w_1)^s(c^\beta w_2)^r}{\Gamma (\delta )s!r!}c^\gamma e^{cx}\bigg ]\\ {}&\qquad +\bigg [\sum _{s=0}^\infty \sum _{r=0}^\infty \frac{\Gamma (\delta +s+r)(-\alpha s- \beta r-\gamma )(c^\alpha w_1)^s(c^\beta w_2)^r}{\Gamma (\delta )s!r!}c^{\gamma +1} e^{cx}\bigg ] \\ {}&\quad =\sum _{s=0}^\infty \sum _{r=0}^\infty \frac{\Gamma (\delta +s+r)(c^\alpha w_1)^s(c^\beta w_2)^r}{\Gamma (\delta )s!r!}c^\gamma e^{cx}(x-c(\alpha s+\beta r+\gamma ), \end{aligned}$$

hence we have computed the integral ${}_{{i\infty }}{I}_{\alpha ,\beta ,\gamma }^{\delta ,w_1,w_2}\big (xe^{cx}\big )$.

The chain rule in classical calculus can be generalized to certain models of fractional calculus. In the following theorem for this particular calculus model we use the proposed series formula (6) to derive a result of the chain rule.

Theorem 7

Let f and g be complex functions such that g is smooth and f(g(x)) is a function of the form $x^\eta \Delta (x)$ with $Re(\eta )>0$ and $\Delta$ holomorphic on a complex domain $U \subset \mathbb {C}$. Under the conditions $Re(\beta )>0, Re(\alpha )>0, Re(\kappa )>0$, the integral operator (2) satisfies the following chain rule.

$$\begin{aligned} \begin{aligned} {}_{{a}}{I}_{\alpha ,\beta ,\gamma }^{\delta ,w_1,w_2}\big (f(g(x))\big ){}&=(x-a)^\gamma \sum _{s=0}^\infty \sum _{r=0}^\infty \frac{\Gamma (\delta +s+r)[w_1(x-a)^\alpha ]^s[w_2(x-a)^\beta ]^{ r}}{\Gamma (\delta )\Gamma (\alpha s+\beta r+\gamma )s!r!} \\ {}&\quad \times \sum _{n=0}^\infty \frac{(a-x)^n}{n!(\beta r+\alpha s+n)}\sum _{r=1}^{m}\frac{d^rf(g(x))}{dg(x)^r} \\ {}&\quad \times \sum _{(P_1,\dots ,P_m)}\prod _{j=1}^{n}\frac{1}{P_j!(j!)^{P_j}}\bigg (\frac{d^{j}g(x)}{d x^{j}}\bigg )^{P_j}, \end{aligned} \end{aligned}$$

(13)

where the summation $(P_1,\dots ,P_m)$ is over the set

$$\begin{aligned} \bigg \{(P_1,\dots ,P_m) \in \mathbb (Z_0^+)^m:\sum _jP_j=r, \sum _j jP_j=m\bigg \}. \end{aligned}$$

Proof

Applying Theorem 6 to the product of $\big (f(g(x))\big )$ and $I(x)=1$, where the Riemann–Liouville fractional differintegrals of I are well known, yields (12),

$$\begin{aligned} \begin{aligned} {}&{}_{{a}}{I}_{\alpha ,\beta ,\gamma }^{\delta ,w_1,w_2}\big (f(g(x)\big ) \\ {}&\quad =\sum _{n=0}^\infty \frac{d^{n}f(g(x))}{n!dx^n}\bigg [\sum _{s=0}^\infty \sum _{r=0}^\infty \frac{\Gamma (\delta +s+r)\Gamma (1-\alpha s- \beta r-\gamma )w_1^sw_2^r}{\Gamma (\delta )\Gamma (1-\alpha s- \beta r-\gamma -n)s!r!} {}_{x}I_{a^+}^{\alpha s+\beta r+\gamma +n}(1)\bigg ] \\ {}&\quad =\frac{d^{n}f(g(x))}{n!dx^n}\bigg [\sum _{s=0}^\infty \sum _{r=0}^\infty \frac{\Gamma (\delta +s+r)\Gamma (1-\alpha s- \beta r-\gamma )w_1^sw_2^r}{\Gamma (\delta )\Gamma (1-\alpha s- \beta r-\gamma -n)s!r!} \frac{(x-a)^{\alpha s+\beta r+\gamma +n}}{\Gamma (\alpha s+\beta r+\gamma +n+1)}\bigg ]. \end{aligned} \end{aligned}$$

We use the reflection formula in order to eliminate some of the gamma functions in this expression,

$$\begin{aligned} \Gamma (\alpha s+\beta r+\gamma +n+1)\Gamma (1-\alpha s- \beta r-\gamma -n) =(-1)^n(\alpha s+\beta r+\gamma +n)\Gamma (\alpha s+\beta r+\gamma )\Gamma (1-\alpha s-\beta r-\gamma ). \end{aligned}$$

$$\begin{aligned} \begin{aligned} {}&{}_{{a}}{I}_{\alpha ,\beta ,\gamma }^{\delta ,w_1,w_2}\big (f(g(x)\big )\\&\quad =(x-a)^\gamma \sum _{n=0}^\infty \frac{d^{n}f(g(x))}{dx^n}\bigg [\sum _{s=0}^\infty \sum _{r=0}^\infty \frac{\Gamma (\delta +s+r)[w_1(x-a)^\alpha ]^s[w_2(x-a)^\beta ]^r}{\Gamma (\delta )(\alpha s+\beta r+\gamma +n)\Gamma (\alpha s+\beta r+\gamma )s!r!} \bigg ]\frac{(a-x)^n}{n!}\\ {}&\quad =(x-a)^\gamma \bigg [\sum _{s=0}^\infty \sum _{r=0}^\infty \frac{\Gamma (\delta +s+r)[w_1(x-a)^\alpha ]^s[w_2(x-a)^\beta ]^r}{\Gamma (\delta )\Gamma (\alpha s+\beta r+\gamma )s!r!} \bigg ] \\ {}&\qquad \times \sum _{n=0}^\infty \frac{d^{n}f(g(x))}{dx^n}\frac{(a-x)^n}{n!(\alpha s+\beta r+\gamma +n)}. \end{aligned} \end{aligned}$$

By use of classical Faa di Bruno formula we get the result (13). $\square$

3 Eigenvalue problem and Cauchy problem for fractional integro-differential equation

Our focus now shifts to the eigenvalue problem corresponding to the Caputo fractional derivative (5). To initiate our investigation, we analyze the following Cauchy problem

$$\begin{aligned}&\big ({}^{{C}}_{{0}}{D}_{\alpha ,\beta ,\gamma }^{\delta ,w_1,w_2}\psi \big )(x)=C_1\psi (x), \nonumber \\&\psi ^{(j)}(0)=\delta _j, \quad (j=0,1,\dots ,n-1) \end{aligned}$$

(14)

where $\alpha , \beta ,\gamma , \delta \in \mathbb {C}, Re(\beta )>0, Re(\alpha )>0, Re(\gamma )>0$ and $n-1\le Re(\gamma )<n$ $(n \in \mathbb {N})$. $C_1, \delta _0 \dots ,\delta _{n-1} \in \mathbb {C}$. Applying Laplace transform to both sides of (14), we obtain

$$\begin{aligned} \begin{aligned} \mathbb {L}\big [{}^{{C}}_{{0}}{D}_{\alpha ,\beta ,\gamma }^{\delta ,w_1,w_2}\psi \big ](p)= {}&p^{\gamma -m}\bigg (1-\frac{w_1}{s^\alpha }-\frac{w_2}{s^\beta }\bigg )^\delta \bigg [p^m\mathbb {L}[\psi (x)](p)-\sum _{k=0}^{n-1} p^{n-k}\delta _k\bigg ] \\ {}&=C_1\mathbb {L}[\psi ](p). \end{aligned} \end{aligned}$$

This leads to

$$\begin{aligned} \begin{aligned} \mathbb {L}[\psi (x)](p) {}&=\frac{\sum _{k=0}^{n-1} p^{n-k}\delta _k}{1-C_1p^{-\gamma }\bigg (1-\frac{w_1}{s^\alpha }-\frac{w_2}{s^\beta }\bigg )^{-\delta }} \\ {}&=\sum _{k=0}^{n-1} p^{-k}\delta _k\sum _{l=0}^\infty C_1^l p^{-\gamma l}\bigg (1-\frac{w_1}{s^\alpha }-\frac{w_2}{s^\beta }\bigg )^{-\delta l} \\ {}&=\sum _{k=0}^{n-1} \delta _k\sum _{l=0}^\infty C_1^l p^{-\gamma l-k}\bigg (1-\frac{w_1}{s^\alpha }-\frac{w_2}{s^\beta }\bigg )^{-\delta l}, \end{aligned} \end{aligned}$$

given that $\bigg |C_1p^{-\gamma }\bigg (1-\frac{w_1}{s^\alpha }-\frac{w_2}{s^\beta }\bigg )^{-\delta }\bigg | < 1.$ After applying inverse Laplace transform to both sides and making use of Eq. (10), we arrive at

$$\begin{aligned} \psi (x)=\sum _{k=0}^{n-1}\sum _{l=0}^\infty \delta _k C_1^l x^{\gamma n+k-1}E_{\alpha ,\beta ,\gamma n+k}^{\delta n}(w_1x^\alpha ,w_2x^\beta ), \end{aligned}$$

which is the eigenfunction of the ${}^{{C}}_{{0}}{D}_{\alpha ,\beta ,\gamma }^{\delta ,w_1,w_2}$ corresponding to the eigenvalue $C_1$.

Our attention now turns to the following Cauchy-type problem:

$$\begin{aligned}&{}_{x}D_{a^+}^{\sigma }\psi (x)= {}_{{a}}{I}_{\alpha ,\beta ,\gamma }^{\delta ,w_1,w_2}\psi (x)+f(x), \end{aligned}$$

(15)

$$\begin{aligned}&\underset{x \rightarrow a^+}{\lim }\frac{d^n}{dx^n}\big ({}_{x}I_{a^+}^{k-\sigma }\big )\psi (x)=\eta _n \quad (n=0,1,\dots ,k-1), \end{aligned}$$

(16)

where $k-1 <\sigma \le k$ ($k \in \mathbb {N}$) and $\delta ,\alpha ,\beta ,\gamma , w_1,w_2 \in \mathbb {R}$ with $\alpha ,\beta >0$.

Theorem 8

The problem (15) and (16) has solution $\psi \in L_1(0,\infty )$ as:

$$\begin{aligned} \begin{aligned} \psi (x)&= {}\sum _{n=0}^{k-1}\eta _n x^{n+\beta -k}\sum _{r=0}^{\infty } x^{(\gamma +\sigma )r}E_{\alpha ,\beta ,(\gamma +\sigma )r-k+n+1}^{\delta r}(w_1x^\alpha ,w_2x^\beta ) \\ {}&\quad +\sum _{r=0}^\infty \big ({}_{{a}}{I}_{\alpha ,\beta ,(\gamma +\sigma )r+\beta }^{\delta r,w_1,w_2}f\big )(x). \end{aligned} \end{aligned}$$

(17)

Proof

Through the application of the Laplace transform to both sides of Eq. (15), we obtain

$$\begin{aligned} s^\sigma \mathbb {L}[\psi (x)](s)-\sum _{n=0}^{k-1}\eta _n s^{k-n-1}=s^{-\gamma }\bigg (1-\frac{w_1}{s^\alpha }-\frac{w_2}{s^\beta }\bigg )^{-\gamma }\mathbb {L}[\psi (x)](s)+\mathbb {L}[f(x)](s). \end{aligned}$$

Consequently,

$$\begin{aligned} \begin{aligned} \mathbb {L}[\psi (x)](s)&{}= \frac{\sum _{n=0}^{k-1}\eta _n s^{k-n-1-\sigma }+\frac{\mathbb {L}[f(x)](s)}{s^\beta }}{1-s^{-\gamma -\sigma } \bigg (1-\frac{w_1}{s^\alpha }-\frac{w_2}{s^\beta }\bigg )^{-\gamma }} \\ {}&=\bigg [\sum _{n=0}^{k-1}\eta _n s^{k-n-1-\sigma }+\frac{\mathbb {L}[f(x)](s)}{s^\beta }\bigg ]\sum _{r=0}^\infty \bigg [ s^{-\gamma -\sigma }\bigg (1-\frac{w_1}{s^\alpha }-\frac{w_2}{s^\beta }\bigg )^{-\gamma }\bigg ]^r \\ {}&=\sum _{n=0}^{k-1}\sum _{r=0}^\infty \eta _n s^{k-n-1-\sigma -(\gamma +\beta )r}\bigg (1-\frac{w_1}{s^\alpha }-\frac{w_2}{s^\beta }\bigg )^{-\gamma r} \\ {}&\quad +\sum _{r=0}^\infty \ s^{-\gamma r-\beta r -\sigma }\bigg (1-\frac{w_1}{s^\alpha }-\frac{w_2}{s^\beta }\bigg )^{-\gamma r}\mathbb {L}[f(x)](s). \end{aligned} \end{aligned}$$

Applying inverse Laplace transform yields the result (17). $\square$

Example 3

Consider the following nonhomogenous Cauchy-type problem:

$$\begin{aligned}&{}_{x}D_{a^+}^{\sigma }\psi (x)= {}_{{a}}{I}_{\alpha ,\beta ,\gamma }^{\delta _1,w_1,w_2}\psi (x)+x^{\mu -1} E_{\alpha ,\beta ,\mu }^{\delta _2}(w_1x^\alpha ,w_2x^\beta ), \nonumber \\&\underset{x \rightarrow a^+}{\lim }\frac{d^n}{dx^n}\big ({}_{x}I_{a^+}^{k-\sigma }\big )\psi (x)=(-1)^n \quad n=0,1\dots ,k-1, \end{aligned}$$

(18)

where $k-1 <\sigma \le k$ ($k \in \mathbb {N}$) and $\delta , \alpha ,\beta ,\mu ,\gamma , w_1,w_2 \in \mathbb {R}$ with $\alpha ,\beta , \mu >0$. Considering (2) and the identity given in Theorem 6 of [13] gives,

$$\begin{aligned} \begin{aligned} {}&\int _0^x(x-\epsilon )^{\gamma r+\sigma r+\beta -1}E_{\alpha ,\beta ,(\gamma +\sigma )r+\beta }^{\delta _1r}(w_1(x-\epsilon )^\alpha ,\\&w_2(x-\epsilon )^\beta )\epsilon ^{\mu -1}E_{\alpha ,\beta ,\mu }^{\delta _2}(w_1\epsilon ^\alpha ,w_2\epsilon ^\beta )d\epsilon \\ {}&=x^{\gamma r+\sigma r+\beta +\mu -1}E_{\alpha ,\beta ,(\gamma +\sigma )r+\beta +\mu }^{\delta _1 r+\delta _2}(w_1x^\alpha ,w_2x^\beta ). \end{aligned} \end{aligned}$$

From Theorem 8 it follows that the problem (18) has solution belonging to $L_1(0,\infty )$ given by

$$\begin{aligned} \begin{aligned} \psi (x){}&=\sum _{n=0}^{k-1}(-1)^n x^{n+\beta -k}\sum _{r=0}^{\infty } x^{(\gamma +\sigma )r}E_{\alpha ,\beta ,(\gamma +\sigma )r-k+n+1}^{\delta r}(w_1x^\alpha ,w_2x^\beta ) \\ {}&+x^{+\beta +\mu -1}\sum _{r=0}^\infty x^{\gamma r+\sigma r}E_{\alpha ,\beta ,(\gamma +\sigma )r+\beta +\mu }^{\delta _1 r+\delta _2}(w_1x^\alpha ,w_2x^\beta ). \end{aligned} \end{aligned}$$

4 Approximation of Caputo type derivative

Numerical approximations of Caputo–Prabhakar derivative including 3 parameters ML function was given in [31]. In this section, we give numerical approximations to Caputo operator containing bivariate ML in the kernel.

For any $t\in \mathbb {R}^+$ and $\mu \in \mathbb {C}$,

$$\begin{aligned} \int _0^t v^{\gamma -1}E_{\alpha ,\beta ,\gamma }^{\delta }(\mu _1 v^\alpha ,\mu _2 v^\beta )dv= t^\gamma E_{\alpha ,\beta ,\gamma +1}^{\delta }(\mu _1 t^\alpha ,\mu _2 t^\alpha ). \end{aligned}$$

(19)

For any $n \in \mathbb {N}$,

$$\begin{aligned} \frac{d^n}{dt^n}\bigg [t^{\gamma -1}E_{\alpha ,\beta ,\gamma }^{\delta }(\mu _1 t^\alpha ,\mu _2 t^\beta )\bigg ]= t^{\gamma -1-n} E_{\alpha ,\gamma -n}^{\delta }(\mu _1 t^\alpha ,\mu _2 t^\beta ). \end{aligned}$$

(20)

The interval [0, T] is divided into $N_1$ subintervals, each of length $\zeta =\frac{T}{N_1}$ with points $0=t_0<t_1<t_2<\dots <t_{N_1}=T$, where $t_i=i\zeta$, $i=0,1,\dots ,N_1$. We approximate the function g(t) by using first degree Lagrange interpolation function. For two points $(t_{i-1},g(t_{i-1}))$ and $(t_{i},f(t_{i}))$ linear interpolation function $P_1(t)$ is defined as,

$$\begin{aligned} P_1(t) =\frac{t_i-t}{t_i-t_{i-1}}g(t_{i-1})+\frac{t-t_{i-1}}{t_i-t_{i-1}}g(t_{i}). \end{aligned}$$

(21)

Also,

$$\begin{aligned} \frac{d}{dt}P_1(t) =\frac{g(t_{i})-g(t_{i-1})}{\zeta } \end{aligned}$$

(22)

and

$$\begin{aligned} g(t)-P_1(t)=\frac{f''(\epsilon _i)}{2}(t-t_i)(t-t_{i-1}), \end{aligned}$$

(23)

where $\epsilon _i\in (t_{i-1},t_i)$.

For $0<\gamma <1$ and $0<t<T$ by choosing $n=1$, the Caputo-type derivative (5) is defined as,

$$\begin{aligned} {}^{{C}}_{{0}}{D}_{\alpha ,\beta ,\gamma }^{\delta ,w_1,w_2}g(t)=\int _0^t (t-r)^{-\gamma } E_{\alpha ,\beta ,1- \gamma }^{-\delta }(w_1(t-r)^\alpha ,w_2(t-r)^\beta )g'(r)dr \end{aligned}$$

(24)

at points $t=t_l$, $l=1,\dots ,N_1$,

$$\begin{aligned} \begin{aligned} {}&\big ({}^{{C}}_{{0}}{D}_{\alpha ,\beta , \gamma ,w_1,w_2}^{\delta ,(1,N_1)}g\big )(t_l)\\&\quad \approx \int _0^{t_l}(t_l-r)^{-\gamma }E_{\alpha ,\beta ,1-\gamma }^{-\delta }(w_1(t_l-r)^\alpha ,w_2(t_l-r)^\beta )P_1'(r)dr\\&\quad =\sum _{i=1}^l \int _{t_{i-1}}^{t_i}(t_l-r)^{-\gamma }E_{\alpha ,\beta ,1-\gamma }^{-\delta }(w_1(t_l-r)^\alpha ,w_2(t_l-r)^\beta )P_1'(r)dr\\&\quad =\sum _{i=1}^l \int _{t_{i-1}}^{t_i}(t_l-r)^{-\gamma }E_{\alpha ,\beta ,1-\gamma }^{-\delta }(w_1(t_l-r)^\alpha ,w_2(t_l-r)^\beta )\frac{g(t_i)-g(t_{i-1})}{\zeta }dr\\&\quad =\zeta ^{-\gamma }\sum _{i=1}^l w_{1,l-i}g(t_i)- w_{1,l-i}g(t_{i-1}), \end{aligned} \end{aligned}$$

where

$$\begin{aligned} \begin{aligned} w_{1,l-i} {}&=(l-i+1)^{1-\gamma }E_{\alpha ,\beta ,2-\gamma }^{-\delta }(w_1(t_l-t_{i-1})^\alpha ,w_2(t_l-t_{i-1})^\beta )\\&\quad -(l-i)^{1-\gamma }E_{\alpha ,\beta ,2-\gamma }^{-\delta }(w_1(t_l-t_{i})^\alpha ,w_2(t_l-t_{i})^\beta ). \end{aligned} \end{aligned}$$

In the next theorem, we determine error bound for the given approximation.

Theorem 9

Let $g(t) \in C^2[0,T]$ for any $0<\gamma <1$, the truncation error $R(g,\zeta ,\gamma )$ satiesfies the following inequality,

$$\begin{aligned} R_1(g,\zeta ,\gamma ) \le \frac{ \zeta ^{2-\gamma }}{8} \underset{0 \le t \le t_l}{\max }\big |g''(t)\big | K_1, \end{aligned}$$

(25)

where $K_1$ is a positive constant given as

$$\begin{aligned} \begin{aligned} K _1 {}&= \bigg | (l-i+1)^{-\gamma }E_{\alpha ,\beta ,1-\gamma }^{-\delta }(w_1(t_l-t_{i-1})^\alpha ,w_2(t_l-t_{i-1})^\beta )\\&\quad -(l-i)^{-\gamma }E_{\alpha ,\beta ,1-\gamma }^{-\delta }(w_1(t_l-t_{i})^\alpha ,w_2(t_l-t_{i})^\beta ) \bigg |. \end{aligned} \end{aligned}$$

Proof

From (23) it follows that

$$\begin{aligned} \begin{aligned}&R_1(g,\zeta ,\gamma ) \\ {}&=\bigg | \sum _{i=1}^l \int _{t_{i-1}}^{t_i}(t_l-r)^{-\gamma }E_{\alpha ,\beta ,1-\gamma }^{-\delta }(w_1(t_l-r)^\alpha ,w_2(t_l-r)^\beta )(g(r)-P_1(r))'dr\bigg | \\ {}&=\bigg | \sum _{i=1}^l \int _{t_{i-1}}^{t_i}(t_l-r)^{-1-\gamma }E_{\alpha ,\beta ,-\gamma }^{-\delta }(w_1(t_l-r)^\alpha ,w_2(t_l-r)^\beta )(g(r)-P_2(r))dr\bigg | \\ {}&=\bigg | \sum _{i=1}^l \frac{g''(\eta _i)}{2}\int _{t_{i-1}}^{t_i}(r-t_l)(r-t_{l-1})(t_l-r)^{-1-\gamma }E_{\alpha ,\beta , -\gamma }^{-\delta }(w_1(t_l-r)^\alpha ,w_2(t_l-r)^\beta )\bigg | \\ {}&\le \bigg | \sum _{i=1}^l \frac{\zeta ^2g''(\eta _i)}{8}\big [(t_l-t_i)^{-\gamma }E_{\alpha ,\beta ,1-\gamma }^{-\delta }(w_1(t_l-t_i)^\alpha ,w_2(t_l-t_i)^\beta )) \\ {}&\quad -(t_l-t_{i-1})^{-\gamma }E_{\alpha ,\beta ,1-\gamma }^{-\delta }(w_1(t_l-t_{i-1})^\alpha ,w_2(t_l-t_{i-1})^\beta ))\big ]\bigg | \\ {}&\le \bigg | \sum _{i=1}^l \frac{\zeta ^{2-\gamma }g''(\eta _i)}{8}\big [(l-i)^{-\gamma }E_{\alpha ,\beta ,1-\gamma }^{-\delta }(w_1(t_l-t_i)^\alpha ,w_2(t_l-t_i)^\beta )) \\ {}&\quad -(l-i+1)^{-\gamma }E_{\alpha ,\beta ,1-\gamma }^{-\delta }(w_1(t_l-t_{i-1})^\alpha ,w_2(t_l-t_{i-1})^\beta ))\big ]\bigg | \\ {}&\le \frac{ \zeta ^{2-\gamma }}{8} \underset{0 \le t \le t_l}{\max }\big |g''(t)\big | K_1. \end{aligned} \end{aligned}$$

$\square$

For the numerical computations we used Matlab and the ML function $E_{\alpha ,\beta ,\gamma }^{\delta }$ as defined in (1) is approximated by truncating the series with 50 terms. We define the absolute error function

$$\begin{aligned} \big (Er_{\alpha ,\beta ,\gamma ,w_{1},w_{2}}^{\delta ,(n,N_{1})}f\big ) (t)=\left| \big (_{0}^{C}D_{\alpha ,\beta ,\gamma ,w_{1},w_{2}}^{\delta }f \big )(t)-\big (_{0}^{C}D_{\alpha ,\beta ,\gamma ,w_{1},w_{2}}^{\delta ,(n,N_{1})}f\big )(t)\right| . \end{aligned}$$

Example 4

We consider the functions $f\left( t\right) =t^{\mu }$ $\mu \ge 1$ for $t\in \Omega =\left[ 0,1\right]$ and the approximation of Caputo type derivative $\big (_{0}^{C}D_{\alpha ,\beta ,\gamma ,w_{1},w_{2}}^{\delta ,\left( 1,N_{1}\right) }f\big )(t)$ by using linear Lagrange interpolation function $P_1(t)$ as defined in (21) and the exact derivative $\big (_{0}^{C}D_{\alpha ,\beta ,\gamma ,w_{1},w_{2}}^{\delta }f\big )(t)$ are considered. We take the set of mesh points $\Omega ^{\zeta }$ as given in (8), also, $w_{1}=w_{2}=1$ and $\alpha =0.3,\beta =0.9$ and $\gamma =0.6$. Figure 3, presents $\big (_{0}^{C}D_{0.3,0.9,0.6,1,1}^{-0.5,\left( 1,64\right) }f\big )(t)$ and the exact derivative $\big ( _{0}^{C}D_{0.3,0.9,0.6,1,1}^{-0.5}f\big )(t)$ denoted by ap and ex respectively for $\mu =1$. Further, Fig. 4, presents $\big ( _{0}^{C}D_{0.3,0.9,0.6,1,1}^{-0.5,\left( 1,64\right) }f\big )(t)$ and the exact derivative $\big (_{0}^{C}D_{0.3,0.9,0.6,1,1}^{-0.5}f\big )(t)$ for $\mu =4$. It can be viewed from Fig. 3 that the approximation is almost exact in double precision for $f\left( t\right) =t.$ Furthermore, Fig. 5 shows the error function $\big ( Er_{0.3,0.9,0.6,1,1}^{-0.5,(1,64)}f\big )(t)$ for $f\left( t\right) =t^{4}.$

Fig. 3

$\big ( _{0}^{C}D_{0.3,0.9,0.6,1,1}^{-0.5,\left( 1,64\right) }f\big )(t)$ and the exact derivative $\big (_{0}^{C}D_{0.3,0.9,0.6,1,1}^{-0.5}f\big )(t)$ denoted by ap and ex respectively for $f\left( t\right) =t.$

Fig. 4

$\big ( _{0}^{C}D_{0.3,0.9,0.6,1,1}^{-0.5,\left( 1,64\right) }f\big )(t)$ and the exact derivative $\big (_{0}^{C}D_{0.3,0.9,0.6,1,1}^{-0.5}f\big )(t)$ denoted by ap and ex respectively for $f\left( t\right) =t^{4}$

Fig. 5

The error function $\big ( Er_{0.3,0.9,0.6,1,1}^{-0.5,(1,64)}f\big )(t)$ for $f\left( t\right) =t^{4}$

Now we construct the quadratic Lagrange interpolation polynomial for the nodes $t_{i-2}$, $t_{i-1}$ and $t_{i}$.

$$\begin{aligned} \begin{aligned} P_2(t) {}&=\frac{(t-t_{i-1})(t-t_{i})(t_i-t_{i-1})}{(t_{i-1}-t_{i-2})(t_i-t_{i-1})(t_i-t_{i-2})}g(t_{i-2})\\&-\frac{(t-t_{i-1})(t-t_{i})(t_{i}-t_{i-2})}{(t_{i-1}-t_{i-2})(t_i-t_{i-1})(t_i-t_{i-2})}g(t_{i-1}) \\ {}&\quad +\frac{(t-t_{i-2})(t-t_{i-1})(t_{i-1}-t_{i-2})}{(t_{i-1}-t_{i-2})(t_i-t_{i-1})(t_i-t_{i-2})}g(t_{i}). \end{aligned} \end{aligned}$$

(26)

Further,

$$\begin{aligned} \frac{d^2}{dt^2}P_2(t) =\frac{2g(t_{i-2})(t_i-t_{i-1})-2g(t_{i-1})(t_{i}-t_{i-2})+2g(t_{i}) (t_{i-1}-t_{i-2})}{(t_{i-1}-t_{i-2})(t_i-t_{i-1})(t_i-t_{i-2})}, \end{aligned}$$

and

$$\begin{aligned} g(t)-P_2(t)=\frac{f'''(\epsilon _i)}{6}(t-t_i)(t-t_{i-1})(t-t_{i-2}), \end{aligned}$$

where $\epsilon _i\in (t_{i-2},t_i)$.

For $1<\gamma <2$ and $0<t<T$ by choosing $n=2$, the Caputo type derivative (5) is defined as

$$\begin{aligned} {}^{{C}}_{{0}}{D}_{\alpha ,\beta ,\gamma }^{\delta ,w_1,w_2}g(t)=\int _0^t (t-r)^{1-\gamma } E_{\alpha ,\beta ,2- \gamma }^{-\delta }(w_1(t-r)^\alpha ,w_2(t-r)^\beta )\frac{d^2}{dr^2}g(r)dr, \end{aligned}$$

at points $t=t_l=l\frac{T}{N_1}$, $l=1,\dots ,N_1$,

$$\begin{aligned} \begin{aligned} {}&\big ({}^{{C}}_{{0}}{D}_{\alpha ,\beta , \gamma ,w_1,w_2}^{\delta ,(2,N_1)}g\big )(t_l) \\ {}&\quad \approx \int _0^{t_l}(t_l-r)^{1-\gamma }E_{\alpha ,\beta ,2-\gamma }^{-\delta }(w_1(t_l-r)^\alpha ,w_2(t_l-r)^\beta )P_2''(r)dr \\ {}&\quad =\sum _{i=1}^l \int _{t_{i-1}}^{t_i}(t_l-r)^{1-\gamma }E_{\alpha ,\beta ,2-\gamma }^{-\delta }(w_1(t_l-r)^\alpha ,w_2(t_l-r)^\beta )P_2''(r)dr \\ {}&\quad =\int _{t_0}^{t_1}(t_l-r)^{1-\gamma }E_{\alpha ,\beta ,2-\gamma }^{-\delta }(w_1(t_l-r)^\alpha ,w_2(t_l-r)^\beta )\\ {}&\qquad \times \bigg (\frac{g(t_0)(t_2-t_1)-g(t_1)(t_2-t_0)+g(t_2)(t_1-t_0)}{\zeta ^3}\bigg )dr \\ {}&\qquad +\sum _{i=2}^l \int _{t_{i-1}}^{t_i}(t_l-r)^{1-\gamma }E_{\alpha ,\beta ,2-\gamma }^{-\delta }(w_1(t_l-r)^\alpha ,w_2(t_l-r)^\beta ) \\ {}&\qquad \times \bigg (\frac{g(t_{i-2})(t_i-t_{i-1})-g(t_{i-1})(t_{i}-t_{i-2})+g(t_{i})(t_{i-1}-t_{i-2})}{\zeta ^3}\bigg )dr \\ {}&\quad =\zeta ^{-\gamma }\bigg [\bar{w}_{1,l-1}g(t_0)-2\bar{w}_{1,l-1}g(t_1)+\bar{w}_{1,l-1}g(t_2) \\ {}&\qquad + \sum _{i=2}^l \bar{w}_{1,l-i}g(t_{i-2})- 2\bar{w}_{1,l-i}g(t_{i-1})+\bar{w}_{1,l-i}g(t_i)\bigg ], \end{aligned} \end{aligned}$$

where

$$\begin{aligned} \begin{aligned} \bar{w}_{1,l-i}{}&=(l-i+1)^{2-\gamma }E_{\alpha ,\beta ,3-\gamma }^{-\delta }(w_1(t_l-t_{i-1})^\alpha ,w_2(t_l-t_{i-1})^\beta )\\&\quad -(l-i)^{2-\gamma }E_{\alpha ,\beta ,3-\gamma }^{-\delta }(w_1(t_l-t_{i})^\alpha ,w_2(t_l-t_{i})^\beta ). \end{aligned} \end{aligned}$$

In the next theorem, we determine error bound for the given approximation.

Theorem 10

Let $g(t) \in C^3[0,T]$ for any $1<\gamma <2$, the truncation error $R_2(g,\zeta ,\gamma )$ is given as follows

$$\begin{aligned} R_2(g,\zeta ,\gamma ) \le \frac{ \zeta ^{3-\gamma }}{48} \underset{0 \le t \le t_l}{\max }\big |g'''(t)\big | K_2, \end{aligned}$$

where $K_2$ is a positive constant

$$\begin{aligned} \begin{aligned} K _2{}&=\bigg | (l-i+1)^{-\gamma }E_{\alpha ,\beta ,1-\gamma }^{-\delta }(w_1(t_l-t_{i-1})^\alpha ,w_2(t_l-t_{i-1})^\beta )\\&\quad -(l-i)^{-\gamma }E_{\alpha ,\beta ,1-\gamma }^{-\delta }(w_1(t_l-t_{i})^\alpha ,w_2(t_l-t_{i})^\beta )\bigg |. \end{aligned} \end{aligned}$$

Proof

$$\begin{aligned} \begin{aligned}&R_2(g,\zeta ,\gamma ) {}\\&=\bigg | \sum _{i=1}^l \int _{t_{i-1}}^{t_i}(t_l-r)^{1-\gamma }E_{\alpha ,\beta ,2-\gamma }^{-\delta }(w_1(t_l-r)^\alpha ,w_2(t_l-r)^\beta )(g(r)-P_2(r))''dr\bigg |\\&=\bigg | \sum _{i=1}^l \int _{t_{i-1}}^{t_i}(t_l-r)^{-1-\gamma }E_{\alpha ,\beta ,-\gamma }^{-\delta }(w_1(t_l-r)^\alpha ,w_2(t_l-r)^\beta )(g(r)-P_2(r))dr\bigg | \\ {}&=\bigg | \sum _{i=1}^l \frac{g'''(\eta _i)}{6}\int _{t_{i-1}}^{t_i}(r-t_l)(r-t_{l-1})(r-t_{l-2})(t_l-r)^{-1-\gamma } \\ {}&\quad \times E_{\alpha ,\beta ,-\gamma }^{-\delta }(w_1(t_l-r)^\alpha ,w_2(t_l-r)^\beta )\bigg | \\ {}&\le \bigg | \sum _{i=1}^l \frac{\zeta ^3g''(\eta _i)}{48}\big [(t_l-t_i)^{-\gamma }E_{\alpha ,\beta ,1-\gamma }^{-\delta }(w_1(t_l-t_i)^\alpha ,w_2(t_l-t_i)^\beta )) \\ {}&\quad -(t_l-t_{i-1})^{-\gamma }E_{\alpha ,\beta ,1-\gamma }^{-\delta }(w_1(t_l-t_{i-1})^\alpha ,w_2(t_l-t_{i-1})^\beta ))\big ]\bigg | \\ {}&\le \bigg | \sum _{i=1}^l \frac{\zeta ^{3-\gamma }g'''(\eta _i)}{48}\big [(l-i)^{-\gamma }E_{\alpha ,\beta ,1-\gamma }^{-\delta } (w_1(t_l-t_i)^\alpha ,w_2(t_l-t_i)^\beta )) \\ {}&\quad -(l-i+1)^{-\gamma }E_{\alpha ,\beta ,1-\gamma }^{-\delta }(w_1(t_l-t_{i-1})^\alpha ,w_2(t_l-t_{i-1})^\beta ))\big ]\bigg | \\ {}&\le \frac{ \zeta ^{3-\gamma }}{48} \underset{0 \le t \le t_l}{\max }\big |g'''(t)\big |K_2. \end{aligned} \end{aligned}$$

$\square$

Example 5

We take the function $f\left( t\right) =t^{4 }$ for $t\in \Omega =\left[ 0,1\right]$ and the approximation of Caputo type derivative $\big (_{0}^{C}D_{\alpha ,\beta ,\gamma ,w_{1},w_{2}}^{\delta ,\left( 2,N_{1}\right) }f\big )(t)$ by using quadratic Lagrange interpolation polynomial $P_2(t)$ given in (26) and the exact derivative $\big (_{0}^{C}D_{\alpha ,\beta ,\gamma ,w_{1},w_{2}}^{\delta }f\big )(t)$ are considered. We take the set of mesh points $\Omega ^{\zeta }$ as given in (8), also, $w_{1}=w_{2}=1$ and $\alpha =0.3,\beta =0.9$ and $\gamma =1.7$. Figure 6 shows the error function $\big ( Er_{0.3,0.9,1.7,1,1}^{-0.5,(2,128)}f\big )(t)$ for $f\left( t\right) =t^{4}$.

Fig. 6

The error function $\big ( Er_{0.3,0.9,1.7,1,1}^{-0.5,(2,128)}f\big )(t)$ for $f\left( t\right) =t^{4}$

5 Considering the operators with respect to functions

The concept of applying fractional integration and differentiation operators to function $\psi (x)$ with respect to a monotonic function h(x) was studied [32, 33] and Osler’s 1970 paper marked the introduction of the complete generality of Riemann–Liouville fractional calculus with respect to functions. Within this section, we aim to generalize the operators with bivariate ML kernel with respect to an arbitrary monotonic function.

Definition 1

Let $\psi \in L^1[a,b]$ and $h \in C^1[a,b]$ be two functions defined on a real interval [a, b] additionally let h be positive and monotonically increasing. For $\alpha , \beta ,\gamma , \delta , w_1, w_2 \in \mathbb {C}, Re(\beta )>0, Re(\alpha )>0, Re(\gamma )>0$ the fractional integral operator with bivariate ML kernel of the function $\psi$ with respect to the function h is defined by

$$\begin{aligned} \begin{aligned} {}&{}_{{a}}{I}_{\alpha ,\beta ,\gamma ;h(x)}^{\delta ,w_1,w_2}\psi (x) \\ {}&=\int _a^x \big (h(x)-h(t)\big )^{\gamma -1} E_{\alpha ,\beta , \gamma }^{\delta }\bigg (w_1(h(x)-h(t))^\alpha ,w_2(h(x)-h(t))^\beta \bigg )\psi (t)h'(t)dt. \end{aligned} \end{aligned}$$

(27)

Definition 2

Let $h \in C^1[a,b]$ be a positive and monotonically increasing function on a real interval [a, b]. For $\alpha , \beta ,\gamma , \delta , w_1, w_2 \in \mathbb {C}, Re(\beta )>0, Re(\alpha )>0, Re(\gamma )>0$ and let n be a natural number satisfying the condition $n-1<Re(\gamma )<n$. The fractional derivative operator with bivariate ML kernel of the function $\psi \in C^n[a,b]$ with respect to the function h is defined by

$$\begin{aligned} {}_{{a}}{D}_{\alpha ,\beta ,\gamma ;h(x)}^{\delta ,w_1,w_2} \psi (x)=\bigg (\frac{1}{h'(x)}.\frac{d}{dx}\bigg )^n\bigg ({}_{{a}} {I}_{\alpha ,\beta ,n-\gamma ;h(x)}^{-\delta ,w_1,w_2}\psi (x)\bigg ). \end{aligned}$$

(28)

Theorem 11

The operators (27) and (28) can be expressed as conjugates of the operators with respect to x

$$\begin{aligned}&{}_{{a}}{I}_{\alpha ,\beta ,\gamma ;h(x)}^{\delta ,w_1,w_2}=Q_h\circ {}_{{h(a)}}{I}_{\alpha ,\beta ,\gamma }^{\delta ,w_1,w_2}\circ Q_h^{-1}, \end{aligned}$$

(29)

$$\begin{aligned}&{}_{{a}}{D}_{\alpha ,\beta ,\gamma ;h(x)}^{\delta ,w_1,w_2}=Q_h\circ {}_{{h(a)}}{D}_{\alpha ,\beta ,\gamma }^{\delta ,w_1,w_2}\circ Q_h^{-1}, \end{aligned}$$

(30)

where the operator $Q_h$ is defined as,

$$\begin{aligned} \big (Q_h(g)\big )(x)=(g\circ h)(x)=g\big (h(x)\big ). \end{aligned}$$

Proof

The proof can be given in a similar manner as the classical Riemann–Liouville fractional calculus with respect to functions. It is well established that the classical derivative satisfies the corresponding operational identity

$$\begin{aligned} {}^{{RL}}D_{h(x)}^1=\frac{1}{h'(x)}\frac{d}{dx}=Q_h\circ {}^{{RL}}D^1 \circ Q_h^{-1}. \end{aligned}$$

(31)

We first select a function $g \in L^1[a,b]$ and then follow the subsequent steps:

$$\begin{aligned}&g:x\mapsto g(x),\\&Q_h^{-1} g: x \mapsto g\big (h^{-1}(x)\big ),\\&{}_{{h(a)}}{I}_{\alpha ,\beta ,\gamma }^{\delta ,w_1,w_2}\circ Q_h^{-1}g : x\mapsto {}_{{h(a)}}{I}_{\alpha ,\beta ,\gamma }^{\delta ,w_1,w_2} \big (g \circ h^{-1}\big ) (x),\\&Q_h\circ {}_{{h(a)}}{I}_{\alpha ,\beta ,\gamma }^{\delta ,w_1,w_2}\circ Q_h^{-1}g : x\mapsto \big ( {}_{{h(a)}}{I}_{\alpha ,\beta ,\gamma }^{\delta ,w_1,w_2} \big (g \circ h^{-1}\big )\big ) (h(x)). \end{aligned}$$

Based on definition (2), we have

$$\begin{aligned} \begin{aligned} {}&{}_{{h(a)}}{I}_{\alpha ,\beta ,\gamma }^{\delta ,w_1,w_2} \big (g \circ h^{-1}\big )(x) \\ {}&= \int _{h(a)}^x (x-s)^{\gamma -1} E_{\alpha ,\beta , \gamma }^{\delta }(w_1(x-s)^\alpha ,w_2(x-s)^\beta )g(h^{-1}(s))ds \\ {}&= \int _{a}^{h^{-1}(x)}(x-h(v))E_{\alpha ,\beta , \gamma }^{\delta }(w_1(x-h(v))^\alpha ,w_2(x-h(v))^\beta )g(v)h'(v)dv, \end{aligned} \end{aligned}$$

where we perform the substitution of $v=h^{-1}(s)$ into the integral. Finally, by substituting x with h(x) throughout this expression, we obtain

$$\begin{aligned} \begin{aligned} {}&Q_h\circ {}_{{h(a)}}{I}_{\alpha ,\beta ,\gamma }^{\delta ,w_1,w_2}\circ Q_h^{-1}g(x) \\ {}&= \int _{a}^x (h(x)-h(v))E_{\alpha ,\beta , \gamma }^{\delta }(w_1(h(x)-h(v))^\alpha ,w_2(h(x)-h(v))^\beta )g(v)h'(v)dv \\ {}&={}_{{h(a)}}{I}_{\alpha ,\beta ,\gamma }^{\delta ,w_1,w_2}g(x), \end{aligned} \end{aligned}$$

this establishes (29). The result in (30) for fractional derivatives directly follows from the composition of (29) with (31) repeated n times. $\square$

Theorem 12

The following relation holds true for the operators (27)

$$\begin{aligned} {}^{{RL}}_{{h(x)}}I_{a}^{\mu }\bigg ({}_{{a}} {I}_{\alpha ,\beta ,\gamma ;h(x)}^{\delta ,w_1,w_2}\psi (x)\bigg )={}_{{a}} {I}_{\alpha ,\beta ,\gamma +\mu ;h(x)}^{\delta ,w_1,w_2}={}_{{a}} {I}_{\alpha ,\beta ,\gamma ;h(x)}^{\delta ,w_1,w_2}\bigg ({}^{{RL}}_{{h(x)}}I_{a }^{\mu }\psi (x)\bigg ), \end{aligned}$$

for any function $\psi \in L^1[a,b]$ and for any monotonic function $h\in C^1[a,b]$ where $\alpha , \beta ,\gamma , \delta , w_1, w_2 \in \mathbb {C}, Re(\beta )>0, Re(\alpha )>0, Re(\gamma )>0, Re(\mu )>0.$

Proof

This is a result of Corollary 4 [13], utilizing the conjugation expressions for all operators. These expressions are provided by Theorem 11 for operators with bivariate ML kernels and by the classical result for Riemann–Liouville integrals

$$\begin{aligned} {}^{{RL}}_{{h(x)}}I_{a}^{\mu }=Q_h\circ {}^{{RL}}I_{h(a)}^{\mu } \circ Q_h^{-1}. \end{aligned}$$

$\square$

Theorem 13

The following relation holds true for the operators (27)

$$\begin{aligned} {}^{{RL}}_{{h(x)}}D_{a}^{\mu } \bigg ({}_{{a}}{I}_{\alpha ,\beta ,\gamma ;h(x)}^{\delta ,w_1,w_2} \psi (x)\bigg )=\bigg ({}_{{a}}{I}_{\alpha ,\beta ,\gamma -\mu ;h(x)}^{\delta ,w_1,w_2}\psi (x)\bigg ), \end{aligned}$$

for any function $\psi \in C^k[a,b]$, $k=\lceil \mu \rceil$, any monotonic function $h \in C^1[a,b]$ and $\alpha , \beta ,\gamma , \delta , w_1, w_2 \in \mathbb {C}, Re(\beta )>0, Re(\alpha )>0, Re(\gamma )>0, Re(\mu )>0.$

Proof

This is a result of Corollary 4 [13], using results of Theorem 11 and

$$\begin{aligned} {}^{{RL}}_{{h(x)}}D_{a}^{\mu }=Q_h\circ {}^{{RL}}D_{h(a)}^{\mu } \circ Q_h^{-1}. \end{aligned}$$

$\square$

Theorem 14

The operator (27) satisfies the following semigroup property given by

$$\begin{aligned} {}_{{a}}{I}_{\alpha ,\beta ,\gamma _1;h(x)}^{\delta _1,w_1,w_2} \bigg ({}_{{a}}{I}_{\alpha ,\beta ,\gamma _2;h(x)}^{\delta _2,w_1,w_2} \psi (x)\bigg )={}_{{a}}{I}_{\alpha ,\beta ,\gamma _1+\gamma _2;h(x)}^{\delta _1+\delta _2,w_1,w_2}, \end{aligned}$$

for any function $\psi \in L^1[a,b]$ and for any monotonic function $h \in C^1[a,b]$, where $\alpha , \beta ,\gamma _1, \gamma _2, \delta _1,\delta _2,w_1, w_2 \in \mathbb {C}$, $Re(\beta )>0, Re(\alpha )>0, Re(\gamma _1)>0, Re(\gamma _2)>0.$

Proof

This is a result of Theorem 9 [13], in combination with results of Theorem 11. $\square$

Theorem 15

The application of the bivariate ML integral operator to another power function can be demonstrated in the following manner

$$\begin{aligned} {}_{{a}}{I}_{\alpha ,\beta ,\gamma ;(x-a)^\mu }^{\delta ,w_1,w_2}(x-a)^\zeta =\Gamma (\frac{\zeta }{\mu }+1)(x-a)^{\gamma \mu +\zeta }E_{\alpha ,\beta ,\gamma +\frac{\zeta }{\mu }}^\delta (w_1(x-a)^{\alpha \mu },w_2(x-a)^{\beta \mu }), \end{aligned}$$

(32)

where $\alpha , \beta ,\gamma , \delta , w_1, w_2 \in \mathbb {C}, Re(\beta )>0, Re(\alpha )>0, Re(\gamma )>0, \mu>0, \zeta >0$.

Proof

From the assumption $\mu >0$ the function $h(x)=(x-a)^\mu$ is a monotonically increasing function on any interval [a, c]. Further, the function to which we apply the operator (27) is $(x-a)^\zeta =h(x)^{\frac{\zeta }{\mu }}$. Thus, utilizing the result of Proposition 1 we have

$$\begin{aligned}&\psi : x \mapsto (x-a)^\zeta ,\\&Q_h^{-1}\psi : x \mapsto x^{\frac{\zeta }{\mu }},\\&{}_{{0}}{I}_{\alpha ,\beta ,\gamma }^{\delta ,w_1,w_2}\circ Q_h^{-1}\psi : x \mapsto \big ({}_{{0}}{I}_{\alpha ,\beta ,\gamma }^{\delta ,w_1,w_2}\big ) x^{\frac{\zeta }{\mu }} =\Gamma (\frac{\zeta }{\mu }+1)x^{\gamma +\frac{\zeta }{\mu }}E_{\alpha ,\beta ,\gamma +\frac{\zeta }{\mu }}^\delta (w_1x^{\alpha },w_2x^{\beta }),\\&\begin{aligned} {}&Q_h\circ {}_{{0}}{I}_{\alpha ,\beta ,\gamma }^{\delta ,w_1,w_2}\circ Q_h^{-1}\psi : x\mapsto \big ( {}_{{0}}{I}_{\alpha ,\beta ,\gamma }^{\delta ,w_1,w_2}\big )(x-a)^\zeta \\ {}&=\Gamma (\frac{\zeta }{\mu }+1)(x-a)^{\gamma \mu +\zeta }E_{\alpha ,\beta ,\gamma +\frac{\zeta }{\mu }}^\delta (w_1(x-a)^{\alpha \mu },w_2(x-a)^{\beta \mu }), \end{aligned} \end{aligned}$$

thus (32) is established. $\square$

Theorem 16

The application of the bivariate ML integral operator with respect to a logarithm function to a power function can be given as

$$\begin{aligned} {}_{{a}}{I}_{\alpha ,\beta ,\gamma ;log(x-a)}^{\delta ,w_1,w_2} (x-a)^\zeta =\zeta ^{-\gamma }\bigg (1-\frac{w_1}{\zeta ^\alpha }-\frac{w_2}{\zeta ^\beta }\bigg )^{-\delta }(x-a)^\zeta , \end{aligned}$$

(33)

where $\alpha , \beta ,\gamma , \delta , w_1, w_2 \in \mathbb {C}, Re(\beta )>0, Re(\alpha )>0, Re(\gamma )>0, \zeta >0$.

Proof

The function $h(x)=log(x-a)$ is monotonically increasing function on any interval (a, c] with $h(x)\rightarrow -\infty$ as $x\rightarrow a^+$, and we apply the operator (27) to $(x-a)^\zeta =e^{\zeta h(x)}$. Using the Proposition 2 it follows that

$$\begin{aligned}&\psi : x \mapsto (x-a)^\zeta ,\\&Q_h^{-1}\psi : x \mapsto e^{\zeta x},\\&{}_{{-\infty }}{I}_{\alpha ,\beta ,\gamma }^{\delta ,w_1,w_2}\circ Q_h^{-1}\psi : x \mapsto \big ({}_{{0}}{I}_{\alpha ,\beta ,\gamma }^{\delta ,w_1,w_2}\big ) e^{\zeta x} =\zeta ^{-\gamma }\bigg (1-\frac{w_1}{\zeta ^\alpha }-\frac{w_2}{\zeta ^\beta }\bigg )^{-\delta } e^{\zeta x},\\&Q_h\circ {}_{{0}}{I}_{\alpha ,\beta ,\gamma }^{\delta ,w_1,w_2}\circ Q_h^{-1}\psi : x\mapsto \big ( {}_{{0}}{I}_{\alpha ,\beta ,\gamma }^{\delta ,w_1,w_2}\big )(x-a)^\zeta =\zeta ^{-\gamma }\bigg (1-\frac{w_1}{\zeta ^\alpha }-\frac{w_2}{\zeta ^\beta }\bigg )^{-\delta }(x-a)^\zeta , \end{aligned}$$

hence we get (33). $\square$

Example 6

In this example we consider the integral operators given in ( 32) and (33) over the interval $\Omega =\left[ 0,1\right]$ at the mesh points (8) for $N_{1}=64$. Taking $\alpha =0.3,\beta =0.9, \gamma =0.6$ and $\delta =0.7$ the integral $_{0}I_{0.3,0.9,0.6;x^ {\mu }}^{0.7,1,1}\left[ x^{\zeta }\right]$ is presented in Fig. 7 for various values of $\zeta$ as 0.5, 0.9, 1.9 when $\mu =1.2$ is fixed and for various values of $\mu$ as 0.5, 0.9, 1.9 when $\zeta =0.8$ is fixed. Further, Taking $\alpha =0.3,\beta =0.9$ and $\delta =0.7$ the integral $_{0}I_{0.3,0.9,0.8;log(x)}^{0.7,1,1}\left[ x^{\zeta }\right]$ is presented in Fig. 8 for various values of $\zeta$ as 0.5, 0.9, 1.9 when $\gamma =0.8$ is fixed and for various values of $\gamma$ as 0.5, 0.9, 1.9 for $\zeta =0.8$.

Fig. 7

The integral $_{0}I_{\alpha ,\beta ,\gamma ;x^ {\mu }}^{\delta ,w_{1},w_{2}}\left[ x^{\zeta }\right]$ for various values of $\zeta$ and $\mu$

Fig. 8

The integral $_{0}I_{\alpha ,\beta ,\gamma ;log(x)}^{\delta ,w_{1},w_{2}}\left[ x^{\zeta }\right]$ for various values of $\gamma$ and $\zeta$

6 Conclusion

In this paper, we have extended the work of [13] by developing a fully formed theory of fractional calculus from integral operators with bivariate ML kernels. Our approach naturally leads to many important results concerning these integral operators, such as the Laplace transform, product rule, and chain rule. We have also extended this model of fractional calculus to a higher level of generality by applying the operators with respect to functions as well as with respect to the independent variable. The composition, semigroup, and inverse properties extend naturally into this more general setting.

As applications of the work in this paper, we considered the solution of an eigenvalue problem involving Caputo type derivative operators with bivariate ML kernels, as well as a Cauchy problem involving both Riemann–Liouville derivatives and the fractional integrals with bivariate ML kernels. We discussed linear and quadratic approximations to approximate the Caputo type derivative operators with bivariate ML kernels.

As a future work we are intending to give numerical approximations for the bivariate and multivariate franctional calculus operators (see [24, 34]) by considering different approaches. Another direction of research will be to give new approximation methods for the solution of fractional differential equations.