Single-step and multi-step methods for Caputo fractional-order differential equations with arbitrary kernels

1.   Introduction

Fractional calculus has been as a mathematical theory of interest over three centuries. However, this theory was not initially applied to any real situation. The scientists and engineers in other fields commonly mentioned how the practical knowledge of fractional calculus has been used and how to operate it to the relevant studies. Fractional calculus and its application can be found in many fields such as physics ^[34,35], neural networks ^[32], mechanics and dynamic systems ^[12,38], biology ^[15,40] and economics ^[37].

The effect of increasing attention on fractional order models has made the investigation and growth of numerical methods for nonlinear fractional integro-differential equations ^{[5,6,7,8,9,10,11]} and partial differential equations with time or space fractional derivatives ^[13,31] become more active in the recent years. Moreover, the computational algorithm for the numerical solution of nonlinear fractional differential equations are found in ^{[17,18,19,23,24,25]}.

Researchers can obtain particular schemes for fractional differential equations based on the corresponding integral equation, for example, the product integration rules ^[39]. Nevertheless, the case of Adams predictor-corrector method ^[20] is different. It has been adjusted to create strong theoretical knowledge for numerical treatment of fractional differential equations. Some numerical schemes for fractional differential equations are developed directly from integral equations approaches, for example, product integration rules ^[39]. In addition, the Adams predictor-corrector approaches ^[20] is a completely different approach for the numerical solution of fractional differential equations.

Moreover, when the fractional order $\alpha$ converges to the nearest integer, the product integration rules and Adams predictor-corrector noticeably became the same methods. In fractional differential equations, the generalizations of the consistent approach for ordinary differential equations are considered as different methods.

The existence of different characteristics in approaches for fractional differential equations, which is established from the similar approach for ordinary differential equations, needs to be studied carefully. Particular distinctions are of theoretical concern, for instance, stability and convergence. However, these methods similarly indicate different types of computational nature, which have impact on the efficiency of the solution process. Therefore, our paper intends to explain various methods for fractional differential equations. A comparison of specific methods is also presented. The strengths and weaknesses will be investigated. Moreover, we compare the appropriateness among these methods.

Product integration rules and Adams predictor-corrector method are restricted to the analysis. These two approaches provide the least possible error constant ^[30], leading to a suitable balance between accuracy stability properties and computational complexity. This paper investigates when those strengths are rooted by the corresponding methods for fractional differential equations.

Additionally, we observe various definitions of fractional calculus ^[28,29] including Riemann-Liouville, Caputo, Hilfer, Riesz, Erdelyi-Kober and Hadamard, among others. In particular, R. Almeida ^[2] suggests some generalizations of fractional operators of a function with arbitrary kernels involving a weight function $\Phi$.

The $\Phi$-Caputo fractional differential equations have gained more attention recently. Numerical schemes for solving the $\Phi$-Caputo fractional differential equations are still under development. There are some research studies on the $\Phi$-Caputo fractional derivative. For example, Almeida et al. ^[4] indicated that mathematical models based on $\Phi$-Caputo fractional derivatives can be more adaptable. The $\Phi$-Caputo fractional derivative has the ability to extract hidden parts of real-world situations. In 2019, Almeida et al. ^[3] demonstrated $\Phi$-shifted Legendre polynomials to solve relaxation-oscillation equations with derivative of $\Phi$-Caputo. In ^[16,36], the monotone iteration of upper and lower solutions will be used to approximate the extremal solutions of $\Phi$-Caputo fractional differential equations.

To the finest of our understanding, numerical schemes for nonlinear fractional differential equations in the sense of $\Phi$-Caputo derivative have never been investigated. On the basis of the above works on fractional derivative with arbitrary kernels, this paper investigates the numerical approaches for the solution of nonlinear fractional differential equations involving $\Phi$-Caputo fractional derivative. In particular, the efficiency of numerical schemes, error and stability analysis are considered.

The paper is set as follows the next section, several preliminary knowledge of fractional derivatives and integrals are presented. The initial value problem involving $\Phi$-Caputo fractional derivative is defined in Section 3. Moreover, we presented four numerical schemes, namely explicit product integration rectangular rule, implicit product integration rectangular rule, implicit product integration trapezoidal rule, and Adams predictor-corrector method to find numerical solutions of the fractional differential equations in sense of the $\Phi$-Caputo fractional derivative. Next, the error estimation of the approximations and stability are obtained in Section 4. In Section 5, the simulation results including numerical convergence order are discussed for four test examples. The approximate solution and the error estimation for the test examples are presented through figures and tables, respectively. Further, a comparative study of these numerical schemes is performed.

2.   Preliminaries

In this section, we will examine basic definitions and theorems, which will be used to declare and verify our essential results. Let $\Phi$ be a continuously differentiable function on $[t_{0}, T]$ such that $\Phi^{\prime}(t) > 0,$ for all $t \in[t_{0}, T]$.

Definition 2.1 (^[14]). The Gamma function is defined as

Definition 2.2 (^[14]). The Euler beta function is defined as

Definition 2.3 (^[2]). The $\Phi$-Riemann-Liouville fractional integral of a real valued function $y$ on $[t_0, T]$ is given by

Definition 2.4 (^[2]). The $\Phi$-Riemann-Liouville fractional derivative of a real valued function $y$ on $[t_0, T]$ is defined by

where $n-1 < \alpha < n$ and $n\in \mathbb{N}$.

Definition 2.5 (^[2]). Let $n \in \mathbb{N}, \alpha \in (n-1, n)$ and $y, \Phi \in \mathscr{C}^{n}([t_0, T])$. The $\Phi$-Caputo fractional derivative of $y$ of order $\alpha$ is given by

Moreover, when $\alpha \in(0, 1)$, we have

In particular, the $\Phi$-Caputo fractional derivative is a generalization of the fractional derivatives as follows:

● the classical Caputo fractional derivative ^[29] when $\Phi(t) = t$,

● the Caputo-Erdélyi-Kober fractional derivative ^[33] when $\Phi(t) = t^{\sigma}$,

● the Caputo-Hadamard fractional derivative ^[21,27] when $\Phi(t) = \ln (t)$.

Theorem 2.1 (^[1]). If $y:[t_0, T] \rightarrow \mathbb{R}$ is a continuous function, then

3.   Numerical schemes

In this work, we study the $\Phi$-Caputo fractional derivative to differential equation as below:

where ${ }^{C} D_{t_0}^{\alpha, \Phi}$ is the $\Phi$-Caputo fractional derivative of order $\alpha \in (0, 1)$, $f : [t_{0}, T] \times \mathbb{R} \rightarrow \mathbb{R}$ is a given continuous function and $y_{0} \in \mathbb{R}$.

By Theorem 2.1, the solution of the problem (3.1) can be written in terms of the integral equation

However, the $\Phi$-Caputo fractional-order differential equation (3.1) is difficult to obtain the exact solution. In order to motivate the construction of our numerical methods, the concept of product integration rules ^[39] can be useed to estimate the integral in (3.2) by the appropriate polynomials.

In order to numerically solve the integral equation (3.2), we consider the approximation solutions $y_{n}, n = 1, 2, \dots, N$. The uniform grid is divided as

where $t_{j} = t_{0} +j h$ and $h = \frac{T-t_{0}}{N}$ for $0 \leq j \leq N$.

Additionally, we assume the approximations $y_{j} \approx y\left(t_{j}\right) (0 \leq j \leq n-1)$ in the basic step. The piecewise decomposition of (3.2) is defined as

3.1.   Explicit product integration rectangular rule (Ex. PI Rec.)

To obtain the approximation of $y_{n} \approx y\left(t_{n}\right)$, the function $f(\nu, y(\nu))$ in the integrand of (3.3) is approximated by the (explicit) forward Euler method, which is the constant values $f\left(t_{j}, y_{j}\right)$. The resulting methods is

where

3.2.   Implicit product integration rectangular rule (Im. PI Rec.)

We use a similar step of the explicit product integration rectangular rule to find the resulting method. However, the (implicit) backward Euler method is used to approximate, the function in the integrand of (3.3), which is

The resulting methods is

where $w^{\alpha, \Phi}_{n, j}$ is defined by (3.5).

3.3.   Implicit product integration trapezoidal rule (Im. PI Trap.)

In this method, we replace the function in each subinterval of the integrand (3.3) by the first-degree polynomial interpolant

The implicit product integration trapezoidal rule is

with $I_{j, j+1}^{\alpha, \Phi} = \int_{t_{j}}^{t_{j+1}}(\Phi(t_{n})-\Phi(\nu))^{\alpha} d\nu$, and

3.4.   Adams predictor-corrector method (Adams PC)

According to (3.7), Newton's method is necessary for approximating $y_{n}$. However, we can predict the approximation $y_{n}$ in (3.7) by using (3.4). The approximation in the predictor-step is called $y_{n}^{p}$. Furthermore, $y_{n}^{p}$ can be used in (3.7). This step is called corrector-step. Overall, this numerical scheme is called Adams predictor-corrector method and is defined by

where $y_{n}^{p}$ is the resulting method at predictor-step and $y_{n}$ is the resulting method at corrector-step.

In order to fulfill stability and error analysis, some lemmas are mentioned below.

Lemma 3.1. For $\alpha \in (0, 1)$ and $j = 1, \ldots, n-1$, we have

and

Proof. Applying the mean value theorem, we can find $\xi \in\left(\Phi(t_{n})-\Phi(t_{j+1}), \Phi(t_{n})-\Phi(t_{j})\right)$ such that

In a similar way, the second inequality can be proved.

Lemma 3.2. If $\alpha \in (0, 1), m \in \mathbb{N}$ and $\beta$ be nonnegative, then

Proof. Let $f(y) = \left(\Phi(t_{m})-\Phi(t_0)-y\right)^{\alpha-1} y^{\beta}$ for $0 < y < \Phi(t_{m})-\Phi(t_0)$.

Then, we have

This implies that $f(y)$ is a monotone increasing function, and then

which completes the proof.

In order to verify the stability and error analysis of our numerical schemes, we present a modified Gronwall's inequality as follows.

Lemma 3.3. Suppose $\alpha \in (0, 1), t_0 \in (0, T)$, and

Let $\sum\limits_{j = n}^{j = m} \tilde{b}^{\alpha, \Phi}_{j, m}\left|e_{j}\right| = 0$ for $1 \leq n < m$. Then, we have

if

where $M$ and $C$ are positive constants.

Proof. Let $c^{\alpha, \Phi}_{j, m} = M \tilde{b}^{\alpha, \Phi}_{j, m}$ be such that

From the inequality (3.9), we have $\left|e_{1}\right| \leq\left|\eta_{0}\right|$ to obtain

According to Lemma 3.2, it yields

Therefore,

Now, we want to show the following result:

In effect, we have

Applying the property of two gamma functions as follows:

Then, we have

Therefore, we obtain

which yields the convergence of $\sum\limits_{q = 1}^{\infty} \Phi_{q}$. Thus, the inequality of (3.10) is bounded, which means

The proof is completed.

Now, we outline the following assumptions on the nonlinearity $f$.

(H₁) The function $f: [t_{0}, T] \times \mathbb{R} \rightarrow \mathbb{R}$ is a continuous function.

(H₂) The function $f: [t_{0}, T] \times \mathbb{R} \rightarrow \mathbb{R}$ is Lipschitz continuous in the second variable, that is, there exists $L > 0$ such that

for all $t \in [t_{0}, T]$ and $y_{1}, y_{2} \in \mathbb{R}$.

(H₃) The function $f: [t_{0}, T] \times \mathbb{R} \rightarrow \mathbb{R}$ is Lipschitz condition in $t$ and $y(t)$ with a Lipschitz constant $K$.

Lemma 3.4. Assume that (H₁) holds. Let $y(t)$ be the solution of problem (3.2) and $h > 0$ be sufficiently small. Then, there exists a constant $C$ which is independent of $h$ such that

Proof. By the assumption (H₁), there is a positive constant $M$ such that $|f(t, y(t))| \leq M$. From the integral equation (3.2), we get

Then, we have

In the following, we claim that

In the fact that the inequality

one immediately gets

Combining (3.11) and (3.12), we obtain

which completes the proof.

Lemma 3.5. Assume that (H₁) and (H₃) hold. Let $y(t)$ be the solution of problem (3.2) and $h > 0$ be sufficiently small. Then, we have

where $\left(w^{\alpha, \Phi}_{n, j}- w^{\alpha, \Phi}_{n, j+1}\right)$ is given by (3.4).

Proof. By the assumption (H₃) and Lemma 3.4, which contributes to

Lemma 3.6. Assume that (H₁) and (H₃) hold. Let $y(t)$ be the solution of problem (3.2) and $h > 0$ be sufficiently small. Then, we have

where $\left(w^{\alpha, \Phi}_{n, j-1}- w^{\alpha, \Phi}_{n, j}\right)$ is given by (3.6).

Proof. The proof follows the same argument as in Lemma 3.5.

Lemma 3.7. If $f(t) \in \mathscr{C}^{1}[t_0, T]$, then

where $u^{\alpha, \Phi}_{n, 0}f(t_{0})+ \sum\limits_{j = 1}^{n} u^{\alpha, \Phi}_{n, j}f(t_{j})$ is given by (3.7).

Proof. Applying the mean value theorem, there exist $\xi_{j}, \eta_{j} \in\left(t_{j}, t_{j+1}\right)$ for $j = 0, 1, \ldots, n$ such that

This completes the proof.

Then, the stability analysis and error estimations of the explicit product integration rectangular rule (3.4), the implicit product integration rectangular rule (3.6), the implicit product integration trapezoidal rule (3.7), and the Adams predictor-corrector method (3.8) are investigated in the next section.

4.   Error estimation of the approximation and stability analysis

Theorem 4.1. Assume that (H₃) holds. Let $y(t) \in \mathscr{C}[t_0, T]$ be the solution of problem (3.2) and $y_{j}(0 \leq j \leq n-1)$ be the solution of the explicit product integration rectangular rule (3.4). Then, the error equation

holds.

Proof. For each $n = 1, 2, \ldots, N$, we let $e_{n-1} = y\left(t_{n-1}\right)-y_{n-1}$ and $e_{0} = 0$. By the integral equation (3.2) and the fractional left rectangle scheme (3.4), we have the error equation as follows:

By the assumption (H₃), and applying Lemmas 3.1 and 3.5, we obtain

It follows that the inequality (4.1) is obtained by using Lemma 3.3.

Theorem 4.2. Assume that (H₃) holds. Let $y(t)\in \mathscr{C}[t_{0}, T]$ be the solution of problem (3.2) and $y_{j}(0 \leq j \leq n-1)$ be the solution of the implicit product integration rectangular rule (3.6). Then, the error equation

holds.

Proof. The proof is essentially similar to the proof of Theorem 4.1.

Theorem 4.3. Assume that (H₁) and (H₃) hold. Let $y(t) \in \mathscr{C}[t_0, T]$ be the solution of problem (3.2) and $y_{j}(1 \leq j \leq N)$ be the solution of the implicit product integration trapezoidal rule (3.7). Then, the error equation

holds.

Proof. The proof of this theorem follows the same technique as Theorem 4.1.

Theorem 4.4. Assume that (H₁) and (H₂) hold. Let $y(t) \in \mathscr{C}[t_0, T]$ be the solution of problem (3.2), $D_{t_0}^{\alpha, \Phi} y(t) \in \mathscr{C}[t_0, T]$, and $y_{j}(1 \leq j \leq N)$ be the solution of the Adams predictor-corrector method (3.8). Then, the error equation

holds.

Proof. For each $n = 0, 1, \ldots, N-1$, we let $e_{n} = y\left(t_{n}\right)-y_{n}$ and $e_{0} = 0$. By the integral equation (3.2) and the Adams predictor-corrector method (3.8), we obtain

Then,

From Lemma 3.7, it follows that

By the assumption (H₂) of $f$, we obtain

As for $u^{\alpha, \Phi}_{n, n}$, one gets

Since $f(t, y(t)) = D_{t_0}^{\alpha, \Phi} y(t)$ is continuous and bounded, it follows Lemma 3.7 that

According to Lemmas 3.1 and 3.3, we obtain

which implies the statement of Theorem 4.4.

Theorem 4.5. Assume that (H₁) and (H₂) are true. Let $y_{j}(j = 1, 2, \ldots, n)$ be the solution to the explicit product integration rectangular rule (3.4) on the existed interval of its unique solution. Then the explicit product integration rectangular rule (3.4) is conditionally stable.

Proof. Suppose that $y_{0}$ and $y_{j} (j = 1, 2, \ldots, n)$ have perturbations $\tilde{y}_{0}$ and $\tilde{y}_{j}$, respectively. From (3.4), it follows

By the assumption (H₂), and the equations (3.4) and (4.3), we get

By Lemmas 3.1 and 3.3, it implies to

where $\eta_{0} = \max\limits_{0 \leq n \leq N-1}\left\{\left|\tilde{y}_{0}\right|+ \frac{L \left(w^{\alpha, \Phi}_{n, 0}-w^{\alpha, \Phi}_{n, 1}\right) \left|\tilde{y}_{0}\right|}{\Gamma(\alpha+1)} \right\}.$ The proof is completed.

In the same way, we obtain the stability results for the implicit rectangular and trapezoidal rules in Theorems 4.6 and 4.7, respectively. Hence, we omit the proof.

Theorem 4.6. Assume that (H₁) and (H₂) hold. Let $y_{j}(j = 1, 2, \ldots, n)$ be the solution to the implicit product integration rectangular rule (3.6) on the existed interval of its unique solution. Then the implicit product integration rectangular rule (3.6) is conditionally stable.

Theorem 4.7. Assume that (H₁) and (H₂) hold. Let $y_{j}(j = 1, 2, \ldots, n)$ be the solution to the implicit product integration trapezoidal rule (3.7) on the existed interval of its unique solution. Then the implicit product integration trapezoidal rule (3.7) is conditionally stable.

Next, we investigate the stability of the fractional predictor-corrector scheme (3.8).

Theorem 4.8. Assume that (H₁) and (H₂) hold. Let $y_{j}(j = 1, 2, \ldots, n)$ be the solution to the Adams predictor-corrector method (3.8) on the existed interval of its unique solution. Then, the Adams predictor-corrector method (3.8) is conditionally stable.

Proof. Assume that $\tilde{y}_{0}, \tilde{y}_{j}(j = 1, 2, \ldots, n)$, and $\tilde{y}_{n}^{p}(n = 1, \ldots, N)$ are perturbation terms of $y_{0}, y_{j}$, and $y_{n}^{p}$, respectively. Then, we construct the perturbation equations as follows:

From the assumption (H₂) and the inequality (4.4), we get

where

By Lemmas 3.1 and 3.3, it yields to

This completes the proof.

5.   Numerical examples

Motivated by ^[2], we assume that $J = [0, 1]$ and $\Phi \in \mathscr{C}(J)$ be a differentiable function such that $\Phi^{\prime}(t) > 0, t \in J$ and $\Phi(J) = [0, 1]$. Moreover, we solve the numerical examples by using MATLAB software and investigate different choices of suitable functions $\Phi$ in the numerical examples as below

Example 5.1. Consider the fractional order initial value problem given by

It is clearly seen that $y(t) = (\Phi(t))^{2}$ is the exact solution.

The exact solution and numerical solutions of (5.1) are plotted for different kernels $\Phi$ in Figure 1. Moreover, the numerical results are closed to the exact solution. Notice that the behaviors of the solutions are similar although we change the different kernels $\Phi$. In Table 1, we display the maximum errors of four numerical schemes for (5.1) when $\Phi(t) = t$ and $\alpha = 0.8$. From the data given in Table 1, the accuracy of numerical solutions corresponds the Theorems 4.1–4.4 when $h$ is sufficiently small. Tables 2 and 3 also present maximum errors for the kernels $\Phi(t) = \sin\left(\frac{t}{10}\right)$ and $\Phi(t) = \frac{t}{2}$, respectively. Overall, the implicit product integration trapezoidal rule and Adams predictor-corrector method have higher accuracy than the other method.

Example 5.2. Consider the nonlinear $\Phi$-Caputo fractional differential equations given by

From Figure 2, the numerical and exact solutions of (5.2) are plotted for kernels $\Phi$. Moreover, the numerical solutions are close to the exact solution. However, the behavior of the solutions is changed when the functions of $\Phi$ are change.

When $\alpha = 0.8$ and $\Phi(t) = t$, the maximum errors of four numerical schemes of (5.2) are presented in Table 4. The accuracy of numerical solutions also corresponds to the Theorems 4.1–4.4 when the step size $h$ is getting small. Next, the maximum errors of (5.2) with the different kernels $\Phi(t) = \sin\left(\frac{t}{10}\right)$ and $\Phi(t) = \frac{t}{2}$, respectively, for different values of $h$ are showed in Tables 5 and 6. Therefore, the implicit product integration trapezoidal rule has high accuracy than the other method. However, the implicit product integration rectangular rule gives higher accuracy than the other method when $\Phi(t) = \sin\left(\frac{t}{10}\right)$. Particularly, the four numerical schemes can be reduced to the numerical schemes in ^[22] when $\Phi(t) = t$. Furthermore, if $\Phi(t) = \ln (t)$, the numerical schemes are agreed in the example of ^[26]. It can be seen that $y(t) = (\Phi(t)-\Phi(t_0))^{2 \alpha}-(\Phi(t)-\Phi(t_0))^{2}$.

6.   Conclusions and future work

Four numerical schemes namely explicit product integration rectangular rule, implicit product integration rectangular rule, implicit product integration trapezoidal rule, and Adams method are extended to solve the nonlinear $\Phi$-Caputo fractional differential problems. We also analyze the error estimation and stability for those numerical schemes. When the exact solutions are known as in Examples 1 and 2, the implicit product integration trapezoidal rule and Adams predictor-corrector method can perform comparatively better than the explicit product integration rectangular and the implicit product integration rectangular rules, where the convergence order depends on the step size $h$. Moreover, these schemes are investigated for the linear and nonlinear $\Phi$-Caputo fractional differential equations. Further studies on numerical methods for fractional differential equations based on $\Phi$-Caputo derivative could be investigated for various types of fractional differential equations such as the fractional differential equations with delay and a nonlocal term, integro-differential equations, and higher order fractional differential equations.

Acknowledgments

The authors would like to thank the referees for their comments and suggestions which improve the quality of the manuscript. This research project is supported by Thailand Science Research and Innovation (TSRI), Basic Research Fund: Fiscal year 2022 under project number FRB650048/0164.

Conflict of interest

The authors declare no conflict of interest in this paper.