Spectral Approach to the Fractional Power of Operator and Its Matrix Approximation

1. Introduction

One of the most important problem in the fractional calculus and its applications is to derive a formula that, for suitable operators A and a range of parameters $\alpha$, generates a family of operators $\left\{A^{\alpha }\right\}$ possessing the properties expected of powers. In particular, $A^n$ should coincide with the iterated power $A^n={\underbrace{A·A·A·...·A}}_n$ when $\alpha$ is a positive integer n, and the index law $A^{\alpha }{A}^{\beta }=A^{\alpha +\beta }$ should hold whenever $A^{\alpha }$, $A\beta$, and $A^{\alpha +\beta }$ are all defined. Several methods now exist for constructing such families [1]. Let us consider some methods of constructing powers of operator A.

The first method involves finding the iterated operator $A^n$ and replacing the natural number n with a real or complex number $\alpha$. This approach works when the functions of natural numbers which are presented in the formula for the iterated operator can be generalized through their extensions to real (complex) numbers. For example, the fractional power of the regular integral arises in this way. Namely, starting from the formula for the iterated integral, we see that it is sufficient to replace the power function with a natural exponent by a power function with a real exponent, and the factorial by the Gamma function. This yields the well-known Riemann-Liouville fractional integral (see [2], p. 33). By finding the inverse operator, we obtain the Riemann-Liouville fractional derivative ([2], p. 30).

Another approach is to begin with the definition of the derivative as the limit of the ratio of the function’s increment to the argument’s increment as the argument’s increment tends to zero, derive a finite-difference formula for the derivative of a positive integer order, and generalize it to a real order. This results in the Grünwald–Letnikov fractional derivative, whose inverse is the corresponding fractional integral in a difference form (see [3], p. 43). However, if the operator is more general than an integral or a derivative, then using the iterative method becomes technically quite difficult.

A further method is based on solving the Cauchy problem for the iterated operator and then replacing the positive integer denoting the order of the operator with a real number (see [4]). Additionally, if the operator can be presented in a factorized form using integral transforms and a symbol, its fractional power can be expressed as a pseudo-differential operator with the symbol raised to a real power (see, for example, the fractional Laplace operator in [2], p. 483).

All listed approaches are not constructive. For this reason, they are often not applicable from an applied point of view. Methods for constructing fractional powers of operators based on theirs properties rather than on replacing a natural number with a real one appeared in functional analysis. When the operator $(-A)$ is the infinitesimal generator of a bounded semigroup, Phillips and Hille (see [5], 1952 and [6], first published in 1957) established that its fractional powers could be constructed using an operational calculus they developed. Balakrishnan [7] in 1959 provided a comprehensive treatment of this approach. Subsequently, Balakrishnan [8] in 1960 introduced a new definition based on resolvents that extended the theory to a broader class of operators. Further study of this approach have been given in [9,10,11,12]. Independently, Krasnosel’skii and Sobolevskii [13] in 1959 (see also [14]) using the splitting method, constructed fractional powers of self-adjoint operators.

The significant role of fractional calculus in applications [15,16,17] drives the need to construct approximations of fractional powers of operators and to quantify their convergence rates. The most useful way to represent, disnamical process, and study linear mappings between finite-dimensional vector spaces are matrices. Also, matrix approximation of operators is crucial for simplifying integrals, differential operators and stochastic processes.

In this paper, we present the convergence estimate for fractional powers of operators. Specifically, we begin with a base operator and its approximation, both satisfying the same resolvent estimate. Using the resolvent approach, we then construct the corresponding fractional powers and derive an estimate for the norm of their difference. As a practical application, we also develop a matrix approximation method.

2. The Spectral Approach to the Fractional Powers of Operators

This section is devoted to the study of the spectral approach to the fractional powers of operators, including rate of convergence and examples.

2.1. Application of the Resolvent to the Construction of a Fractional Power

In [18] the rate of convergence for powers of operators constructed from semigroups theory were obtained. In this article, we use a method for defining fractional powers that applies even when the base operator does not generate a semigroup or have a dense domain. This construction was first given by Balakrishnan in [7]. Namely, if the resolvent of the operator A exists for $\lambda >0$, and is $O(1/\lambda )$ for all $\lambda$, then fractional powers was constructed, using an abstract version of the Stieltjes transform.

Recall the definition of ${(-A)}^{\alpha }$ from [7].

Definition 1.

If A is a closed linear operator with domain $D\left(A\right)$ and range in a Banach space X, for which $\left|\right|\lambda R(\lambda ,A)\left|\right|<M<\infty$, $\lambda >0$, then

$${(-A)}^{\alpha }=-{\displaystyle \frac{sin\left(\alpha \pi \right)}{\pi }}\underset{0}{\overset{\infty }{\int }}{\lambda }^{\alpha -1}R(\lambda ;A)Axd\lambda ,x\in D\left(A\right),0<\alpha <1.$$

(1)

For a such that $n-1<\alpha \le n$, $n\in \mathbb{N}$

$${(-A)}^{\alpha }={(-A)}^{\alpha -n+1}{(-A)}^{n-1}x,x\in D\left(A^{n+1}\right).$$

Some important properties were obtained in [7].

1.

The operators ${(-A)}^{\alpha }$ can be extended to be closed linear.

2.

For $0<\alpha +\beta <1$, a semigroup property is valid ${(-A)}^{\alpha +\beta }={(-A)}^{\alpha }{(-A)}^{\beta }$, $x\in D\left(A^2\right)$.

2.2. Rate of Convergence: Spectral Approach

The precise problem considered in this subsection is as follows. Let B be a Banach space, and let $A,C$ be linear operators with domain and range in B. Suppose the resolvents of $A,C$ exist for $\lambda >0$, and are $O(1/\lambda )$ for all $\lambda$. Consider a subordinate Banach space $D\subset B$ such that ${\parallel (A-C)f\parallel }_B\lesssim {\parallel f\parallel }_D$. We pursue two goals: first, to estimate the closeness of the resolvents $\left|R\right(\lambda ;A)-R(\lambda ;C\left)\right|$, and second, to estimate the closeness of the fractional powers ${|(-A)}^{\alpha }-{(-C)}^{\alpha }|$. Using methods from [19], we obtain the following theorem, which provides effective upper bounds for these quantities.

Next we will consider $\alpha \in (0,1)$.

Theorem 2.1.

Assume $A,C$ are linear operators with domain and range in a Banach space B with the resolvents $R(\lambda ,A)$, $R(\lambda ,C)$ defined for $\lambda >0$ such that

$${\left|\right|\lambda R(\lambda ,A)\left|\right|}_B{<M<\infty ,\left|\right|\lambda R(\lambda ,C)\left|\right|}_B<M<\infty .$$

Let there exist a Banach space $D\subset B$ with norm ${\parallel .\parallel }_D\ge {\parallel .\parallel }_B$ (can be just $D=B$ for bounded $A,C$) such that $R(\lambda ,C)$ is defined and has similar bounds in D for $\lambda >0$:

$${\left|\right|\lambda R(\lambda ,C)\left|\right|}_D<M^{\prime }<\infty .$$

Assume moreover that

$${\parallel (A-C)f\parallel }_B\le \omega {\parallel f\parallel }_D,\omega >0.$$

Then

$${\parallel (R(\lambda ;A)-R(\lambda ;C))f\parallel }_B\le {\displaystyle \frac{M{M}^{\prime }\omega }{{\lambda }^2}}{\parallel f\parallel }_D,$$

(2)

and

$$\parallel ({(-A)}^{\alpha }-{(-C)}^{\alpha }){f\parallel }_B\le \gamma {\parallel f\parallel }_D{\omega }^{\alpha },$$

(3)

where

$$\gamma =M\left(2+M^{\prime }{\displaystyle \frac{\alpha }{1-\alpha }}\right).$$

Proof. 

From the second resolvent equality

$$R(\lambda ;A)-R(\lambda ;C)=R(\lambda ;C)(A-C)R(\lambda ;A),$$

it follows that

$${\parallel R(\lambda ,A)f-R(\lambda ,C)f\parallel }_B\le {\displaystyle \frac{M\omega }{\lambda }}{\parallel R(\lambda ,C)f\parallel }_D\le {\displaystyle \frac{\omega M{M}^{\prime }}{{\lambda }^2}}{\parallel f\parallel }_D.$$

Next, since $R(\lambda ,A)Af=\lambda R(\lambda ,A)f-f$ and the same holds for B, applying (1) we obtain

$${\parallel (-A)}^{\alpha }f-{(-C)}^{\alpha }{f\parallel }_B\le \alpha \underset{0}{\overset{\infty }{\int }}{\lambda }^{\alpha }{\parallel R(\lambda ;A)f-R(\lambda ;C)f\parallel }_{B}d\lambda .$$

Decomposing the domain of integration in two parts, $[0,\omega ]$ and $[\omega ,\infty )$, and using (2) only for the second integral we get the estimate

$$\begin{array}{c}{\parallel (-A)}^{\alpha }f-{(-C)}^{\alpha }{f\parallel }_B\le 2M\alpha \underset{0}{\overset{\omega }{\int }}{\lambda }^{\alpha -1}{\parallel f\parallel }_{B}d\lambda +\alpha \omega M{M}^{\prime }\underset{\omega }{\overset{\infty }{\int }}{\lambda }^{\alpha -2}{\parallel f\parallel }_{D}d\lambda =\hfill \\ \hfill =2M{\omega }^{\alpha }{\parallel f\parallel }_B+{\displaystyle \frac{\alpha }{1-\alpha }}{M}^{\prime }M{\omega }^{\alpha }{\parallel f\parallel }_D\le M\left(2+M^{\prime }{\displaystyle \frac{\alpha }{1-\alpha }}\right){\parallel f\parallel }_D{\omega }^{\alpha }.\end{array}$$

The proof is complete. □

Remark 1.

In the most important case, when $A,C$ generate Feller semigroups, then the above estimate for the norms of their resolvents in B are known to hold with $M=1$. Moreover, if D is the domain of C with ${\parallel f\parallel }_D={\parallel f\parallel }_B+{\parallel Cf\parallel }_B$ (the canonical choice of the norm in the domain of an operator), then the above estimate for the norm of $R(\lambda ,C)$ in D holds as well with $M^{\prime }=1$, because then

$${\parallel CR(\lambda ,C)f\parallel }_B\le {\displaystyle \frac{1}{\lambda }}{\parallel Cf\parallel }_B\le {\displaystyle \frac{1}{\lambda }}{\parallel f\parallel }_D.$$

In this case, γ in (3) simplifies to $\gamma =(2-\alpha )/(1-\alpha )$. In any case it is seen that these estimates are not good for α close to 1.

Next, we analyze Definition 1 in details and show how to construct ${(-A)}^{\alpha }$ in practice, starting from the operator $(-A)$ and its resolvent.

2.3. Riemann-Liouville Fractional Integral

In this subsection we show how this method works with the simplest integral:

$$\left(I_{a+}^{1}f\right)\left(x\right)=\underset{a}{\overset{x}{\int }}f\left(t\right)dt.$$

It is easy to see, that [20]

$$g=R(\lambda ,-A)f={(A+\lambda I)}^{-1}f={\displaystyle \frac{1}{\lambda }}\left(\sum _{k=0}^{\infty }{\displaystyle \frac{{(-1)}^k{A}^k}{{\lambda }^k}}f\right).$$

(4)

It is known that iterated operator $I_{a+}^1$ has a form

$$\left({\left(I_{a+}^1\right)}^{k}f\right)\left(x\right)=\left(I_{a+}^{k}f\right)\left(x\right)=\underset{k}{\underbrace{\underset{a}{\overset{x}{\int }}dx\underset{a}{\overset{x}{\int }}dx...\underset{a}{\overset{x}{\int }}}}f\left(t\right)dt={\displaystyle \frac{1}{(k-1)!}}\underset{a}{\overset{x}{\int }}{(x-t)}^{k-1}f\left(t\right)dt.$$

Then by (4) the resolvent of the operator $(-I_{a+}^1)$ to be

$$R(\lambda ,-I_{a+}^1)f={\displaystyle \frac{1}{\lambda }}\left(\sum _{k=0}^{\infty }{\displaystyle \frac{{(-1)}^k{A}^k}{{\lambda }^k}}f\right)={\displaystyle \frac{1}{\lambda }}f+{\displaystyle \frac{1}{\lambda }}\left(\sum _{k=1}^{\infty }{\displaystyle \frac{{(-1)}^k}{(k-1)!}}{\displaystyle \frac{1}{{\lambda }^k}}\underset{a}{\overset{x}{\int }}{(x-t)}^{k-1}f\left(t\right)dt\right)=$$

$$={\displaystyle \frac{1}{\lambda }}f+{\displaystyle \frac{1}{{\lambda }^2}}\left(\underset{a}{\overset{x}{\int }}\sum _{k=1}^{\infty }{\displaystyle \frac{{(-1)}^k}{(k-1)!}}{\left({\displaystyle \frac{x-t}{\lambda }}\right)}^{k-1}f\left(t\right)dt\right)=$$

$$={\displaystyle \frac{1}{\lambda }}f-{\displaystyle \frac{1}{{\lambda }^2}}\left(\underset{a}{\overset{x}{\int }}\sum _{k=0}^{\infty }{\displaystyle \frac{{(-1)}^k}{k!}}{\left({\displaystyle \frac{x-t}{\lambda }}\right)}^{k}f\left(t\right)dt\right)={\displaystyle \frac{1}{\lambda }}f-{\displaystyle \frac{1}{{\lambda }^2}}\underset{a}{\overset{x}{\int }}{e}^{-{\displaystyle \frac{x-t}{\lambda }}}f\left(t\right)dt.$$

So we can write

$$\left(R(\lambda ,-I_{a+}^1)I_{a+}^{1}f\right)\left(x\right)={\displaystyle \frac{1}{\lambda }}\underset{a}{\overset{x}{\int }}{e}^{-{\displaystyle \frac{x-t}{\lambda }}}f\left(t\right)dt,$$

then by (1) we get

$$\begin{array}{c}\left({\left(I_{a+}^1\right)}^{\alpha }f\right)\left(x\right)={\displaystyle \frac{sin\left(\alpha \pi \right)}{\pi }}\underset{0}{\overset{\infty }{\int }}{\lambda }^{\alpha -2}\left(\underset{a}{\overset{x}{\int }}{e}^{-{\displaystyle \frac{x-t}{\lambda }}}f\left(t\right)dt\right)d\lambda =\hfill \\ ={\displaystyle \frac{sin\left(\alpha \pi \right)}{\pi }}\underset{a}{\overset{x}{\int }}f\left(t\right)\left(\underset{0}{\overset{\infty }{\int }}{\lambda }^{\alpha -2}{e}^{-{\displaystyle \frac{x-t}{\lambda }}}d\lambda \right)dt=\{1/\lambda =z\}=\\ ={\displaystyle \frac{sin\left(\alpha \pi \right)}{\pi }}\underset{a}{\overset{x}{\int }}f\left(t\right)\left(\underset{0}{\overset{\infty }{\int }}{z}^{-\alpha }{e}^{-(x-t)z}dz\right)dt=\\ \hfill ={\displaystyle \frac{sin\left(\alpha \pi \right)\mathrm{\Gamma }(1-\alpha )}{\pi }}\underset{a}{\overset{x}{\int }}{(x-t)}^{\alpha -1}f\left(t\right)dt={\displaystyle \frac{1}{\mathrm{\Gamma }\left(\alpha \right)}}\underset{a}{\overset{x}{\int }}{(x-t)}^{\alpha -1}f\left(t\right)dt=\left(I_{a+}^{\alpha }f\right)\left(x\right).\end{array}$$

Here the Euler’s reflection formula

$$\mathrm{\Gamma }(1-z)\mathrm{\Gamma }\left(z\right)={\displaystyle \frac{\pi }{sin\pi z}},z\notin \mathbb{Z}$$

was used. So, we obtain

$$\left(I_{a+}^{\alpha }f\right)\left(x\right)={\displaystyle \frac{1}{\mathrm{\Gamma }\left(\alpha \right)}}\underset{a}{\overset{x}{\int }}{(x-t)}^{\alpha -1}f\left(t\right)dt.$$

(5)

The expression (5) consides to the well-known Riemann-Liouville fractional integral $I_{a+}^{\alpha }$ (see [2]).

2.4. Erdélyi-Kober Fractional Integral

Looking at the integral

$$\left(I_{a+;x}^{1}f\right)\left(x\right)=\underset{a}{\overset{x}{\int }}f\left(t\right)tdt$$

first, we easily obtain that iterated operator ${\left(I_{a+;x}^1\right)}^k$ has a form

$$\begin{array}{c}\left({\left(I_{a+;x}^1\right)}^{k}f\right)\left(x\right)=\left(I_{a+;x}^{k}f\right)\left(x\right)=\underset{k}{\underbrace{\underset{a}{\overset{x}{\int }}xdx\underset{a}{\overset{x}{\int }}xdx...\underset{a}{\overset{x}{\int }}}}f\left(t\right)tdt=\hfill \\ \hfill \hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}={\displaystyle \frac{1}{2^{k-1}(k-1)!}}\underset{a}{\overset{x}{\int }}{(x^2-t^2)}^{k-1}tf\left(t\right)dt.\end{array}$$

Then by (4) the resolvent of the operator $(-I_{a+}^1)$ to be

$$R(\lambda ,-I_{a+;x}^1)f={\displaystyle \frac{1}{\lambda }}\left(\sum _{k=0}^{\infty }{\displaystyle \frac{{(-1)}^k{A}^k}{{\lambda }^k}}f\right)=$$

$$={\displaystyle \frac{1}{\lambda }}f+{\displaystyle \frac{1}{\lambda }}\left(\sum _{k=1}^{\infty }{\displaystyle \frac{{(-1)}^k}{2^{k-1}(k-1)!}}{\displaystyle \frac{1}{{\lambda }^k}}\underset{a}{\overset{x}{\int }}{(x^2-t^2)}^{k-1}f\left(t\right)tdt\right)=$$

$$={\displaystyle \frac{1}{\lambda }}f+{\displaystyle \frac{1}{{\lambda }^2}}\left(\underset{a}{\overset{x}{\int }}\sum _{k=1}^{\infty }{\displaystyle \frac{{(-1)}^k}{(k-1)!}}{\left({\displaystyle \frac{x^2-t^2}{2\lambda }}\right)}^{k-1}f\left(t\right)tdt\right)=$$

$$={\displaystyle \frac{1}{\lambda }}f-{\displaystyle \frac{1}{{\lambda }^2}}\left(\underset{a}{\overset{x}{\int }}\sum _{k=0}^{\infty }{\displaystyle \frac{{(-1)}^k}{k!}}{\left({\displaystyle \frac{x^2-t^2}{2\lambda }}\right)}^{k}f\left(t\right)tdt\right)={\displaystyle \frac{1}{\lambda }}f-{\displaystyle \frac{1}{{\lambda }^2}}\underset{a}{\overset{x}{\int }}{e}^{-{\displaystyle \frac{x^2-t^2}{2\lambda }}}f\left(t\right)tdt.$$

Since

$$\left((\lambda ,-I_{a+;x}^1)I_{a+;x}^{1}f\right)\left(x\right)={\displaystyle \frac{1}{\lambda }}\underset{a}{\overset{x}{\int }}{e}^{-{\displaystyle \frac{x^2-t^2}{2\lambda }}}f\left(t\right)tdt,$$

then by (1) we get

$$\begin{array}{c}\left({\left(I_{a+;x}^1\right)}^{\alpha }f\right)\left(x\right)={\displaystyle \frac{sin\left(\alpha \pi \right)}{\pi }}\underset{0}{\overset{\infty }{\int }}{\lambda }^{\alpha -2}\left(\underset{a}{\overset{x}{\int }}{e}^{-{\displaystyle \frac{x^2-t^2}{2\lambda }}}f\left(t\right)tdt\right)d\lambda =\hfill \\ ={\displaystyle \frac{sin\left(\alpha \pi \right)}{\pi }}\underset{a}{\overset{x}{\int }}f\left(t\right)t\left(\underset{0}{\overset{\infty }{\int }}{\lambda }^{\alpha -2}{e}^{-{\displaystyle \frac{x^2-t^2}{2\lambda }}}d\lambda \right)dt=\{1/\lambda =z\}=\\ ={\displaystyle \frac{sin\left(\alpha \pi \right)}{\pi }}\underset{a}{\overset{x}{\int }}f\left(t\right)t\left(\underset{0}{\overset{\infty }{\int }}{z}^{-\alpha }{e}^{-{\displaystyle \frac{1}{2}}(x^2-t^2)z}dz\right)dt=\\ ={\displaystyle \frac{sin\left(\alpha \pi \right)\mathrm{\Gamma }(1-\alpha )}{\pi }}\underset{a}{\overset{x}{\int }}{\left({\displaystyle \frac{x^2-t^2}{2}}\right)}^{\alpha -1}f\left(t\right)tdt=\\ \hfill ={\displaystyle \frac{1}{\mathrm{\Gamma }\left(\alpha \right)}}\underset{a}{\overset{x}{\int }}{\left({\displaystyle \frac{x^2-t^2}{2}}\right)}^{\alpha -1}f\left(t\right)tdt=\left(I_{a+;x}^{\alpha }f\right)\left(x\right).\end{array}$$

Since the Erdelýi-Kober integral is defined by formula (see [2]):

$$\left(I_{a+;2}^{\alpha }f\right)\left(x\right)={\displaystyle \frac{2}{\mathrm{\Gamma }\left(\alpha \right)}}{x}^{-2\left(\alpha \right)}\underset{0}{\overset{x}{\int }}{(x^2-t^2)}^{\alpha -1}tf\left(t\right)dt,$$

(6)

operator $I_{a+;x}^{\alpha }$ coincides with (6) up to multiplier $C{x}^{-2\left(\alpha \right)}$ and

$$\left(I_{a+;x}^{\alpha }f\right)\left(x\right)={\displaystyle \frac{1}{\mathrm{\Gamma }\left(\alpha \right)}}\underset{a}{\overset{x}{\int }}{\left({\displaystyle \frac{x^2-t^2}{2}}\right)}^{\alpha -1}f\left(t\right)tdt.$$

2.5. Hadamard-Type Fractional Integral

We will use the function space [1]

$$F_{p,\mu }=\{\phi \in C^{\infty }(0,\infty ):\left|\right|x^k{D}^k\left(x^{-\mu }\phi \right){\left|\right|}_p<\infty ,k=0,1,2,...\},$$

where $1\le p<\infty$, $\mu$ is any real number, $D={\displaystyle \frac{d}{dx}}$, $\left|\right|·{\left|\right|}_p$ is $L_p$-norm.

We consider the integral operator of the form

$$\left(I_m^{\eta ,1}f\right)\left(x\right)={\displaystyle \frac{m}{x^{m(\eta +1)}}}\underset{0}{\overset{x}{\int }}{t}^{m(\eta +1)-1}f\left(t\right)dt,x>0,$$

which is a homeomorphism from $F_{p,\mu }$ onto $F_{p,\mu }$ and has inverse given by

$$\left({\left(I_m^{\eta ,1}\right)}^{-1}f\right)\left(x\right)=(\eta +1)f\left(x\right)+{\displaystyle \frac{1}{m}}x{\displaystyle \frac{df}{dx}}=\left(\eta +1+{\displaystyle \frac{1}{m}}x{\displaystyle \frac{d}{dx}}\right)f\left(x\right).$$

Indeed,

$$\left({\left(I_m^{\eta ,1}\right)}^{-1}{I}_m^{\eta ,1}f\right)\left(x\right)={\left(I_m^{\eta ,1}\right)}^{-1}\left({\displaystyle \frac{m}{x^{m(\eta +1)}}}\underset{0}{\overset{x}{\int }}{t}^{m(\eta +1)-1}f\left(t\right)dt\right)=$$

$$={\displaystyle \frac{m(\eta +1)}{x^{m(\eta +1)}}}\underset{0}{\overset{x}{\int }}{t}^{m(\eta +1)-1}f\left(t\right)dt+{\displaystyle \frac{1}{m}}x{\displaystyle \frac{d}{dx}}{\displaystyle \frac{m}{x^{m(\eta +1)}}}\underset{0}{\overset{x}{\int }}{t}^{m(\eta +1)-1}f\left(t\right)dt=$$

$$={\displaystyle \frac{m(\eta +1)}{x^{m(\eta +1)}}}\underset{0}{\overset{x}{\int }}{t}^{m(\eta +1)-1}f\left(t\right)dt-{\displaystyle \frac{m(\eta +1)}{x^{m(\eta +1)}}}\underset{0}{\overset{x}{\int }}{t}^{m(\eta +1)-1}f\left(t\right)dt+f\left(x\right)=f\left(x\right).$$

If $f\left(0\right)=0$, then

$${\left({\displaystyle \frac{x}{m}}{\displaystyle \frac{d}{dx}}\right)}^{-1}f\left(x\right)=m\underset{0}{\overset{x}{\int }}{\displaystyle \frac{1}{t}}f\left(t\right)dt=\left(I_m^{-1,1}f\right)\left(x\right).$$

Since

$$\left(I_m^{-1,1}{x}^{m(\eta +1)}f\right)\left(x\right)=x^{m(\eta +1)}\left(I_m^{\eta ,1}f\right)\left(x\right),$$

(7)

we can write

$$\lambda I+I_m^{\eta ,1}={\displaystyle \frac{1}{x^{m(\eta +1)}}}(\lambda I+I_m^{-1,1})x^{m(\eta +1)}.$$

(8)

That means that the spectrum of $(-I_m^{\eta ,1})$ on $F_{p,\mu }$ is identical to the spectrum of $(-I_m^{-1,1})$ $F_{p,\mu +m(\eta +1)}$. Bearing this in mind, in order to find the resolvent of $(-I_m^{\eta ,1})$, we consider the equation

$$(\lambda I+I_m^{-1,1})\phi =\psi ,\phi ,\psi \in F_{p,\mu +m(\eta +1)}.$$

(9)

Solution to the equation (9) gives the resolvent of the operator $(-I_m^{-1,1})$:

$$\phi ={(\lambda I+I_m^{-1,1})}^{-1}\psi =R(\lambda ,-I_m^{-1,1})\psi .$$

Starting with finding the resolvent $R(\lambda ,-I_m^{-1,1})$. Applying ${\displaystyle \frac{x}{m}}{\displaystyle \frac{d}{dx}}$ to each side of (9) and dividing by $\lambda$ produces

$$\left({\displaystyle \frac{x}{m}}{\displaystyle \frac{d}{dx}}+{\displaystyle \frac{1}{\lambda }}\right)\phi ={\displaystyle \frac{1}{\lambda }}\left({\displaystyle \frac{x}{m}}{\displaystyle \frac{d}{dx}}\right)\psi .$$

This can be written as

$${\displaystyle \frac{1}{x^{\frac{m}{\lambda }}}}\left({\displaystyle \frac{x}{m}}{\displaystyle \frac{d}{dx}}\right)x^{{\displaystyle \frac{m}{\lambda }}}\phi ={\displaystyle \frac{1}{\lambda }}\left({\displaystyle \frac{x}{m}}{\displaystyle \frac{d}{dx}}\right)\psi \Rightarrow \left({\displaystyle \frac{x}{m}}{\displaystyle \frac{d}{dx}}\right)x^{{\displaystyle \frac{m}{\lambda }}}\phi ={\displaystyle \frac{x^{\frac{m}{\lambda }}}{\lambda }}\left({\displaystyle \frac{x}{m}}{\displaystyle \frac{d}{dx}}\right)\psi ,$$

and we get

$$\phi ={\displaystyle \frac{1}{\lambda x^{\frac{m}{\lambda }}}}{\left({\displaystyle \frac{x}{m}}{\displaystyle \frac{d}{dx}}\right)}^{-1}{x}^{{\displaystyle \frac{m}{\lambda }}}\left({\displaystyle \frac{x}{m}}{\displaystyle \frac{d}{dx}}\right)\psi ={\displaystyle \frac{1}{\lambda x^{\frac{m}{\lambda }}}}\left(I_m^{-1,1}{x}^{{\displaystyle \frac{m}{\lambda }}}\left({\displaystyle \frac{x}{m}}{\displaystyle \frac{d}{dx}}\right)\psi \right)\left(x\right).$$

Integrating by parts we obtain

$$\left(I_m^{-1,1}{x}^{{\displaystyle \frac{m}{\lambda }}}\left({\displaystyle \frac{x}{m}}{\displaystyle \frac{d}{dx}}\right)\psi \right)\left(x\right)=\underset{0}{\overset{x}{\int }}{t}^{{\displaystyle \frac{m}{\lambda }}}{\psi }^{\prime }\left(t\right)dt=x^{{\displaystyle \frac{m}{\lambda }}}\psi \left(x\right)-{\displaystyle \frac{m}{\lambda }}\left(I_m^{-1,1}{x}^{{\displaystyle \frac{m}{\lambda }}}\psi \right)\left(x\right)$$

and

$$\phi ={\displaystyle \frac{1}{\lambda }}\psi -{\displaystyle \frac{m}{{\lambda }^2{x}^{\frac{m}{\lambda }}}}{I}_m^{-1,1}{x}^{{\displaystyle \frac{m}{\lambda }}}\psi ={\displaystyle \frac{1}{\lambda }}\psi -{\displaystyle \frac{m}{{\lambda }^2{x}^{\frac{m}{\lambda }}}}{\left({\displaystyle \frac{x}{m}}{\displaystyle \frac{d}{dx}}\right)}^{-1}{x}^{{\displaystyle \frac{m}{\lambda }}}\psi ,$$

provided that $\mathrm{Re}(\mu +m(\eta +1)+1/\lambda )\ne 1/p$. When $\mathrm{Re}(\mu +m(\eta +1)+1/\lambda )>1/p$ by (7) we obtain

$$\phi ={\displaystyle \frac{1}{\lambda }}\psi -{\displaystyle \frac{m}{{\lambda }^2}}{I}_m^{{\displaystyle \frac{1}{\lambda }}-1,1}\psi$$

and

$$R(\lambda ,-I_m^{-1,1})={(\lambda I+I_m^{-1,1})}^{-1}={\displaystyle \frac{1}{\lambda }}I-{\displaystyle \frac{m}{{\lambda }^2}}{I}_m^{{\displaystyle \frac{1}{\lambda }}-1,1}.$$

Consequently, from (8), for $\lambda >0$, $\mathrm{Re}(\mu +m(\eta +1\left)\right)>1/p$ and $f\in F_{p,\mu }$,

$$(\lambda I+I_m^{\eta ,1})g={\displaystyle \frac{1}{x^{m(\eta +1)}}}(\lambda I+I_m^{-1,1})x^{m(\eta +1)}g=f\Rightarrow$$

$$\Rightarrow g={\displaystyle \frac{1}{x^{m(\eta +1)}}}{(\lambda I+I_m^{-1,1})}^{-1}{x}^{m(\eta +1)}f\Rightarrow$$

$$\Rightarrow g={\displaystyle \frac{1}{\lambda }}f-{\displaystyle \frac{m}{{\lambda }^2}}{\displaystyle \frac{1}{x^{m(\eta +1)}}}{I}_m^{{\displaystyle \frac{1}{\lambda }}-1,1}{x}^{m(\eta +1)}f.$$

We have

$${\displaystyle \frac{1}{x^{m(\eta +1)}}}\left(I_m^{{\displaystyle \frac{1}{\lambda }}-1,1}{x}^{m(\eta +1)}f\right)\left(x\right)={\displaystyle \frac{m}{x^{m\left(\eta +\frac{1}{\lambda }+1\right)}}}\underset{0}{\overset{x}{\int }}{t}^{m\left(\eta +{\displaystyle \frac{1}{\lambda }}+1\right)-1}f\left(t\right)dt=\left(I_m^{\eta +{\displaystyle \frac{1}{\lambda }},1}f\right)\left(x\right).$$

Finally,

$$R(\lambda ,-I_m^{\eta ,1})={\displaystyle \frac{1}{\lambda }}I-{\displaystyle \frac{m}{{\lambda }^2}}{I}_m^{\eta +{\displaystyle \frac{1}{\lambda }},1}\iff I-\lambda R(\lambda ,-I_m^{\eta ,1})={\displaystyle \frac{m}{\lambda }}{I}_m^{{\displaystyle \frac{1}{\eta +\lambda }},1}.$$

We obtain

$$R(\lambda ,-I_m^{\eta ,1})I_m^{\eta ,1}=I-\lambda R(\lambda ,-I_m^{\eta ,1})={\displaystyle \frac{m}{\lambda }}{I}_m^{{\displaystyle \frac{1}{\eta +\lambda }},1}.$$

(10)

We get by (1), taking into account (10)

$$\left({\left(I_m^{\eta ,1}\right)}^{\alpha }f\right)\left(x\right)={\displaystyle \frac{sin\left(\alpha \pi \right)}{\pi }}\underset{0}{\overset{\infty }{\int }}{\lambda }^{\alpha -1}R(\lambda ,-I_m^{\eta ,1})I_m^{\eta ,1}f\left(x\right)d\lambda =$$

$$={\displaystyle \frac{sin\left(\alpha \pi \right)m}{\pi }}\underset{0}{\overset{\infty }{\int }}{\lambda }^{\alpha -2}{I}_m^{{\displaystyle \frac{1}{\lambda }}+\eta ,1}f\left(x\right)d\lambda =$$

$${\displaystyle \frac{sin\left(\alpha \pi \right)m^2}{\pi }}\underset{0}{\overset{\infty }{\int }}{\lambda }^{\alpha -2}{x}^{-m\left(\eta +{\displaystyle \frac{1}{\lambda }}+1\right)}\left(\underset{0}{\overset{x}{\int }}{t}^{m\left(\eta +{\displaystyle \frac{1}{\lambda }}+1\right)-1}f\left(t\right)dt\right)d\lambda =$$

$$={\displaystyle \frac{sin\left(\alpha \pi \right)m^2}{\pi x^{m(\eta +1)}}}\underset{0}{\overset{x}{\int }}f\left(t\right)t^{n(\eta +1)-1}dt\underset{0}{\overset{\infty }{\int }}{\lambda }^{\alpha -2}{\left({\displaystyle \frac{x}{t}}\right)}^{-{\displaystyle \frac{m}{\lambda }}}d\lambda =$$

$$={\displaystyle \frac{sin\left(\alpha \pi \right)m^2}{\pi x^{m(\eta +1)}}}\underset{0}{\overset{x}{\int }}f\left(t\right)t^{n(\eta +1)-1}dt\underset{0}{\overset{\infty }{\int }}{\lambda }^{\alpha -2}{e}^{-{\displaystyle \frac{m}{\lambda }}log{\displaystyle \frac{x}{t}}}d\lambda =\left\{u={\displaystyle \frac{m}{\lambda }}log{\displaystyle \frac{x}{t}}\right\}=$$

$$={\displaystyle \frac{sin\left(\alpha \pi \right)m^2}{\pi x^{m(\eta +1)}}}\underset{0}{\overset{x}{\int }}f\left(t\right)t^{n(\eta +1)-1}{\left(mlog{\displaystyle \frac{x}{t}}\right)}^{\alpha -1}dt\underset{0}{\overset{\infty }{\int }}{u}^{-\alpha }{e}^{-u}du=$$

$$={\displaystyle \frac{sin\left(\alpha \pi \right)\mathrm{\Gamma }(1-\alpha )m^2}{\pi x^{m(\eta +1)}}}\underset{0}{\overset{x}{\int }}f\left(t\right)t^{n(\eta +1)-1}{\left(mlog{\displaystyle \frac{x}{t}}\right)}^{\alpha -1}dt=$$

$$={\displaystyle \frac{m^2}{\mathrm{\Gamma }\left(\alpha \right)x^{m(\eta +1)}}}\underset{0}{\overset{x}{\int }}f\left(t\right)t^{n(\eta +1)-1}{\left(mlog{\displaystyle \frac{x}{t}}\right)}^{\alpha -1}dt.$$

So we get generalisation of Hadamard fractional integral

$$\left({\left(I_m^{\eta ,1}\right)}^{\alpha }f\right)\left(x\right)={\displaystyle \frac{m^2}{\mathrm{\Gamma }\left(\alpha \right)x^{m(\eta +1)}}}\underset{0}{\overset{x}{\int }}f\left(t\right)t^{n(\eta +1)-1}{\left(mlog{\displaystyle \frac{x}{t}}\right)}^{\alpha -1}dt.$$

(11)

If we put $m=1$ in (11) we get the following operator

$$\left({\left(I_1^{\eta ,1}\right)}^{\alpha }f\right)\left(x\right)={\displaystyle \frac{1}{\mathrm{\Gamma }\left(\alpha \right)}}\underset{0}{\overset{x}{\int }}f\left(t\right){\left(log{\displaystyle \frac{x}{t}}\right)}^{\alpha -1}{\left({\displaystyle \frac{t}{x}}\right)}^{\eta +1}{\displaystyle \frac{dt}{t}}=\left(H^{\eta ,\alpha }f\right)\left(x\right).$$

This example in the particular case, when $m=1$ was given in [1]. Such integral was considered in [21,22]. If additionally we put $\eta =-1$ in (11) we obtain classical Hadamard fractional integral (see [2]) which is realised negative power of the Cauchy-Euler operator $x{\displaystyle \frac{d}{dx}}$

$$\left({\left(I_1^{-1,1}\right)}^{\alpha }f\right)\left(x\right)={\displaystyle \frac{1}{\mathrm{\Gamma }\left(\alpha \right)}}\underset{0}{\overset{x}{\int }}f\left(t\right){\left(log{\displaystyle \frac{x}{t}}\right)}^{\alpha -1}{\displaystyle \frac{dt}{t}}={\left(x{\displaystyle \frac{d}{dx}}\right)}^{-\alpha }f\left(x\right).$$

It is known that for $\alpha >0$ and $\eta >-1$ and for an infinitely differentiable function $f\left(x\right)$ defined for $x>0$ such that its Taylor series converges at every point $x>0$ the following equality holds (see [21])

$$\left({\left(I_1^{\eta ,1}\right)}^{\alpha }f\right)\left(x\right)={\displaystyle \frac{1}{\mathrm{\Gamma }\left(\alpha \right)}}\underset{0}{\overset{x}{\int }}{\left({\displaystyle \frac{t}{x}}\right)}^{\eta +1}{\left(ln{\displaystyle \frac{x}{t}}\right)}^{\alpha -1}f\left(t\right){\displaystyle \frac{dt}{t}}=\sum _{k=0}^{\infty }{\left\{{\displaystyle \genfrac{}{}{0pt}{}{-\alpha }{k}}\right\}}_{\eta +1}{x}^k{f}^{\left(k\right)}\left(x\right),$$

where

$${\left\{{\displaystyle \genfrac{}{}{0pt}{}{\alpha }{k}}\right\}}_c={\displaystyle \frac{1}{k!}}\sum _{j=1}^k{(-1)}^{k+j}\left({\displaystyle \genfrac{}{}{0pt}{}{k}{j}}\right){(c+j)}^{\alpha },\alpha ,c\in \mathbb{R},k\in \mathbb{N}\cup \left\{0\right\}$$

(12)

are the generalised Stirling numbers of the second kind.

3. The Matrix Approach to the Fractional Powers of Operator

This section is devoted to the matrix approach for calculating fractional powers of operators. We start from the illustrative example giving the fractional power of the simplest integral operator, then we construct a power of a diagonalizable triangle matrix and its resolvent.

3.1. Illustrative Example

To illustrate the concept of using matrices for introducing new forms of fractional integro-differentiation, starting with the simplest case and consider a Riemann integral $\left({\mathcal{I}}_{a+}^{1}f\right)\left(x\right)=\underset{a}{\overset{x}{\int }}f\left(t\right)dt$, where $a<x\le b$, $f\in C[a,b]$.

Consider a uniform partition $P=(x_0,x_1,\dots ,x_n)$ of the interval $[a,x]$, defined by points $x_p=a+ph$ for $p=0,1,2,...,n$, where the step size is $h={\displaystyle \frac{x-a}{n}}$. Then the integral ${\mathcal{I}}_{a+}^1$ is defined as the limit of the corresponding integral sums as $n\to \infty$:

$$\left({\mathcal{I}}_{a+}^{1}f\right)\left(x\right)=\underset{a}{\overset{x}{\int }}f\left(t\right)dt=\underset{n\to +\infty }{lim}\sum _{p=0}^{n-1}f\left({\xi }_p\right)h=\underset{h\to +0}{lim}\sum _{p=0}^{n-1}f\left({\xi }_p\right)h,$$

where ${\xi }_p$ is an arbitrary number ${\xi }_p\in [x-(p+1)h,x-ph]$, $p=0,...,n-1$.

If we take ${\xi }_p=x-ph$ and denote $f_p=f(x-ph)$, $p=0,...,n-1$, $f_n=0$, we obtain

$$\left({\mathcal{I}}_{a+}^{1}f\right)\left(x\right)=\underset{a}{\overset{x}{\int }}f\left(t\right)dt=\underset{n\to +\infty }{lim}h\sum _{p=0}^{n-1}f(x-ph)=\underset{h\to +0}{lim}h\sum _{p=0}^{n-1}{f}_p,h={\displaystyle \frac{x-a}{n}}.$$

Approximate operator is

$$\left({\mathbb{I}}_{a+,n}^{1}f\right)\left(x\right)=h\sum _{p=0}^{n-1}f(x-ph),h={\displaystyle \frac{x-a}{n}}.$$

(13)

Let ${\mathbf{f}}_n^{\prime }=(f_0,f_1,...,f_{n-1})$. Operator ${\mathcal{I}}_{a+}$ we associate with the $\left(n\times n\right)$ matrix

$${\mathcal{A}}_n=h\left(\begin{array}{ccccc}1& 1& \dots & 1& 1\\ 0& 1& \dots & 1& 1\\ ⋮& ⋮& \dots & ⋮& ⋮\\ 0& 0& \dots & 1& 1\\ 0& 0& \dots & 0& 1\end{array}\right).$$

(14)

in the sense

$$\underset{h\to +0}{lim}{\mathcal{A}}_n·{\mathbf{f}}_n=$$

$$=\underset{h\to +0}{lim}h\left(\begin{array}{ccccc}1& 1& \dots & 1& 1\\ 0& 1& \dots & 1& 1\\ ⋮& ⋮& \dots & ⋮& ⋮\\ 0& 0& \dots & 1& 1\\ 0& 0& \dots & 0& 1\end{array}\right)·\left(\begin{array}{c}{f}_0\\ f_1\\ ⋮\\ f_{n-2}\\ f_{n-1}\end{array}\right)=\left(\begin{array}{c}\underset{h\to +0}{lim}h\sum _{p=0}^{n-1}{f}_p\\ \underset{h\to +0}{lim}h\sum _{p=1}^{n-1}{f}_p\\ ⋮\\ \underset{h\to +0}{lim}h\sum _{p=n-2}^{n-1}{f}_p\\ \underset{h\to +0}{lim}h{f}_{n-1}\end{array}\right)=\left(\begin{array}{c}\underset{a}{\overset{x}{\int }}f\left(x\right)dx\\ \underset{a_1}{\overset{x}{\int }}f\left(x\right)dx\\ ⋮\\ \underset{a_{n-2}}{\overset{x}{\int }}f\left(x\right)dx\\ \underset{a_{n-1}}{\overset{x}{\int }}f\left(x\right)dx\end{array}\right).$$

We can see that ${\mathcal{I}}_{a+}^1$ can be represented by

$$\left({\mathcal{I}}_{a+}^{1}f\right)\left(x\right)=\underset{n\to +\infty }{lim}({\mathcal{A}}_n·{\mathbf{f}}_n)·I_{1n},h={\displaystyle \frac{x-a}{n}},$$

(15)

where $I_{1n}={\underbrace{(1,0,...,0)}}_n$. We have so called "matrix operator" ${\mathcal{A}}_n$ and the limit when $h\to +0$ of the first row of ${\mathcal{A}}_n·{\mathbf{f}}_n$ gives $\left({\mathcal{I}}_{a+}^{1}f\right)\left(x\right)$.

We explore the concept of associating integral and differential operators with their corresponding matrices. Our goal is to construct the fractional power ${\mathcal{A}}_n^{\alpha }$, $\alpha \in \mathbb{R}$, of the matrix operator ${\mathcal{A}}_n$. Subsequently, we will derive the fractional power of an operator as a limit by the formula

$$\left({\mathcal{I}}_{a+}^{\alpha }f\right)\left(x\right)=\underset{n\to \infty }{lim}({\mathcal{A}}_n^{\alpha }·{\mathbf{f}}_n)·I_{1n}=\underset{h\to +0}{lim}({\mathcal{A}}_n^{\alpha }·{\mathbf{f}}_n)·I_{1n},$$

(16)

where $A_n^{\alpha }$ is first row of ${\mathcal{A}}_n^{\alpha }$ for all $\alpha \in \mathbb{R}$.

In order to construct the power ${\mathcal{A}}_n^{\alpha }$, $\alpha \in \mathbb{R}$ we will use the formula ${\mathcal{A}}_n^{\alpha }=e^{\alpha ln{A}_n}$. Choosing for uniqueness the main branch of the logarithm $ln{\mathcal{A}}_n$, we obtain

$$ln{\mathcal{A}}_n=lnh+\left(\begin{array}{cccccc}0& 1& {\displaystyle \frac{1}{2}}& \dots & {\displaystyle \frac{1}{n-2}}& {\displaystyle \frac{1}{n-1}}\\ 0& 0& 1& \dots & {\displaystyle \frac{1}{n-3}}& {\displaystyle \frac{1}{n-2}}\\ ⋮& ⋮& ⋮& \dots & ⋮& ⋮\\ 0& 0& 0& \dots & 0& 1\\ 0& 0& 0& \dots & 0& 0\end{array}\right),$$

$$\alpha ln{\mathcal{A}}_n=\alpha lnh+\left(\begin{array}{cccccc}0& \alpha & {\displaystyle \frac{\alpha }{2}}& \dots & {\displaystyle \frac{\alpha }{n-2}}& {\displaystyle \frac{\alpha }{n-1}}\\ 0& 0& \alpha & \dots & {\displaystyle \frac{\alpha }{n-3}}& {\displaystyle \frac{\alpha }{n-2}}\\ ⋮& ⋮& ⋮& \dots & ⋮& ⋮\\ 0& 0& 0& \dots & 0& \alpha \\ 0& 0& 0& \dots & 0& 0\end{array}\right),$$

$${\mathcal{A}}_n^{\alpha }=h^{\alpha }\left(\begin{array}{cccccc}1& \alpha & {\displaystyle \frac{\alpha (\alpha +1)}{2}}& \dots & {\displaystyle \frac{\alpha (\alpha +1)...(\alpha +n-3)}{(n-2)!}}& {\displaystyle \frac{\alpha (\alpha +1)...(\alpha +n-2)}{(n-1)!}}\\ 0& 1& \alpha & \dots & {\displaystyle \frac{\alpha (\alpha +1)...(\alpha +n-4)}{(n-3)!}}& {\displaystyle \frac{\alpha (\alpha +1)...(\alpha +n-3)}{(n-2)!}}\\ ⋮& ⋮& ⋮& \dots & ⋮& ⋮\\ 0& 0& 0& \dots & 1& \alpha \\ 0& 0& 0& \dots & 0& 1\end{array}\right).$$

(17)

We get general power (17) of matrix operator ${\mathcal{A}}_n$ for all $\alpha \in \mathbb{R}$ which automatically has semigroup property.

It is easy to see that when $\alpha =-1$

$${\mathcal{A}}_n^{-1}={\displaystyle \frac{1}{h}}\left(\begin{array}{cccccc}1& -1& 0& \dots & 0& 0\\ 0& 1& -1& \dots & 0& 0\\ ⋮& ⋮& ⋮& \dots & ⋮& ⋮\\ 0& 0& 0& \dots & 1& -1\\ 0& 0& 0& \dots & 0& 1\end{array}\right),$$

we get the formula of matrix ${\mathcal{A}}_n$ inversion and

$$\underset{h\to +0}{lim}{\mathcal{A}}_n^{-1}·{\mathbf{f}}_n=$$

$$=\underset{h\to +0}{lim}{\displaystyle \frac{1}{h}}\left(\begin{array}{cccccc}1& -1& 0& \dots & 0& 0\\ 0& 1& -1& \dots & 0& 0\\ ⋮& ⋮& ⋮& \dots & ⋮& ⋮\\ 0& 0& 0& \dots & 1& -1\\ 0& 0& 0& \dots & 0& 1\end{array}\right)·\left(\begin{array}{c}{f}_0\\ f_1\\ ⋮\\ f_{n-2}\\ f_{n-1}\end{array}\right)=\underset{h\to +0}{lim}{\displaystyle \frac{1}{h}}\left(\begin{array}{c}{f}_0-f_1\\ f_1-f_2\\ ⋮\\ f_{n-2}-f_{n-1}\\ f_{n-1}-f_n\end{array}\right)=$$

$$=\left(\begin{array}{c}\underset{h\to +0}{lim}{\displaystyle \frac{f\left(x\right)-f(x-h)}{h}}\\ \underset{h\to +0}{lim}{\displaystyle \frac{f(x-h)-f(x-2h)}{h}}\\ ⋮\\ \underset{h\to +0}{lim}{\displaystyle \frac{f(x-(n-2\left)h\right)-f(x-(n-1\left)h\right)}{h}}\\ \underset{h\to +0}{lim}{\displaystyle \frac{f(x-(n-1\left)h\right)-f(x-nh)}{h}}\end{array}\right)=\left(\begin{array}{c}{f}_{-}^{\prime }\left(x\right)\\ f_{-}^{\prime }\left(x_1\right)\\ ⋮\\ f_{-}^{\prime }\left(x_{n-2}\right)\\ f_{-}^{\prime }\left(x_{n-1}\right)\end{array}\right)$$

and

$$\left({\mathcal{I}}_{a+}^{-1}f\right)\left(x\right)=\underset{n\to +\infty }{lim}({\mathcal{A}}_n^{-1}·{\mathbf{f}}_n)·I_{1n}=\underset{h\to +0}{lim}{\displaystyle \frac{f\left(x\right)-f(x-h)}{h}}=f_{-}^{\prime }\left(x\right)$$

(18)

is a left derivative of f at x.

Next we notice that when $\alpha =-m$, $m\in \mathbb{N}$, $m\le n$ we get that in the first row all elements with a column number greater than $(m+2)$ are equal to zero

$${\mathcal{A}}_n^{-m}={\displaystyle \frac{1}{h^m}}\left(\begin{array}{cccccccc}1& -m& {\displaystyle \frac{m(m-1)}{2}}& \dots & {(-1)}^m{\displaystyle \frac{m(m-1)...1}{m!}}& 0& \dots & 0\\ 0& 1& -m& \dots & {(-1)}^{m-1}{\displaystyle \frac{m(m-1)...2}{(m+1)!}}& 0& \dots & 0\\ ⋮& ⋮& ⋮& \dots & ⋮& ⋮& ⋮& ⋮\\ 0& 0& 0& \dots & 0& 1& -m& {\displaystyle \frac{m(m-1)}{2}}\\ 0& 0& 0& \dots & 0& 0& 1& -m\\ 0& 0& 0& \dots & 0& 0& 0& 1\end{array}\right),$$

(19)

so for $\alpha =-m$, $m\in \mathbb{N}$, $m\le n$ by (16)

$$\begin{array}{c}\left({\mathcal{I}}_{a+}^{-m}f\right)\left(x\right)=\underset{h\to +0}{lim}({\mathcal{A}}_n^{-m}·{\mathbf{f}}_n)·I_{1n}=\underset{h\to +0}{lim}\sum _{p=0}^m{\displaystyle \frac{{(-1)}^p}{h^m}}\left({\displaystyle \genfrac{}{}{0pt}{}{m}{p}}\right)f_p=\hfill \\ \hfill =\underset{h\to +0}{lim}\sum _{p=0}^m{\displaystyle \frac{{(-1)}^p}{h^m}}\left({\displaystyle \genfrac{}{}{0pt}{}{m}{p}}\right)f(x-ph)=\underset{h\to +0}{lim}{\displaystyle \frac{{\nabla }_h^m\left[f\right]\left(x\right)}{h^m}}=f_{-}^{\left(m\right)}\left(x\right),\end{array}$$

where $\left({\displaystyle \genfrac{}{}{0pt}{}{m}{p}}\right)={\displaystyle \frac{m!}{p!(m-p)!}}$ are binomial coefficients, ${\nabla }_h^m\left[f\right]\left(x\right)=\sum _{p=0}^m{(-1)}^p\left({\displaystyle \genfrac{}{}{0pt}{}{m}{p}}\right)f(x-ph)$ is a backward difference, $f_{-}^{\left(m\right)}\left(x\right)$ is a left derivative of f at x of order $m\in \mathbb{N}$.

It $\alpha \in \mathbb{R}$, $\alpha \ne -m$, $m\in \mathbb{N}$, then, Taking into account that

$$\alpha (\alpha +1)...(\alpha +p-2)={\alpha }^{(p-1)}$$

is the rising factorial which can be written using the Gamma function by the formula

$${\alpha }^{(p-1)}={\displaystyle \frac{\mathrm{\Gamma }(\alpha +p-1)}{\mathrm{\Gamma }\left(\alpha \right)}},$$

we get

$${\mathcal{A}}_n^{\alpha }=h^{\alpha }\left(\begin{array}{ccccc}1& {\displaystyle \frac{\mathrm{\Gamma }(\alpha +1)}{\mathrm{\Gamma }\left(\alpha \right)\mathrm{\Gamma }\left(2\right)}}& {\displaystyle \frac{\mathrm{\Gamma }(\alpha +2)}{\mathrm{\Gamma }\left(\alpha \right)\mathrm{\Gamma }\left(3\right)}}& \dots & {\displaystyle \frac{\mathrm{\Gamma }(\alpha +n-1)}{\mathrm{\Gamma }\left(\alpha \right)\mathrm{\Gamma }\left(n\right)}}\\ 0& 1& {\displaystyle \frac{\mathrm{\Gamma }(\alpha +1)}{\mathrm{\Gamma }\left(\alpha \right)\mathrm{\Gamma }\left(2\right)}}& \dots & {\displaystyle \frac{\mathrm{\Gamma }(\alpha +n-2)}{\mathrm{\Gamma }\left(\alpha \right)\mathrm{\Gamma }(n-1)}}\\ ⋮& ⋮& ⋮& \dots & ⋮\\ 0& 0& 0& \dots & 1\end{array}\right).$$

(20)

Therefore, by (16) we can write

$$\left({\mathcal{I}}_{a+}^{\alpha }f\right)\left(x\right)=\underset{h\to +0}{lim}({\mathcal{A}}_n^{\alpha }·{\mathbf{f}}_n)·I_{1n}=\underset{h\to +0}{lim}{\displaystyle \frac{h^{\alpha }}{\mathrm{\Gamma }\left(\alpha \right)}}\sum _{p=0}^{n-1}{\displaystyle \frac{\mathrm{\Gamma }(\alpha +p)}{\mathrm{\Gamma }(p+1)}}{f}_p,$$

(21)

where $f_p=f\left(x-ph\right)$, $h={\displaystyle \frac{x-a}{n}}$. Formula (21) gives the left-sided Grünwald-Letnikov fractional integro-differentiation of order $\alpha \in \mathbb{R}$, $\alpha \ne -m$, $m\in \mathbb{N}$ (see [2]).

It is knows that (see [3]) for $f\left(x\right)\in L_1(a,b)$ limit (21) exists for almost everyone x and

$$\left({\mathcal{I}}_{a+}^{\alpha }f\right)\left(x\right)=\left(I_{a+}^{\alpha }f\right)\left(x\right),\left({\mathcal{D}}_{a+}^{\alpha }f\right)\left(x\right)=\left(D_{a+}^{\alpha }f\right)\left(x\right)$$

where

$$\left(I_{a+}^{\alpha }f\right)\left(x\right)={\displaystyle \frac{1}{\mathrm{\Gamma }\left(\alpha \right)}}\underset{a}{\overset{x}{\int }}{\displaystyle \frac{f\left(t\right)dt}{{(x-t)}^{1-\alpha }}}$$

is the left-sided Riemann-Liouville integral,

$$\left(D_{a+}^{\alpha }f\right)\left(x\right)={\displaystyle \frac{1}{\mathrm{\Gamma }(p-\alpha )}}{\left({\displaystyle \frac{d}{dx}}\right)}^p\underset{a}{\overset{x}{\int }}{\displaystyle \frac{f\left(t\right)dt}{{(x-t)}^{\alpha -p+1}}},p=\left[\alpha \right]+1$$

is the left-sided Riemann-Liouville derivative.

The finite-difference operator

$$\left({\mathbb{I}}_{a+,n}^{\alpha }f\right)\left(x\right)={\displaystyle \frac{h^{\alpha }}{\mathrm{\Gamma }\left(\alpha \right)}}\sum _{p=0}^{n-1}{\displaystyle \frac{\mathrm{\Gamma }(\alpha +p)}{\mathrm{\Gamma }(p+1)}}{f}_p,f_p=f\left(x-ph\right),h={\displaystyle \frac{x-a}{n}}$$

(22)

which is a prelimiting form of (21) will be used for approximation of the left-sided Riemann-Liouville integral. Let us note, that in [3] this approximation was used for the numerical calculation for concrete engineering problems.

3.2. Power of a Diagonalizable Triangle Matrix

In order to construct fractional power of the operator $\mathcal{A}$ we shall construct the power of matrix $A^{\alpha }=e^{\alpha lnA}$, where we choose the main branch of the logarithm for uniqueness. In Section 3.1 we constructed the $A^{\alpha }$ when all $a_k=1$, $k=0,1,...,N-1$. In this section we construct $A^{\alpha }=e^{\alpha lnA}$ when $a_0\ne a_1\ne ...\ne a_{N-1}$. In this case a matrix A is diagonalizable.

Theorem 3.1.

Let $\det A\ne 0$, $a_1\ne a_2\ne ...\ne a_n$, then the matrix

$$A=\left(\begin{array}{cccc}{a}_1& a_2& \dots & a_n\\ 0& a_2& \dots & a_n\\ ⋮& ⋮& \ddots & ⋮\\ 0& 0& \dots & a_n\end{array}\right)$$

is diagonalizable and eigenvalue decomposition for matrix A is

$$A=PD{P}^{-1},$$

(23)

where $P={\left(b_{sm}\right)}_{m=1,s=1}^{n,n}$, $P^{-1}={\left(c_{sm}\right)}_{m=1,s=1}^{n,n}$, $D={\left(d_{sm}\right)}_{m=1,s=1}^{n,n}$,

$$b_{sm}=\left\{\begin{array}{cc}1\hfill & ifs=m;\hfill \\ 0\hfill & ifm<s;\hfill \\ {\displaystyle \frac{a_m^{m-s}}{\prod _{i=s}^{m-1}(a_m-a_i)}}\hfill & ifs<m,\hfill \end{array}\right.c_{sm}=\left\{\begin{array}{cc}1\hfill & ifs=m;\hfill \\ 0\hfill & ifm<s;\hfill \\ {\displaystyle \frac{a_s^{m-s-1}{a}_m}{\prod _{i=s+1}^m(a_s-a_i)}}\hfill & ifs<m,\hfill \end{array}\right.$$

$$d_{sm}=\left\{\begin{array}{cc}{a}_m\hfill & ifm=s;\hfill \\ 0\hfill & ifs\ne m.\hfill \end{array}\right.$$

Proof. Eigenvalues for A are $a_1,a_2,...,a_n$. Since $a_1\ne a_2\ne ...\ne a_n$ then matrix A is diagonalizable and can be written in the form $A=PD{P}^{-1}$, where

$$D=\left(\begin{array}{cccc}{a}_1& 0& \dots & 0\\ 0& a_2& \dots & 0\\ ⋮& ⋮& \ddots & ⋮\\ 0& 0& \dots & a_n\end{array}\right)$$

is a diagonal matrix and

$$P=\left(\begin{array}{cccc}{b}_{11}& b_{12}& \dots & b_{1n}\\ 0& b_{22}& \dots & b_{2n}\\ ⋮& ⋮& \ddots & ⋮\\ 0& 0& \dots & b_{nn}\end{array}\right)$$

is an invertible matrix. Now, we can find P and $P^{-1}$.

Since

$$PD=\left(\begin{array}{cccc}{b}_{11}& b_{12}& \dots & b_{1n}\\ 0& b_{22}& \dots & b_{2n}\\ ⋮& ⋮& \ddots & ⋮\\ 0& 0& \dots & b_{nn}\end{array}\right)·\left(\begin{array}{cccc}{a}_1& 0& \dots & 0\\ 0& a_2& \dots & 0\\ ⋮& ⋮& \ddots & ⋮\\ 0& 0& \dots & a_n\end{array}\right)=\left(\begin{array}{cccc}{a}_1{b}_{11}& a_2{b}_{12}& \dots & a_n{b}_{1n}\\ 0& a_2{b}_{22}& \dots & a_n{b}_{2n}\\ ⋮& ⋮& \ddots & ⋮\\ 0& 0& \dots & a_n{b}_{nn}\end{array}\right),$$

then an equality $AP=PD$ gives

$$\left(\begin{array}{cccc}{a}_1& a_2& \dots & a_n\\ 0& a_2& \dots & a_n\\ ⋮& ⋮& \ddots & ⋮\\ 0& 0& \dots & a_n\end{array}\right)·\left(\begin{array}{cccc}{b}_{11}& b_{12}& \dots & b_{1n}\\ 0& b_{22}& \dots & b_{2n}\\ ⋮& ⋮& \ddots & ⋮\\ 0& 0& \dots & b_{nn}\end{array}\right)=\left(\begin{array}{cccc}{a}_1{b}_{11}& a_2{b}_{12}& \dots & a_n{b}_{1n}\\ 0& a_2{b}_{22}& \dots & a_n{b}_{2n}\\ ⋮& ⋮& \ddots & ⋮\\ 0& 0& \dots & a_n{b}_{nn}\end{array}\right).$$

For the first row:

$$\left\{\begin{array}{cc}{a}_1·b_{11}\hfill & =a_1·b_{11}\hfill \\ a_1·b_{12}+a_2·b_{22}\hfill & =a_2·b_{12}\hfill \\ ...\hfill \\ a_1·b_{1n}+a_2·b_{2n}+...+a_n·b_{nn}\hfill & =a_n·b_{1n}\hfill \end{array}\right.$$

for the second row:

$$\left\{\begin{array}{cc}{a}_2·b_{22}\hfill & =a_2·b_{22}\hfill \\ a_2·b_{23}+a_3·b_{33}\hfill & =a_3·b_{23}\hfill \\ ...\hfill \\ a_2·b_{2n}+a_3·b_{3n}+...+a_n·b_{nn}\hfill & =a_n·b_{2n}\hfill \end{array}\right.$$

for the n-th row

$$a_n·b_{nn}=a_n·b_{nn}.$$

Let s be the number of row, $s=1,...,n$, $m=1,..,n$, $s\le m$ we we can write the system

$$\sum _{i=s}^m{a}_i{b}_{im}=a_m{b}_{sm}.$$

For $s=m-1$

$$a_{m-1}{b}_{m-1m}+a_m{b}_{mm}=a_m{b}_{m-1m}$$

then

$$b_{m-1m}={\displaystyle \frac{a_m}{a_m-a_{m-1}}}{b}_{mm}.$$

If $b_{11}=b_{22}=...=b_{nn}=1$, then

$$b_{m-1m}={\displaystyle \frac{a_m}{a_m-a_{m-1}}}.$$

For $s=m-2$

$$a_{m-2}{b}_{m-2m}+a_{m-1}{b}_{m-1m}+a_m{b}_{mm}=a_m{b}_{m-2m},$$

then

$$b_{m-2m}={\displaystyle \frac{1}{a_m-a_{m-2}}}(a_{m-1}{b}_{m-1m}+a_m{b}_{mm})=$$

$$={\displaystyle \frac{1}{a_m-a_{m-2}}}\left({\displaystyle \frac{a_{m-1}{a}_m}{a_m-a_{m-1}}}+a_m\right)={\displaystyle \frac{a_m^2}{(a_m-a_{m-1})(a_m-a_{m-2})}}.$$

For $s=m-3$

$$a_{m-3}{b}_{m-3m}+a_{m-2}{b}_{m-2m}+a_{m-1}{b}_{m-1m}+a_m{b}_{mm}=a_m{b}_{m-3m},$$

then

$$b_{m-3m}={\displaystyle \frac{1}{a_m-a_{m-3}}}(a_{m-2}{b}_{m-2m}+a_{m-1}{b}_{m-1m}+a_m{b}_{mm})=$$

$$={\displaystyle \frac{1}{a_m-a_{m-3}}}\left({\displaystyle \frac{a_{m-2}{a}_m^2}{(a_m-a_{m-1})(a_m-a_{m-2})}}+{\displaystyle \frac{a_{m-1}{a}_m}{a_m-a_{m-1}}}+a_m\right)=$$

$$={\displaystyle \frac{a_m^3}{(a_m-a_{m-1})(a_m-a_{m-2})(a_m-a_{m-3})}}.$$

Therefore for $s<m$

$$b_{m-sm}={\displaystyle \frac{a_m^s}{\prod _{i=1}^s(a_m-a_{m-i})}}$$

or for $s<m$

$$b_{sm}={\displaystyle \frac{a_m^{m-s}}{\prod _{i=s}^{m-1}(a_m-a_i)}}$$

and

$$P=\left(\begin{array}{ccccc}1& {\displaystyle \frac{a_2}{a_2-a_1}}& {\displaystyle \frac{a_3^2}{(a_3-a_1)(a_3-a_2)}}& \dots & {\displaystyle \frac{a_n^{n-1}}{\prod _{i=1}^{n-1}(a_n-a_i)}}\\ 0& 1& {\displaystyle \frac{a_3}{a_3-a_2}}& \dots & {\displaystyle \frac{a_n^{n-2}}{\prod _{i=2}^{n-1}(a_n-a_i)}}\\ 0& 0& 1& \dots & {\displaystyle \frac{a_n^{n-3}}{\prod _{i=3}^{n-1}(a_n-a_i)}}\\ ⋮& ⋮& ⋮& \ddots & ⋮\\ 0& 0& \dots & \dots & 1\end{array}\right).$$

Finally, $P={\left(b_{sm}\right)}_{m=1,s=1}^{n,n}$, $s=1,...,n$, $m=1,...,n$ where

$$b_{sm}=\left\{\begin{array}{cc}1\hfill & ifs=m;\hfill \\ 0\hfill & ifm<s;\hfill \\ {\displaystyle \frac{a_m^{m-s}}{\prod _{i=s}^{m-1}(a_m-a_i)}}\hfill & ifs<m.\hfill \end{array}\right.$$

We now proceed to find the inverse matrix $P^{-1}$. Inverse of an lower triangular matrix P is another upper triangular matrix

$$P^{-1}=\left(\begin{array}{cccc}{c}_{11}& c_{12}& \dots & c_{1n}\\ 0& c_{22}& \dots & c_{2n}\\ ⋮& ⋮& \ddots & ⋮\\ 0& 0& \dots & c_{nn}\end{array}\right).$$

The elements $c_{sm}$ of $P^{-1}$ can be found from the system

$$P{P}^{-1}=\left(\begin{array}{cccc}{b}_{11}& b_{12}& \dots & b_{1n}\\ 0& b_{22}& \dots & b_{2n}\\ ⋮& ⋮& \ddots & ⋮\\ 0& 0& \dots & b_{nn}\end{array}\right)·\left(\begin{array}{cccc}{c}_{11}& c_{12}& \dots & c_{1n}\\ 0& c_{22}& \dots & c_{2n}\\ ⋮& ⋮& \ddots & ⋮\\ 0& 0& \dots & c_{nn}\end{array}\right)=\left(\begin{array}{cccc}1& 0& \dots & 0\\ 0& 1& \dots & 0\\ ⋮& ⋮& \ddots & ⋮\\ 0& 0& \dots & 1\end{array}\right).$$

It is clear that

$$b_{mm}{c}_{mm}=1,m=1,...,n$$

so diagonal elements of $P^{-1}$ are reciprocal of the elements of P:

$$c_{mm}={\displaystyle \frac{1}{b_{mm}}}=1,m=1,...,n.$$

For other elements we get

$$\sum _{i=s}^k{b}_{si}{c}_{ik}=0,s=1,...,n,k=s+1,...,n.$$

Putting $k=s+1$ we obtain

$$\sum _{i=s}^{s+1}{b}_{si}{c}_{is+1}=b_{ss}{c}_{ss+1}+b_{ss+1}{c}_{s+1s+1}=0,s=1,...,n-1$$

or

$$c_{ss+1}=-{\displaystyle \frac{b_{ss+1}}{b_{ss}}}{c}_{s+1s+1}=-b_{ss+1}=-{\displaystyle \frac{a_{s+1}^{s+1-s}}{\prod _{i=s}^{s+1-1}(a_{s+1}-a_i)}}={\displaystyle \frac{a_{s+1}}{a_s-a_{s+1}}},s=1,...,n-1.$$

For $k=s+2$ we get

$$\sum _{i=s}^{s+2}{b}_{si}{c}_{is+2}=b_{ss}{c}_{ss+2}+b_{ss+1}{c}_{s+1s+2}+b_{ss+2}{c}_{s+2s+2}=0,s=1,...,n-2,$$

$$c_{ss+2}=-{\displaystyle \frac{1}{b_{ss}}}(b_{ss+1}{c}_{s+1s+2}+b_{ss+2}{c}_{s+2s+2})=-(b_{ss+1}{c}_{s+1s+2}+b_{ss+2})=$$

$$=-\left(-{\displaystyle \frac{a_{s+1}{a}_{s+2}}{(a_{s+1}-a_s)(a_{s+2}-a_{s+1})}}+{\displaystyle \frac{a_{s+2}^2}{(a_{s+2}-a_s)(a_{s+2}-a_{s+1})}}\right)=$$

$$=-{\displaystyle \frac{a_s{a}_{s+2}(a_{s+1}-a_{s+2})}{(a_{s+1}-a_s)(a_{s+2}-a_s)(a_{s+2}-a_{s+1})}}={\displaystyle \frac{a_s{a}_{s+2}}{(a_{s+1}-a_s)(a_{s+2}-a_s)}}.$$

For $k=s+3$ we get

$$\begin{array}{c}\sum _{i=s}^{s+3}{b}_{si}{c}_{is+3}=b_{ss}{c}_{ss+3}+b_{ss+1}{c}_{s+1s+3}+b_{ss+2}{c}_{s+2s+3}+b_{ss+3}{c}_{s+3s+3}=0,s=1,...,n-3,\hfill \\ c_{ss+3}=-{\displaystyle \frac{1}{b_{ss}}}(b_{ss+1}{c}_{s+1s+3}+b_{ss+2}{c}_{s+2s+3}+b_{ss+3}{c}_{s+3s+3})=\\ =-(b_{ss+1}{c}_{s+1s+3}+b_{ss+2}{c}_{s+2s+3}+b_{ss+3})=\\ =-\left({\displaystyle \frac{a_{s+1}^2{a}_{s+3}}{(a_{s+1}-a_s)(a_{s+2}-a_{s+1})(a_{s+3}-a_{s+1})}}-{\displaystyle \frac{a_{s+2}^2{a}_{s+3}}{(a_{s+2}-a_s)(a_{s+2}-a_{s+1})(a_{s+3}-a_{s+2})}}+\right.\\ \left.+{\displaystyle \frac{a_{s+3}^3}{(a_{s+3}-a_s)(a_{s+3}-a_{s+1})(a_{s+3}-a_{s+2})}}\right)=\\ \hfill ={\displaystyle \frac{a_s^2{a}_{s+3}}{(a_s-a_{s+1})(a_s-a_{s+})(a_s-a_{s+3})}}.\end{array}$$

Therefore for $s=1,...,n$,

$$c_{ss+m}={\displaystyle \frac{a_s^{m-1}{a}_{s+m}}{\prod _{i=1}^m(a_s-a_{s+i})}},$$

$$c_{sm}={\displaystyle \frac{a_s^{m-s-1}{a}_m}{\prod _{i=1}^{m-s}(a_s-a_{s+i})}}={\displaystyle \frac{a_s^{m-s-1}{a}_m}{\prod _{i=s+1}^m(a_s-a_i)}}$$

and

$$P^{-1}=\left(\begin{array}{ccccc}1& {\displaystyle \frac{a_2}{a_1-a_2}}& {\displaystyle \frac{a_1{a}_3}{(a_1-a_2)(a_1-a_3)}}& \dots & {\displaystyle \frac{a_1^{n-2}{a}_n}{\prod _{i=2}^n(a_1-a_i)}}\\ 0& 1& {\displaystyle \frac{a_3}{a_2-a_3}}& \dots & {\displaystyle \frac{a_2^{n-3}{a}_n}{\prod _{i=3}^n(a_2-a_i)}}\\ 0& 0& 1& \dots & {\displaystyle \frac{a_3^{n-4}{a}_n}{\prod _{i=4}^n(a_3-a_i)}}\\ ⋮& ⋮& ⋮& \ddots & ⋮\\ 0& 0& \dots & \dots & 1\end{array}\right).$$

Finally, $P^{-1}={\left(c_{sm}\right)}_{m=1,s=1}^{n,n}$, where

$$c_{sm}=\left\{\begin{array}{cc}1\hfill & ifs=m;\hfill \\ 0\hfill & ifm<s;\hfill \\ {\displaystyle \frac{a_s^{m-s-1}{a}_m}{\prod _{i=s+1}^m(a_s-a_i)}}\hfill & ifs<m.\hfill \end{array}\right.$$

Theorem 3.2.

Let $\det A\ne 0$, $a_1\ne a_2\ne ...\ne a_n$, then the power $\alpha \in \mathbb{R}$ of matrix A is

$$A^{\alpha }={\left(p_{sm}\right)}_{m=1,s=1}^{n,n},$$

(24)

where

$$p_{sm}=\left\{\begin{array}{cc}0\hfill & ifm<s;\hfill \\ \sum _{i=s}^m{\displaystyle \frac{a_i^{\alpha +m-s-1}{a}_m}{\prod _{j=s,j\ne i}^m(a_j-a_i)}}\hfill & ifs<m;\hfill \\ a_m^{\alpha }\hfill & ifm=s.\hfill \end{array}\right.$$

(25)

Proof. By Theorem 3.1 we get $A=PD{P}^{-1}$. Therefore, for any $\alpha \in \mathbb{R}$ we can write $A^{\alpha }=P{D}^{\alpha }{P}^{-1}$ or

$$\begin{array}{c}{A}^{\alpha }=\left(\begin{array}{cccc}{p}_{11}& p_{12}& \dots & p_{1n}\\ 0& p_{22}& \dots & p_{2n}\\ ⋮& ⋮& \ddots & ⋮\\ 0& 0& \dots & p_{nn}\end{array}\right)=\hfill \\ =\left(\begin{array}{cccc}{b}_{11}& b_{12}& \dots & b_{1n}\\ 0& b_{22}& \dots & b_{2n}\\ ⋮& ⋮& \ddots & ⋮\\ 0& 0& \dots & b_{nn}\end{array}\right)·\left(\begin{array}{cccc}{a}_1^{\alpha }& 0& \dots & 0\\ 0& a_2^{\alpha }& \dots & 0\\ ⋮& ⋮& \ddots & ⋮\\ 0& 0& \dots & a_n^{\alpha }\end{array}\right)·\left(\begin{array}{cccc}{c}_{11}& c_{12}& \dots & c_{1n}\\ 0& c_{22}& \dots & c_{2n}\\ ⋮& ⋮& \ddots & ⋮\\ 0& 0& \dots & c_{nn}\end{array}\right)=\\ =\left(\begin{array}{cccc}{b}_{11}& b_{12}& \dots & b_{1n}\\ 0& b_{22}& \dots & b_{2n}\\ ⋮& ⋮& \ddots & ⋮\\ 0& 0& \dots & b_{nn}\end{array}\right)·\left(\begin{array}{cccc}{a}_1^{\alpha }{c}_{11}& a_1^{\alpha }{c}_{12}& \dots & a_1^{\alpha }{c}_{1n}\\ 0& a_2^{\alpha }{c}_{22}& \dots & a_2^{\alpha }{c}_{2n}\\ ⋮& ⋮& \ddots & ⋮\\ 0& 0& \dots & a_n^{\alpha }{c}_{nn}\end{array}\right)=\\ \hfill =\left(\begin{array}{cccc}{a}_1^{\alpha }{b}_{11}{c}_{11}& a_1^{\alpha }{b}_{11}{c}_{12}+a_2^{\alpha }{b}_{12}{c}_{22}& \dots & a_1^{\alpha }{b}_{11}{c}_{1n}+a_2^{\alpha }{b}_{12}{c}_{2n}+...+a_n^{\alpha }{b}_{1n}{c}_{nn}\\ 0& a_2^{\alpha }{b}_{22}{c}_{22}& \dots & a_2^{\alpha }{b}_{22}{c}_{2n}+a_3^{\alpha }{b}_{23}{c}_{3n}+...+a_n^{\alpha }{b}_{2n}{c}_{nn}\\ ⋮& ⋮& \ddots & ⋮\\ 0& 0& \dots & a_n^{\alpha }{b}_{nn}{c}_{nn}\end{array}\right).\end{array}$$

Since $c_{mm}={\displaystyle \frac{1}{b_{mm}}}=1$, $m=1,...,n$ on the main diagonal we obtain $p_{mm}=a_m^{\alpha }$, $m=1,...,n$. For $s<m$ we get $p_{sm}=\sum _{j=s}^m{a}_j^{\alpha }{b}_{sj}{c}_{jm}$. For $b_{sj}$ and $c_{jm}$ we have representations by Theorem 3.1

$$b_{sj}=\left\{\begin{array}{cc}1\hfill & ifs=j;\hfill \\ 0\hfill & ifj<s;\hfill \\ {\displaystyle \frac{a_j^{j-s}}{\prod _{i=s}^{j-1}(a_j-a_i)}}\hfill & ifs<j,\hfill \end{array}\right.c_{jm}=\left\{\begin{array}{cc}1\hfill & ifj=m;\hfill \\ 0\hfill & ifm<j;\hfill \\ {\displaystyle \frac{a_j^{m-j-1}{a}_m}{\prod _{i=j+1}^m(a_j-a_i)}}\hfill & ifj<m.\hfill \end{array}\right.$$

When $j=s<m$ we get

$$b_{ss}=1,c_{sm}={\displaystyle \frac{a_s^{m-s-1}{a}_m}{\prod _{i=s+1}^m(a_s-a_i)}}={\displaystyle \frac{a_s^{m-s-1}{a}_m}{\prod _{i=s,i\ne s}^m(a_s-a_i)}},$$

when $j=m>s$ we get

$$c_{mm}=1,b_{sm}={\displaystyle \frac{a_m^{m-s}}{\prod _{i=s}^{m-1}(a_m-a_i)}}={\displaystyle \frac{a_m^{m-s}}{\prod _{i=s,i\ne m}^m(a_m-a_i)}}.$$

For $s<j<m$, $s<m$ we have

$$b_{sj}{c}_{jm}={\displaystyle \frac{a_j^{j-s}}{\prod _{i=s}^{j-1}(a_j-a_i)}}·{\displaystyle \frac{a_j^{m-j-1}{a}_m}{\prod _{i=j+1}^m(a_j-a_i)}}={\displaystyle \frac{a_j^{m-s-1}{a}_m}{\prod _{i=s,j\ne i}^m(a_j-a_i)}},$$

therefore for $s<m$

$$p_{sm}=\sum _{i=s}^m{a}_i^{\alpha }{\displaystyle \frac{a_i^{m-s-1}{a}_m}{\prod _{i=s,j\ne i}^m(a_j-a_i)}}.$$

Taking into account the form of $A^{\alpha }$ we get that $p_{sm}$ can be written as (25).

3.3. Resolvent of a Diagonalizable Triangle Matrix

The resolvent serves as a foundational tool for our considerations. Here, we derive the resolvent of a diagonalizable triangular matrix A. In this case, the resolvent takes on a particularly simple and computationally efficient form.

Lemma 3.1.

Let $\det A\ne 0$, $a_1\ne a_2\ne ...\ne a_n$, then the matrix

$$A=\left(\begin{array}{cccc}{a}_1& a_2& \dots & a_n\\ 0& a_2& \dots & a_n\\ ⋮& ⋮& \ddots & ⋮\\ 0& 0& \dots & a_n\end{array}\right)$$

has the resolvent im the form

$$R(\lambda ,A)={\left(r_{sm}\right)}_{m=1,s=1}^{n,n},$$

(26)

where

$$r_{sm}=\left\{\begin{array}{cc}0\hfill & ifm<s;\hfill \\ \sum _{i=s}^m{\displaystyle \frac{a_i^{m-s-1}{a}_m}{(\lambda -a_i)\prod _{j=s,j\ne i}^m(a_j-a_i)}}\hfill & ifs<m;\hfill \\ {\displaystyle \frac{1}{\lambda -a_m}}\hfill & ifm=s.\hfill \end{array}\right.$$

(27)

Proof. By Theorem 3.1 we obtain

$$\begin{array}{c}R(\lambda ,A)=PR(\lambda ,D)P^{-1}=\hfill \\ =\left(\begin{array}{cccc}{b}_{11}& b_{12}& \dots & b_{1n}\\ 0& b_{22}& \dots & b_{2n}\\ ⋮& ⋮& \ddots & ⋮\\ 0& 0& \dots & b_{nn}\end{array}\right)·\left(\begin{array}{cccc}{\displaystyle \frac{1}{\lambda -a_1}}& 0& \dots & 0\\ 0& {\displaystyle \frac{1}{\lambda -a_2}}& \dots & 0\\ ⋮& ⋮& \ddots & ⋮\\ 0& 0& \dots & {\displaystyle \frac{1}{\lambda -a_n}}\end{array}\right)·\left(\begin{array}{cccc}{c}_{11}& c_{12}& \dots & c_{1n}\\ 0& c_{22}& \dots & c_{2n}\\ ⋮& ⋮& \ddots & ⋮\\ 0& 0& \dots & c_{nn}\end{array}\right)=\\ =\left(\begin{array}{cccc}{b}_{11}& b_{12}& \dots & b_{1n}\\ 0& b_{22}& \dots & b_{2n}\\ ⋮& ⋮& \ddots & ⋮\\ 0& 0& \dots & b_{nn}\end{array}\right)·\left(\begin{array}{cccc}{\displaystyle \frac{c_{11}}{\lambda -a_1}}& {\displaystyle \frac{c_{12}}{\lambda -a_1}}& \dots & {\displaystyle \frac{c_{1n}}{\lambda -a_1}}\\ 0& {\displaystyle \frac{c_{22}}{\lambda -a_2}}& \dots & {\displaystyle \frac{c_{2n}}{\lambda -a_2}}\\ ⋮& ⋮& \ddots & ⋮\\ 0& 0& \dots & {\displaystyle \frac{c_{nn}}{\lambda -a_n}}\end{array}\right)=\\ \hfill =\left(\begin{array}{cccc}{\displaystyle \frac{b_{11}{c}_{11}}{\lambda -a_1}}& \sum _{i=1}^2{\displaystyle \frac{b_{1i}{c}_{i2}}{\lambda -a_i}}& \dots & \sum _{i=1}^n{\displaystyle \frac{b_{1i}{c}_{in}}{\lambda -a_i}}\\ 0& {\displaystyle \frac{b_{22}{c}_{22}}{\lambda -a_2}}& \dots & \sum _{i=2}^n{\displaystyle \frac{b_{2i}{c}_{in}}{\lambda -a_i}}\\ ⋮& ⋮& \ddots & ⋮\\ 0& 0& \dots & {\displaystyle \frac{b_{nn}{c}_{nn}}{\lambda -a_n}}\end{array}\right).\end{array}$$

Therefore

$$R(\lambda ,A)=\left(\begin{array}{cccc}{\displaystyle \frac{1}{\lambda -a_1}}& \sum _{i=1}^2{\displaystyle \frac{a_1}{(\lambda -a_i)\prod _{j=1,j\ne i}^2(a_i-a_j)}}& \dots & \sum _{i=1}^n{\displaystyle \frac{a_i^{n-2}{a}_n}{(\lambda -a_i)\prod _{j=1,j\ne i}^n(a_i-a_j)}}\\ 0& {\displaystyle \frac{1}{\lambda -a_2}}& \dots & \sum _{i=2}^n{\displaystyle \frac{a_i^{n-3}{a}_n}{(\lambda -a_i)\prod _{j=2,j\ne k}^n(a_i-a_j)}}\\ ⋮& ⋮& \ddots & ⋮\\ 0& 0& \dots & {\displaystyle \frac{1}{\lambda -a_n}}\end{array}\right)$$

that gives (26).

For constructing ${(-A)}^{\alpha }$, we assume now that $a_i>0$, $i=1,...,n$. Then since

$$\underset{0}{\overset{\infty }{\int }}{\displaystyle \frac{{\lambda }^{\alpha -1}}{\lambda +a_i}}d\lambda ={\displaystyle \frac{\pi a_i^{\alpha -1}}{sin\left(\alpha \pi \right)}},$$

by (1) we obtain (24).

3.4. Definition of Fractional Power of the Riemann–Stieltjes Integral Using Matrix Approach

It is clear that the constructions of the Section 3.1 are intended primarily to generalize the concept of a fractional integral of a function f with respect to another function g.

Let us consider a function $g\left(t\right)$ which is of bounded variation on $[a,b]$. It is known that in this case for $f\in C\left(\right[a,b\left]\right)$ there is exists the Riemann–Stieltjes integral

$$\underset{a}{\overset{x}{\int }}f\left(t\right)dg\left(t\right),x\le b.$$

If $g\left(t\right)\in C^1\left([a,b]\right)$, $w\left(x\right)=g^{\prime }\left(x\right)$ then

$$\left({\mathcal{I}}_{a+;w}^{1}f\right)\left(x\right)=\underset{a}{\overset{x}{\int }}f\left(t\right)dg\left(t\right)=\underset{a}{\overset{x}{\int }}f\left(t\right)g^{\prime }\left(t\right)dt=\underset{a}{\overset{x}{\int }}f\left(t\right)w\left(t\right)dt,x\le b.$$

Iterated operator $I_{a+,w}$ is

$$\left({\mathcal{I}}_{a+;w}^{k}f\right)\left(x\right)={\displaystyle \frac{1}{(k-1)!}}\underset{a}{\overset{x}{\int }}f\left(t\right){(g\left(x\right)-g\left(t\right))}^{k-1}w\left(t\right)dt.$$

(28)

Let $P=(x_0,x_1,\dots ,x_n)$ be a partition of $[a,x]$, $a<x\le b$ that is $a=x_0<x_1<\dots <x_n=x$, $x_k-x_{k-1}=h={\displaystyle \frac{x-a}{n}}>0$, so $x_0=a$, $x_1=x_0+h$, $x_2=x_0+2h$,..., $x_{n-1}=x_0+(n-1)h$, $x_n=x$, $f\left(a\right)=0$, $g\left(a\right)=0$.

The total variation of a continuous real-valued (or more generally complex-valued) function g is

$$V_a^b\left[g\right]=\underset{P}{sup}\sum _{i=0}^{n-1}|g\left(x_{i+1}\right)-g\left(x_i\right)|,$$

where the supremum is taken over the set of all partitions of the interval considered.

For the integral ${\mathcal{I}}_{a+;w}^1$ an estimate is valid

$$|\left({\mathcal{I}}_{a+;w}^{1}f\right)\left(x\right)|=\left|\underset{a}{\overset{x}{\int }}f\left(x\right)dg\left(x\right)\right|\le \underset{t\in (a,x]}{sup}\left|f\left(t\right)\right|V_a^x\left[g\right],$$

where $V_a^x\left[g\right]$ is the total variation of a function g on $[a,x]$.

Let ${\xi }_p\in [x-ph,x-(p+1)h]$, $w\left({\xi }_p\right)=w_p$, $f\left({\xi }_p\right)=f_p$, $p=0,...,n-1$, then we can write $I_{a+,w}$ in the form

$$\left({\mathcal{I}}_{a+;w}^{1}f\right)\left(x\right)=\underset{n\to \infty }{lim}\sum _{p=0}^{n-1}{f}_p{w}_{p}h=\underset{n\to \infty }{lim}\sum _{p=0}^{n-1}f\left({\xi }_p\right)w\left({\xi }_p\right)h.$$

We associate the operator $I_{a+;w}$ with the $(n\times n)$ matrix

$${\mathcal{A}}_{w;n}=h\left(\begin{array}{cccc}{w}_0& w_1& \dots & w_{n-1}\\ 0& w_1& \dots & w_{n-1}\\ ⋮& ⋮& \ddots & ⋮\\ 0& 0& \dots & w_{n-1}\end{array}\right)$$

(29)

and

$$\left({\mathcal{I}}_{a+;w}^{1}f\right)\left(x\right)=\underset{a}{\overset{x}{\int }}f\left(t\right)w\left(t\right)dt=\underset{n\to +\infty }{lim}({\mathcal{A}}_{w;n}·{\mathbf{f}}_n)·I_{1n},h={\displaystyle \frac{x-a}{n}},$$

(30)

where ${\mathbf{f}}_n^{\prime }=(f_0,f_1,...,f_{n-1})$, $I_{1n}={\underbrace{(1,0,...,0)}}_n$.

Taking ${\xi }_p=x-ph$, by the action of the approximating operator ${\mathbb{I}}_{n,a+;w}^1$ to f we will understand now

$$\left({\mathbb{I}}_{n,a+;w}^{1}f\right)\left(x\right)=\sum _{p=0}^{n-1}{f}_p{w}_{p}h=\sum _{p=0}^{n-1}{f}_p{w}_{p}h,f_p=(x-ph),w_p=w(x-ph).$$

Iterated operator ${\mathbb{I}}_{n,a+;w}^k$ has the form

$$\left({\mathbb{I}}_{n,a+;w}^{k}f\right)\left(x\right)={\displaystyle \frac{h^k}{(k-1)!}}\sum _{p=0}^{n-1}{\displaystyle \frac{(p+k-1)!}{p!}}f(x-ph)w(x-ph).$$

(31)

If $\det {\mathcal{A}}_{w;n}\ne 0$ we can construct the power ${\mathcal{A}}_{w;n}^{\alpha }=e^{\alpha ln{\mathcal{A}}_{w;n}}$, choosing for uniqueness the main branch of the logarithm and define ${\mathcal{A}}_{w;n}^{\alpha }$ by the Theorem 3.2. Therefore, we obtain the approximating operator ${\mathbb{I}}_{n,a+;w}^{\alpha }$ in the form

$$\left({\mathbb{I}}_{n,a+;w}^{\alpha }f\right)\left(x\right)=h^{\alpha }\sum _{p=0}^{n-1}\left(\sum _{i=0}^p{\displaystyle \frac{h^p{w}_i^{\alpha +p-1}{w}_p}{\prod _{j=0,j\ne i}^p(w_j-w_i)}}\right)·f_p.$$

(32)

Now we can apply Theorem 2.1 to operators $A={\mathcal{I}}_{a+;w}^1$ and $C={\mathbb{I}}_{n,a+;w}^1$. Let $B=C_w\left([a,b]\right)$ with the norm

$${\parallel f\parallel }_{C_w\left([a,b]\right)}=\underset{x\in [a,b]}{max}\left\{\right|f\left(x\right)w\left(x\right)\left|\right\}$$

and $D=C_w^1\left([a,b]\right)$ with norm

$${\parallel f\parallel }_{C_w^1\left([a,b]\right)}=\underset{x\in [a,b]}{max}{\left\{\right|f\left(x\right)w\left(x\right)|,|\left(f\left(x\right)w\left(x\right)\right)}^{\prime }\left|\right\},$$

$n\in \mathbb{N}$. It is known that for $f,w\in C\left(\right[a,b\left]\right)$, $\lambda >0$

$$\left|\right|\lambda R(\lambda ,{\mathcal{I}}_{a+;w}^1){\left|\right|}_{C_w\left([a,b]\right)}<M<\infty ,\left|\right|\lambda R(\lambda ,{\mathbb{I}}_{n,a+;w}^1){\left|\right|}_{C_w\left([a,b]\right)}<M<\infty .$$

and

$$\left|\right|\lambda R(\lambda ,{\mathbb{I}}_{n,a+;w}^1){\left|\right|}_{C_w^1\left([a,b]\right)}<M^{\prime }<\infty .$$

Also, it is easy to see that

$$\parallel ({\mathcal{I}}_{a+;w}^1-{\mathbb{I}}_{n,a+;w}^1){f\parallel }_{C_w\left([0,1]\right)}\le {\displaystyle \frac{(b-a)\mathbf{h}}{2}}{\parallel f\parallel }_{C_w^1\left([0,1]\right)},\mathbf{h}={\displaystyle \frac{b-a}{n}}.$$

Then

$$\parallel (R(\lambda ;{\mathcal{I}}_{a+;w}^1)-R(\lambda ;{\mathbb{I}}_{n,a+;w}^1)){f\parallel }_{C_w\left([0,1]\right)}\le {\displaystyle \frac{M{M}^{\prime }(b-a)\mathbf{h}}{2{\lambda }^2}}{\parallel f\parallel }_{C_w^1\left([0,1]\right)},$$

and

$$\parallel ({(-{\mathcal{I}}_{a+;w}^1)}^{\alpha }-{(-{\mathbb{I}}_{n,a+;w}^1)}^{\alpha }){f\parallel }_{C_w\left([0,1]\right)}\le \mathbf{c}{\parallel f\parallel }_{C_w^1\left([0,1]\right)}{\mathbf{h}}^{\alpha },\mathbf{h}\to 0,n\to \infty ,$$

$$\mathbf{c}=M\left(2+M^{\prime }{\displaystyle \frac{\alpha }{1-\alpha }}\right){\left({\displaystyle \frac{b-a}{2}}\right)}^{\alpha }.$$

3.5. Rate of Convergence for Concrete Fractional Integrals

1. Let us find the rate of convergence for the Riemann-Liouville integral. We consider $B=C\left(\right[0,1\left]\right)$, $D=C^1\left([0,1]\right)$ with norms

$${\parallel f\parallel }_{C\left(\right[0,1\left]\right)}=\underset{x\in [0,1]}{max}{\left\{\right|f\left(x\right)\left|\right\}\mathrm{and}\parallel f\parallel }_{C^1\left([0,1]\right)}=\underset{x\in [0,1]}{max}\left\{\right|f\left(x\right)|,|f^{\prime }\left(x\right)\left|\right\},$$

respectivy, $n\in \mathbb{N}$, $h={\displaystyle \frac{x}{n}}$

$$\left(Af\right)\left(x\right)=\left(I_{0+}^{1}f\right)\left(x\right)=\underset{0}{\overset{x}{\int }}f\left(t\right)dt,\left(Cf\right)\left(x\right)=\left({\mathbb{I}}_{n,0+}^{1}f\right)\left(x\right)=h\sum _{p=0}^{n-1}f\left(x-ph\right)$$

(see subSection 2.3). Then

$${\parallel (A-C)f\parallel }_{C\left(\right[0,1\left]\right)}=\parallel I_{0+}^{1}f-{\mathbb{I}}_{n,0+}^1{f\parallel }_{C\left(\right[0,1\left]\right)}\le {\displaystyle \frac{\mathbf{h}}{2}}{\parallel f\parallel }_{C^1\left([0,1]\right)},\mathbf{h}={\displaystyle \frac{1}{n}},$$

$$R(\lambda ,-I_{0+}^1)f={\displaystyle \frac{1}{\lambda }}f-{\displaystyle \frac{1}{{\lambda }^2}}\underset{0}{\overset{x}{\int }}{e}^{-{\displaystyle \frac{x-t}{\lambda }}}f\left(t\right)dt,$$

$$\begin{array}{c}R(\lambda ,-{\mathbb{I}}_{n,0+}^1)f={\displaystyle \frac{1}{\lambda }}f+\sum _{k=1}^{\infty }{\displaystyle \frac{{(-1)}^k}{{\lambda }^k}}{\displaystyle \frac{h^k}{\mathrm{\Gamma }\left(k\right)}}\sum _{p=0}^{n-1}{\displaystyle \frac{\mathrm{\Gamma }(k+p)}{\mathrm{\Gamma }(p+1)}}f\left(x-ph\right)=\hfill \\ ={\displaystyle \frac{1}{\lambda }}f+{\displaystyle \frac{1}{\lambda }}\sum _{p=0}^{n-1}{\displaystyle \frac{1}{\mathrm{\Gamma }(p+1)}}f\left(x-ph\right)\sum _{k=1}^{\infty }{\displaystyle \frac{{(-1)}^k\mathrm{\Gamma }(k+p)}{\mathrm{\Gamma }\left(k\right)}}{\left({\displaystyle \frac{h}{\lambda }}\right)}^k=\\ \hfill ={\displaystyle \frac{1}{\lambda }}f-{\displaystyle \frac{h}{\lambda }}\sum _{p=0}^{n-1}{\displaystyle \frac{{\lambda }^p}{{\left(\lambda +h\right)}^{p+1}}}f\left(x-ph\right).\end{array}$$

We get

$$\begin{array}{c}\parallel \lambda R(\lambda ,-I_{0+}^1){f\parallel }_{C\left(\right[0,1\left]\right)}\le \underset{x\in [0,1]}{max}\left\{\lambda \left|{\displaystyle \frac{1}{\lambda }}+{\displaystyle \frac{1-e^{-\frac{x}{\lambda }}}{\lambda }}\right|\right\}{\parallel f\parallel }_{C\left(\right[0,1\left]\right)}\le 2{\parallel f\parallel }_{C\left(\right[0,1\left]\right)},\hfill \\ \parallel \lambda R(\lambda ,-{\mathbb{I}}_{n,0+}^1){f\parallel }_{C\left(\right[0,1\left]\right)}\le \underset{x\in [0,1]}{max}\left\{\lambda \left|{\displaystyle \frac{1}{\lambda }}+{\displaystyle \frac{\left(1-{\lambda }^n{\left(\lambda +\frac{x}{n}\right)}^{-n}\right)}{\lambda }}\right|\right\}{\parallel f\parallel }_{C\left(\right[0,1\left]\right)}\le 2{\parallel f\parallel }_{C\left(\right[0,1\left]\right)},\\ \parallel \lambda R(\lambda ,-{\mathbb{I}}_{n,0+}^1){f\parallel }_{C^1\left([0,1]\right)}\le \\ \hfill \le \underset{x\in [0,1]}{max}\left\{\left|2-{\lambda }^n{\left(\lambda +{\displaystyle \frac{x}{n}}\right)}^{-n}\right|,\left|{\lambda }^n{\left(\lambda +{\displaystyle \frac{x}{n}}\right)}^{-n-1}\right|\right\}{\parallel f\parallel }_{C^1\left([0,1]\right)}=M^{\prime }{\parallel f\parallel }_{C^1\left([0,1]\right)},\end{array}$$

and all conditions of the Theorem 2.1 are satisfied, therefore

$$\parallel (R(\lambda ;-I_{0+}^1)-R(\lambda ;-{\mathbb{I}}_{n,0+}^1)){f\parallel }_{C\left(\right[0,1\left]\right)}\le {\displaystyle \frac{M^{\prime }\mathbf{h}}{{\lambda }^2}}{\parallel f\parallel }_{C^1\left([0,1]\right)},$$

and

$$\parallel I_{0+}^{\alpha }f-{\mathbb{I}}_{n,0+}^{\alpha }{f\parallel }_{C\left(\right[0,1\left]\right)}\le \gamma 2^{-\alpha }{\mathbf{h}}^{\alpha }{\parallel f\parallel }_{C^1\left([0,1]\right)},\mathbf{h}\to 0,n\to \infty .$$

2. Let us find the rate of convergence for the Erdélyi-Kober integral. We consider $B=C_x\left([0,1]\right)$, $D=C_x^1\left([0,1]\right)$, with norms

$${\parallel f\parallel }_{C_x\left([0,1]\right)}=\underset{x\in [0,1]}{max}{\left\{\right|xf\left(x\right)\left|\right\}\mathrm{and}\parallel f\parallel }_{C_x^1\left([0,1]\right)}=\underset{x\in [0,1]}{max}{\left\{\right|xf\left(x\right)|,|\left(xf\left(x\right)\right)}^{\prime }\left|\right\},$$

respectivy, $n\in \mathbb{N}$, $h={\displaystyle \frac{x}{n}}$,

$$\left(Af\right)\left(x\right)=\left(I_{0+;x}^{1}f\right)\left(x\right)=\underset{0}{\overset{x}{\int }}f\left(t\right)tdt,\left(Cf\right)\left(x\right)=\left({\mathbb{I}}_{n,0+;x}^{1}f\right)\left(x\right)=h\sum _{p=0}^{n-1}f(x-ph)·(x-ph)$$

(see subSection 2.4). Then

$${\parallel (A-C)f\parallel }_{C_x\left([0,1]\right)}=\parallel I_{0+;x}^{1}f-{\mathbb{I}}_{n,0+;x}^1{f\parallel }_{C\left(\right[0,1\left]\right)}\le {\displaystyle \frac{\mathbf{h}}{2}}{\parallel f\parallel }_{C_x^1\left([0,1]\right)},\mathbf{h}={\displaystyle \frac{1}{n}},$$

$$\parallel \lambda R(\lambda ,-I_{0+;x}^1){f\parallel }_{C_x\left([0,1]\right)}={∥f-{\displaystyle \frac{1}{\lambda }}\underset{0}{\overset{x}{\int }}{e}^{-{\displaystyle \frac{x^2-t^2}{2\lambda }}}f\left(t\right)tdt∥}_{C_x\left([0,1]\right)}\le 2{\parallel f\parallel }_{C_x\left([0,1]\right)},$$

$$\parallel \lambda R(\lambda ,-{\mathbb{I}}_{n,0+;x}^1){f\parallel }_{C_x\left([0,1]\right)}={∥f-\sum _{p=0}^{n-1}{\displaystyle \frac{{\lambda }^{p}f\left(x-ph\right)·(x-ph)}{{\left(\lambda +h\right)}^{p+1}}}∥}_{C_x\left([0,1]\right)}\le 2{\parallel f\parallel }_{C_x\left([0,1]\right)},$$

$$\parallel \lambda R(\lambda ,-{\mathbb{I}}_{n,0+;x}^1){f\parallel }_{C_x^1\left([0,1]\right)}\le M^{\prime }{\parallel f\parallel }_{C_x^1\left([0,1]\right)},$$

and all conditions of the Theorem 2.1 are satisfied, therefore

$$\parallel (R(\lambda ;-I_{0+;x}^1)-R(\lambda ;-{\mathbb{I}}_{n,0+;x}^1)){f\parallel }_{C_x\left([0,1]\right)}\le {\displaystyle \frac{M^{\prime }\mathbf{h}}{{\lambda }^2}}{\parallel f\parallel }_{C_x^1\left([0,1]\right)},$$

and

$$\parallel I_{0+;x}^{\alpha }f-{\mathbb{I}}_{n,0+;x}^{\alpha }{f\parallel }_{C_x\left([0,1]\right)}\le \gamma 2^{-\alpha }{\mathbf{h}}^{\alpha }{\parallel f\parallel }_{C_x^1\left([0,1]\right)},\mathbf{h}\to 0,n\to \infty .$$

The case of Hadamard-type fractional integral is considered similarly.

4. Discussion

The method presented in this article is universal, applicable to a broad class of operators and their fractional powers. A natural future application is to integrals with signed measures. The study of fractional powers of operators with alternating measures was initiated in the article [23].