Human basilar artery: morphology & variations

1.   Introduction

Fractional calculus (FC) is not a new research area; in reality, it has almost the same history as classical calculus. It motivates the study of derivatives and integrals of fractional order (FO). The history of FC was started in 1695 and appreciated during last few decades, when L'Hopital's question to Leibniz about the differentiation of $\dfrac{d^{\frac{1}{2}}}{dx}(x)$ ^[1]. For this, at that time, dated 30 September 1695, Leibniz responded to L'Hopital that "This is an apparent Paradox from which, one day, useful consequences will be drawn" and that was the birth of FC. After that, this concept of fractional differentiation, fractional integration, and its application has found in many numerous areas of the research field of sciences and engineering, especially in control engineering, electromagnetism, signal processing, fluid mechanics, diffusion process, biosciences, statistical, and continuum mechanics and many more. Also, from time to time, FC has been generalized by many researchers and mathematicians, namely Euler, Laplace, Lagrange, Fourier, Riemann, and many others. The differentiation of FO $\beta > 0$ has various definitions. But Riemann Liouville (RL) derivative and Caputo's derivative are the old and most commonly used. The CFD is used for our work as it has some advantage in dealing with the IVP of fractional differential equation of non-integer order $\beta > 0$. Also, many authors ^[2,3] have given the existence and uniqueness conditions of the IVP for FDEs.

In the study, it has been found that most of the IVP of FDEs don't have an exact scheme to find the solution, especially for non-linear FDEs. So, it becomes a challenging situation for researchers to establish some methods for finding the analytical solutions of FDEs. Therefore, many researchers have suggested several methods numerically for extended approximate solutions of integer differential equations into fractional differential equations. These schemes incorporates: fractional differential transform scheme, Adomain decomposition scheme ^[4], variational iteration method ^[5], spectral collocation scheme ^[6], fractional finite difference scheme ^[7], fractional Adams scheme ^[8], homotopy perturbation scheme ^[9], homotopy analysis scheme ^[10], extrapolation method ^[11], and many others.

In current decades, the IVP of FDEs used as a weapon to solve the various mathematical models, the epidemic model, the disease model, the dynamical system model, and many others. In recent times, Tong et al. ^[12] proposed fractional EM and fractional IEM, which are the generalization of classical EM and IEM for first-order IVP of FDEs. That EM has a linear convergence rate while IEM has a quadratic convergence rate. In ^[13] (Chapter 06), C. Milici et al. study several numerical methods for FO systems. In which, they proposed variational iteration, least squares, Euler's, and Runge-Kutta methods for the system of FDE in the CFD sense. In ^[14], Kumar et al. suggested a numerical scheme to demonstrate the numerical behavior of the IVP of FDEs in which one of the methods is midpoint point whose convergence rate is quadratic. In ^[15], Muhammad et al. developed a two-stage generalized Rk2 scheme of second order in the CFD sense. These all referred works motivate us to establish more accurate schemes to solve the IVP of FDE in the CFD sense. Also, our objective is to show, based on a few concrete examples and a few application models, that FDEs can model the physical problem more effectively than ODEs. Recently, many research article has been found that solve the real-world phenomenon in FDEs ^[16,17,18]. This work proposes a fractional RCM scheme for IVP of linear and non-linear FDEs of order $\beta \in [0, 1]$. This scheme has a cubic convergence rate. It is the generalization of classical RCM, developed by Anthony Ralston ^[19,20]. This proposed scheme has a more accurate approximation compared to existing fractional EM, IEM, Midpoint method (MPM), and many others for the IVP of FDE:

Here, $_{\mathcal C}\mathcal D^{\beta}_{x_0^+}$ indicates the CFD of arbitrary order $\beta$ where $\beta \in [0, 1]$.

By using the proposed work, our main goal is to establish a novel study which is more accurate and appropriate in order to derive the approximate solution of IVP of FDEs (1.1). This scheme incorporates the algorithm, order of convergence, stability, and few numerical examples including the application to WPG model. The scheme also recognize as a type of Runge-Kutta third order method (RK3) and this is explicitly familiar with the term \enquote{One-Step} method.

The presentation of the paper is designed as follows. In Section 2, we provide some preliminary definitions and properties of fractional derivatives and integrals. Next, in Section 3, we suggest our methodology briefly and its order of convergence by following some essential lemma and theorem. Then, in Section 4, we establish the stability of the concerned scheme. After that, in Section 5, we implemented the proposed method on a few examples of linear and nonlinear IVP of FDEs in the CFD frame. In Section 6, we solved World Population Growth (WPG) model via the suggested scheme with the comparison of EM and IEM. Finally, in Section 7, we conclude our methodology with some essential annotations.

2.   Preliminaries

This section focuses on some basic definitions of fractional derivatives and integrals, properties, and valuable results in the RL and Caputo derivative sense. This section will be helpful in our this advancement work as this will arise the generalization of ordinary calculus ^{[21,22,23,24,25]}.

Definition 2.1. ^[26] The FO integral in the sense of RL derivative for the function $\kappa:[a, b]\rightarrow \mathbb{R}$ of arbitrary order $\beta>0$ are

called the left and right RL fractional integral respectively. Here $\Gamma(\beta)$ denotes the Euler's Gamma function.

Definition 2.2. ^[26,27] The FO derivatives in the RL sense for the function $\kappa:[a, b]\rightarrow \mathbb{R}$ of order $\beta>0$ are

called the left and right RL fractional derivative respectively, where $n = 1+[{\beta}]$ and $[\beta]$ indicate the integral part of $\beta$. Particularly, if we take $0 < \beta < 1$, then

are called the left and right RL derivatives of order $\beta$, where $0 < \beta < 1$.

Definition 2.3. ^[26] The FO derivatives in the Caputo sense for the function $\kappa:[a, b]\rightarrow \mathbb{R}$ of order $\beta>0$ are

called the left and right Caputo derivative respectively, where $n = 1+[{\beta}]$.

Particularly, if we take $0 < \beta < 1$, then

are called the left and right FO Caputo derivatives of order $\beta$, where $0 < \beta < 1$.

The relation between the FO derivative in Caputo fractional derivative and Riemann-Liouville fractional derivative is

where $n = 1+[{\beta}]$.

Definition 2.4. ^[26] The one and two parameter Mittag-Leffler function are defined by,

respectively. If $\beta \in \mathbb{C}$ with $Re(\beta)>0$, then the series $E_\beta(\zeta)$ is convergent for all $\zeta\in \mathbb{C}$. Similarly, if $\beta, \gamma \in \mathbb{C}$ with $Re(\beta), Re(\gamma)>0$, then the series $E_{\beta, \gamma}(\zeta)$ is convergent for all $\zeta\in \mathbb{C}$.

Lemma 2.1. ^[26] If $\beta$, $\gamma \geq 0$, and $\Phi \in L_1\left[a, b\right]$, then

holds everywhere on the interval $\left[a, b\right]$. If $\Phi(x)\in C\left[a, b\right]$ or $1\leq\beta+\gamma$, then identity holds everywhere on the interval $\left[a, b\right]$.

Lemma 2.2. ^[28] If $\Phi\in C^n[a, b]$, $a < b$ and $n\in \mathbb{N}$. Moreover, If $\beta_1, \beta_2 > 0$ be such that, $\exists$ some $k\in \mathbb{N}$ with $k\leq n$ and $\beta_1$, $\beta_1+\beta_2 \in [k-1, k]$. Then,

Theorem 2.1. (Existence of IVP of FDE) ^[12] Let $g(x, u)$ be a function that hold the condition $g(x_0, u(x_0)) = 0$ and also the $g(x, u)$ is continuous on the domain $R:0\leq x-x_0\leq d$, $\left| u-u_0\right| \leq e$, then the FDEs:

has at least one solution in the interval $0\leq x-x_0\leq \lambda$ with $\lambda = \min\left\{d, \dfrac{e}{M}\right\}$ and $\max\limits_{(x, u)\in \mathbb{R}}{}_{\mathcal C}\mathcal D^{1-\beta}_{x_0^+}g(x, u) < M$.

Theorem 2.2. (Uniqueness of IVP of FDE) ^[12] Under the hypotheses of Theorem 2.1, and if $g_x(x, u)$ holds the Lipschitz condition in the variable $u$ with Lipschitz constant $0 < L$,

then the FDEs (2.1) have an unique solution.

3.   Proposed scheme

In our work, we are concerned about the approximate solution of the IVP for the linear and non-linear FDEs:

Here, we assume the derivative is in the CFD sense of FO $\beta$ where $\beta \in [0, 1]$. Now, with the help of Lemma 2.2, we apply suitable analogous transformation so that it becomes the classical differential equation. Then, we get the CFD of order $(1-\beta)$ and the revised IVP of FDEs:

Now, to discover the accurate numerical scheme to find the solution of the IVP of FDEs (3.1) is equivalent to locating the accurate numerical scheme for the IVP of FDEs (3.2). For such a precise solution of (3.2), we designed the RCM in fractional derivative operator, which we rename as fractional RCM, and it is more accurate and faster than all other linear and quadratic convergence rates like EM and IEM ^[12]. Below, we describe the fractional RCM.

Fractional RCM of FDE

For obtaining the approximate solution of FDE (3.2), we consider $(x_k, u_k)$ be the set points and we make this points accordingly so that the mesh are equally distribute in the whole interval $[a, b]$ where we set $x_0 = a$ and $x_{End} = b$. This idea will be good by selecting an integer which is non-negative say $N$ and assuming the mesh points. So, By following RCM, we construct the following algorithm:

Also with the use of Matlab, the RCM algorithm is proven to be an effective and much more accurate compared to linear and quadratic schemes.

Before showing the convergence of our suggested methodology, first we establish few relevant results which will be essential in the establishment of the convergence of the methods.

Lemma 3.1. ^[12] If $g_x(x, u)$ be a function that hold the condition of Lipschitz in the variable of $u$, for some Lipschitz constant $L>0$,

and also fulfil the conditions of Theorem 2.1, then $G(x, u) = {}_{{\mathcal C}}\mathcal D^{1-\beta}_{x_0^+}g(x, u)$ also holds the condition of Lipschitz in the of variable $u$, for some another Lipschitz constant $M>0$,

Lemma 3.2. If the function $G(x, y)$ follows the condition of Lipschitz for the variable $u$ and also satisfies the condition of Theorem 2.1, then

also satisfies the condition of Lipschitz in the variable of $u$.

Proof.

So, $\left|\tau(x, u_1)- \tau(x, u_2)\right|\leq L_\tau \left|u_1- u_2 \right|$, where $L_\tau = M\left(1+\frac{hM}{2!}+\frac{h^2M^2}{3!}\right).$

Theorem 3.1. Consider $g_x(x, u)$ be the function that holds the condition of Lipschitz in the variable of $u$ with Lipschitz constant $L>0$,

and $u(x)$ be the unique solution of IVP of FDEs (3.2).

Let $u_k$ be the generated solution approximation by Ralston's Cubic method for some non-negative integer $N$. Then for each $k = 0, 1, 2, \dots N$,

Proof. Let us take RCM iterative formula which is based on $u_k = u(x_k)$, then we get

Assuming, $G(x, u) = {}_{{\mathcal C}}\mathcal D^{1-\beta}_{x_0^+} g(x, u)$, then

By using the Taylor's series, the exact form of the solution will be:

Now, from the expression (3.4) and (3.5), we obtained $\left|u(x_{k+1})-\bar u_{k+1}\right| = O(h^4)$. So, we get

Let us assume that,

then using the stated Lemmas 3.1 and 3.2, we have

Now,

So, we get the estimation, $\left| E_{k+1}\right| = (1+hL_\tau)\left|E_k\right|+Kh^4.$

Thus, we get the recursion relation,

As, $x_k-x_0 = kh \; {\rm{and}}\; E_0 = 0\; {\rm{then}}, \; \; (1+hL_\tau)^k\leq e^{khL_\tau} = \phi_\tau.$

So, we have $\left|E_k \right|\leq \frac{Kh^3}{L_\tau}\; (\phi_\tau -1).$

Therefore, $\left|u(x_k)-u_k\right| = O(h^3).$

This conclude that RCM has a cubic convergence rate.

4.   Stability of proposed scheme

In this section, let us look at the numerical stability of our proposed scheme. Consider the IVP of FDE 3.2 in the simplest form by assuming $G(x, u) = {}_{{\mathcal C}}\mathcal D^{1-\beta}_{x_0^+}g(x, u)$ as,

The numerical solution of the RCM scheme is given by the formula:

In the simplest form, the concept of absolute stability ^[29,30] is based on the analysis of the behavior, according to the values of the step h, of the numerical solutions of the model equation:

The linearized equation uses $G(x, u) = \mu u$. In this case,

These combine to form,

Let us put $h\mu = z$, then the absolute stability region is the set

Let us examine the stability region of the numerical scheme RCM. Here, the stability region is the set of points such that $|\xi(h\mu)|\leq 1$, which is shown in the Figure 1. Note that for stability, the choice of $h$ must guarantee that $|h\mu|$ is inside the region. Furthermore, the real parts must be nonnegative for stability (or marginal stability). Therefore the regions to the left of the imaginary axis are the only ones of relevance.

5.   Numerical results and discussion

In this section, we illustrate how our proposed scheme operates in practice. We consider few examples of linear as well as non-linear IVP for FDEs in CFD form and solved numerically using the RCM scheme. First we compare the numerical solution with the analytical ones and then compare with the existing EM and IEM schemes. All the numerical computation have been carried out in MATLAB R2016a version. Now, before proceeding of numerical examples, we define few terminology.

The absolute error used in the table is defined as, $e_N = \max\limits_{j\in\mathbb{Z}}|u(x_j)-u_N(x_j)|.$ The used estimated oder of convergence (EOC) is defined by the quantity as,

where $u_N$ and $u_{2N}$ are the approximate solution at two distinct grids, with step length $h$ and $\frac{h}{2}$, respectively.

Numerical examples

In this subsection, we consider four examples of fractional IVPs of FDEs. In these four examples, the Examples 5.1 and 5.2 are linear, Examples 5.3 and 5.4 are non-linear. Here, we consider in Examples 5.1–5.3, the FDE has exact solution and Example 5.4, the FDE has no exact solution. These all FDEs are solved using the proposed fractional RCM including the proper comparison with the existing fractional EM and IEM.

Example 5.1. Consider the following fractional linear initial value problem (IVP) of FDE:

For $\beta = \frac{1}{2}$, the exact solution of (5.1) is,

Using our suggested scheme for $\beta = \frac{1}{2}$ and with step length $h = \frac{1}{10}$, the numerical solution of (5.1) are graphically represented in the Figure 2 and their exact solution, approximate solution are illustrated in the Table 1. Together with this, the absolute error visualization is indicated in Figure 3. The EOC and CPU performance of the schemes is tabulated in Table 2, and its plot is illustrated in Figure 4.

In the Figure 2, it represents the comparison between the exact and numerical result of our proposed scheme RCM with EM and IEM. Table 1 indicates the exact value, approximate value and absolute error of EM, IEM, and our proposed scheme RCM in which we observe that the numerical solution of RCM are more accurate to the exact solution. Figure 3 represents the absolute error of EM, IEM and our proposed scheme RCM, in which we notice that RCM has minimum absolute error in the comparison of EM and IEM, while IEM has second order of convergence. The order of convergence of our suggested scheme RCM, EM, IEM is tabulated in the Table 2 and graphically shown in the Figure 4. From there, it is clear that the order of EM is linear, the order of IEM is quadratic and the order of RCM is cubic. So, RCM is better than the EM and IEM. Next, in the Figure 4, the blue line indicates the linear convergence of EM, the magenta line indicates the quadratic convergence of IEM, and the green line indicates the cubic convergence of RCM respectively.

From our stated Example 5.1, the conclusion is that RCM is much more accurate in the comparison of the EM and IEM. Similar conclusion can also be drawn in the Examples 5.2 and 5.3.

Example 5.2. Consider the following linear IVP of FDE ^[31]:

For $\beta = \frac{1}{2}$, the exact solution of (5.2) is,

Using our suggested scheme for $\beta = \frac{1}{2}$ and with step length $h = \frac{1}{10}$, the numerical solution of (5.2) are graphically represented in the Figure 5 and their exact solution, approximate solution are illustrated in the Table 3. Together with this, the absolute error visualization is indicated in Figure 6. The EOC and CPU performance of the schemes is tabulated in Table 4, and its plot is illustrated in Figure 7.

Example 5.3. Consider the following non-linear IVP of FDE:

For $\beta = \frac{1}{2}$, the exact solution of (5.3) is,

Using our suggested scheme for $\beta = \frac{1}{2}$ and with step length $h = \frac{1}{10}$, the numerical solution of (5.3) are graphically represented in the Figure 8 and their exact solution, approximate solution are illustrated in the Table 5. Together with this, the absolute error visualization is indicated in Figure 9. The EOC and CPU performance of the schemes is tabulated in Table 6, and its plot is illustrated in Figure 10.

Example 5.4. Consider the following IVP of FDE:

With the help of our suggested scheme for $\beta = \frac{1}{10}$ and with step length $h = \frac{1}{10}$ and $h = \frac{1}{20}$, the numerical solution of (5.4) is graphically shown in the Figure 11 and their approximate solutions for $h = \frac{1}{20}$ is illustrated in the Table 7. In this example, the absolute error is calculated as the difference between the approximate solution at $N$ grid point and $2N$ grid point. The absolute error visualization is indicated in Figure 12 for $h = \frac{1}{20}$. The EOC and CPU performance of the schemes is tabulated in Table 8, and its plot is illustrated in Figure 13 for $h = \frac{1}{20}$.

6.   Application example and simulation

In this section, we consider one application of the real-world phenomenon WPG model in the form of CFD. We solve the model numerically using our proposed scheme. Also, we discuss the benefits of FC using the WPG model.

Example 6.1. Consider the following linear IVP of FDE of WPG model ^[32],

Here $\mathcal{N}(t)$ suggest the number of individuals population at any time $t$ and here our $P$ represents the production rate where $P = B-M$, $B$ indicates the rate of birth, and $M$ indicates the rate of mortality. Now, if we assume the value of $\beta = 1$ then in this case, our corresponding model will be linear population world growth model which is also famous as a classical population world growth model. The exact solution of (6.1) for $\beta = 1$ is,

Where $\mathcal{N}_0$ indicates the initial population at the initial time $t = t_0$. Our fractional model is better than the classical model from the numerical perspective. We have taken the population database from the year 1920 to 2018, that is around one century from the world population sites https://www.census.gov/data/tables/time-series/demo/international-programs/historical-est-worldpop.html or https://datacommons.org/place/Earth (provided by world bank), and also one is taken by United Nations ^[33]. These statistical population data match our fractional model scheme for $\beta = 1.393298754843208$. Also, The exact solution of (6.1) for the fractional model of population world growth model will be,

Now, in our classical population model scheme, the estimated value of our production rate is $P \approx 0.013501$ and for the fractional model scheme, the production rate $P \approx 0.0034399$ ^[34]. So, it has been found that the statistical population date value fit with the our fractional model for $\beta = 1.3932987548432$.

The world population data from the year 1920 to 2018 is graphically represented in the Figure 14 and the numerical value of (6.1) for $\beta = 1$ and $\beta = 1.393298754843$ with the step size $h = 1$ year is graphically shown in the Figure 15, tabulated in the Table 9. In addition to this, the absolute error visualization represented in Figure 16. From all the tables and figures, we conclude that our suggested scheme RCM is much more accurate and faster than EM and IEM.

From the year 1920 to June 2018, as per \enquote{The Census Bureau's International Data Base} indicates that the world population data reached around 7.5 billion. This data is very accurate and near to our fractional model for the FO $\beta = 1.3932987548432$. So, by our scheme, the numerical population data gives a more precise and appropriate solution than the fractional EM and fractional IEM.

7.   Conclusions

In this paper, the fractional RCM scheme is established for the IVP of FDE in CFD sense for the first time. Here, we do some analogous conversion of CFD of order $\beta$ into an ODE of integer order one, and then we operate our proposed scheme in the revised problem. The numerical scheme was used directly without consuming the perturbation, linearization, or other assumptions. The convergence analysis and stability analysis of the scheme has been proved. Also, in this work, we demonstrated a comparative numerical study of our proposed scheme with the comparison of the existing scheme fractional EM and fractional IEM for various examples of linear and non-linear FDEs. The scheme also solves one real-world phenomenon: The fractional WPG model. Now, here we conclude the significant benefits of our scheme:

● The fractional RCM scheme has a cubic convergence rate which is slightly faster than the other linear and quadratic convergence methods for the IVP of FDEs.

● The idea of computation of the proposed scheme is better, and we can get the desired approximation with the increment of mesh points.

● With the help of the WPG model, we discovered that FDEs fit the model better than ODE.

In the future, this work may be helpful to solve the IVP of FDE more accurately and effectively.

Acknowledgments

The first and third author extended their appreciation to Distinguished Scientist Fellowship Program (DSFP) at King Saud University, Saudi Arabia.

Conflict of interest

All authors declare no conflicts of interest.