The construction of solutions to ${}^{C} \mathit{\boldsymbol{{D}}}^{(1/n)}$ type FDEs via reduction to $\left({}^{C} \mathit{\boldsymbol{{D}}}^{(1/n)}\right)^n$ type FDEs

1.   Introduction

In the past few decades, fractional differential equations (FDEs) have gone from being a niche area of mathematical analysis to the forefront of mathematical modeling. Finding applications in a myriad of areas ranging from the classical FDEs in viscoelasticity ^[12], to more novel physical fields ^[9] and beyond to biology and medicine ^[13]. A review on more recent applications of fractional differential equations in a variety of research fields is given below.

One of the foremost fields of research to feature fractional derivatives in recent years is biomedicine. A type of fractional logistic differential equation used to model the COVID-19 pandemic is discussed in ^[4,17]. Continuous glucose monitoring is analyzed via a fractional differential equation model constructed from a noisy time series in ^[5]. Fractional differential equations have been used in a scheme to detect tea moisture content that was introduced in ^[27]. The memory property of fractional derivatives is exploited to study a combined drug treatment for the Human Immunodeficiency Virus (HIV) in ^[24]. The Gompertz law, used in many areas of biophysics, is generalized using fractional derivatives in ^[7]. An FDE model for the interaction of nutrient phytoplankton and its predator zooplankton is considered in ^[3].

Models of financial and economic processes have also recently featured fractional derivatives. A review of fractional differential equations used in economic growth models is given in ^[11]. Systems of FDEs are used in ^[30] to construct an indicator for the evaluation of economic development of a given region. The evolution of fractional-order chaotic financial systems is studied using the Adams-Bashforth-Moulton method in ^[28]. A financial crisis model represented by a system of fractional differential equations is analyzed in ^[18].

In physics and engineering, optics is a field where fractional differential equations find many uses. Semi-analytical solutions to the fractional Eikonal equation, a problem in optics, are constructed in ^[1]. The Caudrey-Dodd-Gibbon equation, used in laser optics, is analyzed in its fractional form in ^[23]. Optical soliton solutions to the conformable fractional Benjamin-Bona-Mahony equations are constructed in ^[31]. A fractional order model studying light distribution from the main fiber into other branch fibers in optical meta-materials is analyzed in ^[2].

Techniques for integrating FDEs can be classified into two large categories: numerical and analytical methods. Recently, there has been a surge of interest in numerical methods due to the increased reliance on FDE in fields of applied research. A review of classical methods is given in ^[6], while more recent algorithms are discussed in ^[14].

Analytical or semi-analytical techniques for the construction of solutions to FDEs have also experienced recent developments. The natural transform method was applied to contruct analytical solutions to a fractional oscillator in a resisting medium model in ^[10]. The Laplace-Adomian decomposition method is used to obtain the analytical solutions to a class of fractional-order dispersive partial differential equations in ^[20]. The same approach yields the solutions of fractional Zakharov-Kuznetsov equations in ^[21]. The q-homotopy analysis transform method is applied to solve a class of fractional diffusion equations in ^[22].

A particular class of techniques based on fractional power series has been presented in ^[15,16,26]. This approach considers the $\left(\, {}^{C} \mathit{\boldsymbol{{D}}}^{(1/n)}\right)^{n}$-type fractional equation:

where $\, {}^{C} \mathit{\boldsymbol{{D}}}^{(1/n)}$ denotes the Caputo derivative of order $\frac{1}{n}$ with respect to independent variable $x$; $F$ is an analytic function. Note that in the operator sense, the expression $\left(\, {}^{C} \mathit{\boldsymbol{{D}}}^{(1/n)}\right)^{n}$ is not equivalent to the integer-order derivative $\frac{{\rm{d}}}{{{\rm{d}}x}}$ – while the set of solutions to (1.1) does include solutions of the ordinary differential equation $y' = F(x, y)$, it is a much wider set ^[16].

It was demonstrated in ^[15] that (1.1) can be mapped to an equivalent ordinary differential equation (ODE) via the use of fractional power series. The solution to the obtained ODE can then be transformed into a solution to the original FDE (1.1). The main objective of this paper is to extend this approach to $\, {}^{C} \mathit{\boldsymbol{{D}}}^{(1/n)}$-type FDEs:

where $G(x, y)$ is an analytic function. It is demonstrated that FDE (1.2) can be transformed into (1.1) if specific conditions hold true, which can then be solved via the integration scheme presented in ^[25].

Note that while $\left(\, {}^{C} \mathit{\boldsymbol{{D}}}^{(1/n)}\right)^{n}$-type equations (1.1) do not necessarily have a physical interpretation, they are a vital part of the scheme presented in this paper for solving $\, {}^{C} \mathit{\boldsymbol{{D}}}^{(1/n)}$-type FDEs (1.2), which have a wide range of physical applications ^[9].

The paper is organized as follows: Section 2 contains preliminary results; Section 3 contains main definitions and derivations that demonstrate the transformation of $\, {}^{C} \mathit{\boldsymbol{{D}}}^{(1/n)}$-type FDEs into $\left(\, {}^{C} \mathit{\boldsymbol{{D}}}^{(1/n)}\right)^n$-type FDEs via the Riccati equation; Section 4 contains numerical experiments demonstrating the efficacy of the presented scheme.

2.   Preliminaries

2.1.   Fractional power series

In this paper, all functions $f(x)$ are represented via power series consisting of fractional-order powers of the independent variable. If a fractional derivative of order $\alpha = \frac{k}{n}; \gcd(k, n) = 1$ ($\gcd(k, n)$ denotes the greatest common divisor of integers $k$ and $n$) is considered, then the series parameter is set to $n$:

Series (2.1) is required to converge in the neighbourhood $0 \le x < R, \quad R > 0$. The series can be rewritten for a more convenient approach with regards to the Caputo fractional derivative in the following form:

where $w^{(n)}_j, j = 0, 1, \dots$ are the basis elements of series $f(x)$:

The following equality relates coefficients $c_j$ and $v_j$:

As mentioned previously, the series (2.2) and all subsequent fractional power series are required to converge in the neighbourhood $0 \le x < R, R > 0$.

Note that the substitution $t = x^{\frac{1}{n}}$ can be used to convert (2.1) (and (2.2)) into an integer-order power series $\widehat{f}(t)$:

The set of series given by (2.1) is denoted as ${}^C \mathbb{F}$. Multiplication between two elements $f, g \in {}^C \mathbb{F}$ is defined in the Cauchy sense:

since $w^{(n)}_j w^{(n)}_k = \binom{\left(j+k\right)/n}{k/n} w^{(n)}_{j+k}$ for any $j, k \in \mathbb{Z}_0$.

Note that the following property of the binomial coefficient is used in further sections:

where $\alpha, \beta \ge 0$.

More details on the properties of fractional power series are given in ^[15,16].

2.2.   Caputo fractional derivative operator

The Caputo fractional derivative will be considered in this paper. Let $\left(\, {}^{C} \mathit{\boldsymbol{{D}}}^{(1/n)}\right)^{n}$ denote the Caputo derivative of order $\frac{1}{n}$. The Caputo derivative acts on the basis elements (2.3):

The Caputo derivative of order $\alpha = \frac{k}{n}, \gcd(k, n) = 1$ is realized via taking the $k$th power the operator $\, {}^{C} \mathit{\boldsymbol{{D}}}^{(1/n)}$.

2.3.   The construction of analytical solutions to $\left(\, {}^{C} \mathit{\boldsymbol{{D}}}^{(1/n)}\right)^n$ type FDEs

A summary of the scheme for the construction of analytical solutions to $\left(\, {}^{C} \mathit{\boldsymbol{{D}}}^{(1/n)}\right)^{n}$ type FDEs is presented in this section. This scheme relies on the construction of an equivalent ODE via a characteristic function. The proof that a solution obtained using this scheme does satisfy the original FDE is given in ^[26].

Consider the following type $\left(\, {}^{C} \mathit{\boldsymbol{{D}}}^{(1/n)}\right)^n$ FDE:

where $F(x, y)$ is bivariate analytic function. The solution to (2.9) is constructed in the form of a fractional power series (as defined in Section 2.1):

Series (2.10) is convergent for $0 \le x < R, \quad R > 0$.

Inserting (2.10) into (2.9) yields the following relation:

Setting $t = x^{\frac{1}{n}}$ and rearranging (2.11) results in:

where $\widehat{y}$ is the integer power series that corresponds to the fractional power series (2.10):

Note that (2.12) is equivalent to the following ODE:

As shown in ^[26], inserting $t = x^{\frac{1}{n}}$ into the solution of the above equation yields a solution to the following Cauchy problem on (2.9):

The initial condition of fractional derivatives at $x = 0$ is due to (2.10) and the relation:

The algorithm for solving FDE (2.9) is depicted in Figure 1. Note that ^[25] outlines the algorithm for numerical integration of FDE (2.9) based on the extension of fractional power series via the use of generalized differential operators.

Example: The Riccati equation

Consider the following Cauchy problem on the Riccati fractional differential equation with constant coefficients:

where $\alpha, \beta \in \mathbb{R}$. The solution to (2.17) is a fractional power series (2.1) with $n = 2$:

The initial condition ${\, {}^{C} \mathit{\boldsymbol{{D}}}^{(1/2)} y}|_{x = 0} = \beta$ yields $c_1 = \frac{\beta}{\varGamma\left(\frac{3}{2}\right)}$. Furthermore, noting that:

and inserting the series (2.18) into (2.17) yields:

Using the substitution $t = \sqrt{x}$ and denoting $\widehat{y}(t) = y\left(t^2\right)$ transforms the above equation into:

where $P\left(\widehat{y}\right) = \widehat{y}^2 + \widehat{y} - 2$. Multiplying both sides of (2.21) by $2$ yields:

Rearranging the sum on the left-hand side of (2.22) and multiplying by $t$ results in:

Finally, adding $c_1 = \frac{\beta}{\varGamma\left(\frac{3}{2}\right)}$ to both sides results in the following ODE:

Note that the $\beta$, which is an initial condition to FDE (2.17) is a parameter in ODE (2.24).

The kink solitary solution to (2.24) is obtained for $\beta = 0$ in ^[29]. However, this case leads to coefficients $c_{2j+1} = 0, j = 0, 1, \dots$ which in its turn results in a solution to the ODE:

While the kink solitary solution does indeed satisfy (2.17), the entire set of solutions to the FDE is much wider. Every solution to (2.24) for some $\beta \in \mathbb{R}$ also satisfies {(2.17)} after the transformation $t = \sqrt{x}$.

For $\beta \neq 0$, ODE (2.24) can only be solved in series form, via expression of solutions by confluent hypergeometric series ^[19]. Solutions to both (2.24) and (2.17) are depicted in parts (a) and (b) of Figure 2 respectively. Note that the scale of the $x$-axis changes for the FDE and ODE respectively due to the substitution $t = \sqrt{x}$. This also shifts the singularity point from its position in Figure 2 (a) to that in Figure 2 (b).

3.   Main results

The main goal of the following derivations is to provide analytical techniques for the conversion of the type $\, {}^{C} \mathit{\boldsymbol{{D}}}^{(1/n)}$ problem into a problem of the type $\left(\, {}^{C} \mathit{\boldsymbol{{D}}}^{(1/n)}\right)^n$. Without the loss of generality and for the clarity of the presentation, the denominator of the fractional derivative order will be set to $n = 2$. The presented steps can be readily generalized for different values of $n$.

For clarity of presentation, subsequent sections discuss the application of the described scheme on the paradigmatic example of Riccati-type FDEs. However, the analytical and numerical computations can be performed for a general FDE of type (1.2).

3.1.   Auxiliary lemmas

In this section, three lemmas on the series solutions of the Riccati-type FDEs are given. The results presented here define auxiliary functions $\varTheta_z, \varPhi$ and $\varPsi$, which are essential for the transformation of the type $\, {}^{C} \mathit{\boldsymbol{{D}}}^{(1/n)}$ problem into a type $\left(\, {}^{C} \mathit{\boldsymbol{{D}}}^{(1/n)}\right)^n$ problem and, in its turn, for the construction of analytical solutions to (1.2).

Lemma 3.1. Let $z = \sum_{j = 0}^{+\infty} v_j w_j^{(2)} \in {}^C \mathbb{F}^{(1/2)}$ be any fractional power series. The Caputo derivative of $z^2$ reads:

where $\varTheta_z(x) = \sum_{j = 0}^{+\infty} \theta_j w_j^{(2)}$, $\theta_0 = 0$ and

Proof. Inserting the fractional power series $z = \sum_{j = 0}^{+\infty} v_j w_j^{(2)}$ into the left hand side of (3.1) yields:

Analogously, inserting $z = \sum_{j = 0}^{+\infty} v_j w_j^{(2)}$ and $\varTheta_z(x) = \sum_{j = 0}^{+\infty} \theta_j w_j^{(2)}$ into the right hand side of (3.1) results in:

Note that the function $\varTheta_z(x) = \, {}^{C} \mathit{\boldsymbol{{D}}}^{(1/2)} z^2 - 2z\, {}^{C} \mathit{\boldsymbol{{D}}}^{(1/2)} z$ quantifies the effect of fractional differentiation of $z^2$. If $\, {}^{C} \mathit{\boldsymbol{{D}}}^{(1/2)}$ is replaced by an integer order derivative $\frac{{\rm{d}}}{{{\rm{d}}x}}$, the function $\varTheta_z(x)$ becomes equal to zero.

The two following lemmas yield results on coefficients of the fractional power series solutions of the Riccati-type problems. Note that while the solution coefficients can be directly computed using these results, the evaluation of the solution does not readily follow (different numerical algorithms, such as described in ^[25] could be used for the evaluation).

Lemma 3.2. Consider the following Cauchy problem with respect to the Riccati fractional differential equation:

where $a_2, a_1, a_0, \gamma_0 \in \mathbb{R}$, and $\varPhi(x)$ is a given fractional power series with coefficients $\phi_j \in \mathbb{R}$:

The solution to (3.5) reads:

where,

Proof. Coefficients (3.8) are obtained by inserting the fractional power series (3.7) into (3.5).

Lemma 3.3. Consider the following Cauchy problem:

where $b_3, b_2, b_1, \lambda_0, \lambda_1 \in \mathbb{R}$ and $\varPsi(x)$ is a given fractional power series with coefficients $\psi_j \in \mathbb{R}$:

The solution to (3.9) reads:

where,

Proof. Coefficients (3.12) are obtained by inserting the fractional power series (3.11) into (3.9).

3.2.   The construction of solutions to the fractional Riccati equation

The results obtained in section 3.1 are now applied to derive the relationship between problems (3.5) and (3.9) as well as their respective solutions.

Consider the fractional Riccati equation (3.5). Differentiating (3.5) via the operator $\, {}^{C} \mathit{\boldsymbol{{D}}}^{(1/2)}$ yields:

Applying Lemma 1 to the first term of the right hand side of (3.13) yields:

Inserting $\, {}^{C} \mathit{\boldsymbol{{D}}}^{(1/2)} y_1 = a_2 y_1^2 + a_1 y_1 + a_0 + \varPhi(x)$ transforms (3.14) into:

Let us consider the following notation:

The function $\varPsi(x)$ is utilized in constructing solutions to FDE (3.5), while functions $\varPhi(x)$ and $\varTheta_{y_1}(x)$ are given in (3.5) and obtained from (3.1) respectively. Note that $\varPsi(x)$ simplifies to a linear function of $\varTheta_{y_1}(x)$ if $\varPhi(x) = 0$.

Comparing (3.15) and (3.9) yields the following relationship between coefficients $a_0, a_1, a_2$ and $b_0, \dots, b_3$:

Applying (3.16) and (3.17) transforms (3.15) into:

Note, that (3.18) has the same form as (3.9).

Moreover, (3.16) induces the following relations between the coefficients $\psi_j, \phi_j, \theta_j$ and $\gamma_j$, $j = 0, 1, \dots$:

The above derivations result in the following theorem.

Theorem 3.1. Cauchy problems (3.5) and (3.9) have the same solution $y_1 = y_2 = y = \sum_{j = 0}^{+\infty} \gamma_j w_j^{(2)}$ if relations (3.16), (3.17) and the following equalities:

do hold true.

Having derived the relationship between these problems, existing algorithms can be applied to solve (3.9), as detailed in Figure 3.

4.   Computational experiments

Consider the following Cauchy problem with respect to the Riccati fractional differential equation:

Note that in this case, the function $\varPhi(x) = 0$, thus, $\phi_j = 0; j = 0, 1, \dots$.

Using the Theorem 3.1, the values of the parameters $b_1$, $b_2$, $b_3$ and the initial conditions $\lambda_0$, $\lambda_1$ can be computed as follows:

Thus, (4.1) can be transformed into the following Cauchy problem:

where the coefficients of the function $\varPsi(x) = \sum_{j = 0}^{+\infty} \psi_j w_j^{(2)}$ are obtained using relations (3.19).

Following the technique outlined in Section 2.3, (4.3) can be converted into the following ODE:

where $t = \sqrt{x}$ and $\widehat{y} = \widehat{y}(t) = y(x)$. Note that the function $\varPsi(x)$ is changed into $\varPsi\left(t^2\right)$ due to the independent variable substitution. The function $\varPsi(t^2)$ can only be represented by an infinite power series (a known closed form of $\varPsi\left(t^2\right)$ does not exist). Thus, the above ODE cannot be solved directly. To integrate (4.4), $\varPsi(t^2)$ is approximated taking the first $N+1$ terms:

where $\widetilde{y}$ tends to $\widehat{y}$ as $N$ tends to infinity.

It is clear that the approximation of the series $\varPsi$ via the polynomial $\sum_{j = 0}^{N} \psi_j \frac{t^{2j}}{\varGamma\left(1+\frac{j}{n}\right)}$ introduces errors into the solution. Exact expressions and approximate numerical values of the coefficients $\psi_j \left(j = 0, 1, \dots, 8\right)$ are given in the Appendix A.

The solution to FDE (4.1) can now be obtained from the solution to the ODE (4.5) via the algorithm described in ^[25]. Figure 4 (a, b) depicts the solutions to (4.5) and (4.1) respectively for different values of $N$. These solutions are compared with a direct numerical solution computed via Garrappa's method ^[8] to (4.1) in Figure 4 (b). It can be seen from Figure 4 that increasing $N$ does cause the solution to converge, although that convergence is not monotonous.

In general, any numerical method can be used to construct solutions to (4.5). However, using the semi-analytical scheme presented in ^[25] makes it easier to perform the transformation of the time-axis, since the solution to (4.5) is given as a piecewise-polynomial function. If a purely numerical method is used (such as the classical Euler method, or any Runge-Kutta class technique), the nonlinear time axis transformation needs to be taken into consideration when selecting the numerical integration step-size.

Consider the following Cauchy problem with initial condition being equal to $\gamma_0$:

The root mean square error (RMSE) between solutions computed via the presented algorithm (denoted $y(x)$) and Garrappa's method (denoted $y_G(x)$) is defined as:

where $h$ denotes the integration step size; $M$ is the number of integration steps.

It can be seen in Figure 5 that for the initial condition $\gamma_0 \in [0.1, 0.3]$, RMSE between Garrappa's solution and solutions obtained by truncating $\varPsi(x)$ at $N = 1, \dots, 8$ significantly decreases up to $N = 4$. Using a higher-order approximation for $\varPsi(x)$ than $N = 4$ does not yield a significant improvement RMSE-wise.

5.   Conclusions

This paper proposes a new approach for solving $\, {}^{C} \mathit{\boldsymbol{{D}}}^{(1/n)}$-type FDEs. The construction of analytical solutions to a general form FDE without a direct evaluation of Caputo type integrals is a demanding mathematical problem. It has been demonstrated that some $\left(\, {}^{C} \mathit{\boldsymbol{{D}}}^{(1/n)}\right)^n$-type FDEs can be solved by transforming them into ODEs and applying a numerical algorithm ^[26].

The main contribution of this paper is the extension of the class of FDEs where similar fractional power series can be applied: The scheme is no longer limited to $\left(\, {}^{C} \mathit{\boldsymbol{{D}}}^{(1/n)}\right)^n$-type FDEs, but can be applied to $\, {}^{C} \mathit{\boldsymbol{{D}}}^{(1/n)}$-type FDEs. It opens new possibilities for the generation of solutions to such FDEs which previously could be analyzed using only approximate numerical techniques. Difficulties related to the application of the proposed technique are discussed in the paper and the presented numerical examples demonstrate the efficacy of the proposed technique.

Conflict of interest

The authors declare that they have no competing interests.

Appendix A

Exact expressions of the coefficients $\psi_j \left(j = 0, 1, \dots, 8\right)$ in (4.5) are as follows: