A third-order numerical method for solving fractional ordinary differential equations

1.   Introduction

Fractional calculus, as the name suggests, extends traditional integer-order calculus to fractional-order differential and integral calculus ^[1]. In recent decades, applications of fractional differential equations (FDEs) have received widespread attention in various disciplines, including physics ^[2], engineering ^[3] and biology ^[4]. Obtaining analytical solutions to FDEs is often challenging due to the nonlocal and complex nature of fractional derivatives. Therefore, it is crucial to develop efficient numerical methods for solving FDEs.

In the literature, two main numerical approximation strategies for solving fractional ordinary differential equations (FODEs) have been extensively studied. The first strategy involves a direct approximation of the fractional derivatives in the FODEs, while the second strategy focuses on solving the equivalent Volterra integral equations. For the first strategy, a second-order numerical method for solving FODEs was proposed in ^[5], where the fractional derivatives are approximated by a weighted sum of function values at specified points. A finite difference method of order $(2-\alpha)$ with L1, L2, and L2C formulas was established in ^[6], where $\alpha\in (0, 1)$ is a scalar representing the fractional order of the derivative. The L1-2 formula, an enhanced version of the L1 formula, was proposed in ^[7] to achieve a convergence rate of $3-\alpha$. Other approximation methods with the same convergence rate as obtained in ^[7] can be found in ^[8,9]. It is worth noting that the direct discretization method for fractional derivatives described in ^[8] is specifically designed for linear FODEs. Also, the convergence rates of most of these methods depend on the fractional order $\alpha$.

For the second strategy for solving FODEs, the classic $(1+\alpha)$th-order predicting-correction method was introduced in ^[10], and a detailed error analysis of this method was carried out in ^[11]. In ^[12], a method with the convergence order of $\min\{2, 1+2\alpha\}$ was developed. In ^[13], a one-step numerical integration method with second-order convergence rate, which is independent of $\alpha$, was proposed. The second method in ^[8] used piecewise quadratic interpolation polynomials to approximate the Volterra integral equation, achieving a convergence rate of $1+2\alpha$. In addition, the perturbed quadratic predictor-corrector (Q-PCP) and decomposed quadratic predictor-corrector (Q-PCD) methods were developed in ^[14], where the convergence rates of $3-\alpha$ and $3-\frac{\alpha}{2}$ were achieved, respectively. However, the convergence rates of almost all of the aforementioned methods depend on the fractional order $\alpha$.

In this paper, we present a third-order numerical scheme for the following general nonlinear FODEs:

where $\mathit{\boldsymbol{\alpha}} = (\alpha_{1}, \alpha_{2}, \ldots, \alpha_{n})^{\top}\in (0, 1]^{n}$ is a given vector of fractional orders; $\mathit{\boldsymbol{x}}(t) = (x_{1}(t), x_{2}(t), \ldots, x_{n}(t))^{\top}\in\mathbb{R}^{n}$ is the state vector at time $t$; $\mathit{\boldsymbol{f}}:\mathbb{R}\times\mathbb{R}^{n}\rightarrow \mathbb{R}^{n}$ is a given thrice continuously differentiable function with respect to all its arguments; $t_{0}$ is a given initial time; $t_{f}$ is a given terminal time; $\mathit{\boldsymbol{x}}^{0} = (x_{1}(t_{0}), x_{2}(t_{0}), \ldots, x_{n}(t_{0}))^{\top}\in \mathbb{R}^{n}$ is a given initial state vector; and ${_{t_{0}}^{C}}{D}_{t}^{{\mathbf{\alpha}}}{\mathit{\boldsymbol{x}}}(t) = \big({_{t_{0}}^{C}}{D}_{t}^{{\alpha_{1}}} {x_{1}}(t), {_{t_{0}}^{C}}{D}_{t}^{{\alpha_{2}}}{x_{2}}(t), \ldots, {_{t_{0}}^{C}}{D}_{t}^{{\alpha_{n}}}{x_{n}}(t) \big)^{\top}$ with ${_{t_{0}}^{C}}{D}_{t}^{{\alpha_{i}}}{x_{i}}(t)$ denoting the $\alpha_{i}$th Liouville-Caputo fractional derivative of $x_{i}(t)$.

We will develop a novel time-stepping numerical method to solve FODEs (1). Specifically, we first transform the FODEs into a set of equivalent Volterra integral equations. Then, we approximate these integral equations by using a third-order Taylor expansion (for the first two mesh points by a second-order Taylor expansion). This approximation leads to implicit nonlinear algebraic equations at each given mesh point. Then, Newton's method is employed to iteratively solve the nonlinear equations. A rigorous convergence analysis of the proposed method is carried out, which shows that the convergence rate from the numerical solution to the exact solution is third order and independent of the fractional order. Finally, the effectiveness and convergence results of the proposed method are illustrated through solving four numerical examples.

The organization of the rest of the paper is as follows. Section 2 proposes the numerical method for solving nonlinear FODEs. Section 3 provides the convergence and error analysis. Section 4 gives numerical results of four numerical examples. Finally, conclusions are drawn in Section 5.

2.   Numerical method

Applying the Riemann-Liouville integral to both sides of FODEs (1) yields the following Volterra integral equations ^[1]:

for $i = 1, 2, \ldots, n$, where $x_{i}^{0} = x_{i}(t_{0})$. For a given positive integer $N$, define a mesh on $[t_{0}, t_{f}]$ with $N+1$ mesh points $t_{q}, q = 0, 1, \ldots, N$, satisfying $t_{0} < t_{1} < \cdots < t_{N} = t_{f}$. For any $i = 1, 2, \ldots, n$, and $q = 0, 1, \ldots, N$, Eq (2) can be rewritten as

We now consider the approximation of the integrand $f_{i}(\tau, \mathit{\boldsymbol{x}}(\tau))$ on the RHS of Eq (3) in $(t_{0}, t_{1}]$ by the following Taylor expansion:

where $a_{i0}$ and $b_{i0}$ are coefficients to be determined, and $c_{i0}(\tau-t_{1})^{2}$ is the remainder. To determine $a_{i0}$ and $b_{i0}$, we omit the remainder and force the following equations to hold:

Solving these coupled equations gives

where $h_{1} = t_{1}-t_{0}$. Using these divided differences, Eq (4) can be expressed as

Substituting the RHS of Eq (5) into the first term of the sum in Eq (3) yields

where

For any $j = 1, 2, \ldots, q-1$, and $q = 2, 3, \ldots, N$, we approximate $f_{i}(\tau, \mathit{\boldsymbol{x}}(\tau))$ on $(t_{j}, t_{j+1}]$ by the following third-order Taylor expansion:

where $a_{ij}$, $b_{ij}$ and $c_{ij}$ are coefficients to be determined, and $d_{ij}(\tau-t_{j+1})^{3}$ is the remainder. Omitting the remainder, we can use the values of $f_{i}(\tau, \mathit{\boldsymbol{x}}(\tau))$ at points $t_{j-1}, t_{j}$ and $t_{j+1}$ to determine $a_{ij}, b_{ij}$ and $c_{ij}$, yielding the following divided differences:

where $h_{j+1} = t_{j+1}-t_{j}$; and $h_{j} = t_{j}-t_{j-1}$. Using these differences, we rewrite Eq (9) as

Substituting the RHS of Eq (10) into the $(j+1)$th term of the sum in Eq (3) gives

where

Combining Eqs (3), (6), and (11), we obtain the following equations, which are equivalent to the equations appearing in Eq (3),

for $i = 1, 2, \ldots, n,$ and $q = 1, 2, \ldots, N$, where $R_{i}^{q} = \sum\limits_{j = 0}^{q-1}R_{ij}$ is the cumulative truncation error up to the point $t_{q}$.

Omitting the truncation error $R_{i}^{q}$ in Eq (15), we define the following numerical scheme for Eq (3):

for $i = 1, 2, \ldots, n$, and $q = 1, 2, \ldots, N$, where $x_{i}^{q}$ represents an approximation of $x_{i}(t_{q})$ for each feasible $i$ and $q$.

We comment that Eq (16) defines an implicit time-stepping scheme for the numerical solution of Eq (1). Since Eq (16) is implicit, we need to use a technique, such as a Newton's method, to solve the nonlinear algebraic equations at each point $t_{q}, q = 1, 2, \ldots, N$.

3.   Convergence analysis

In this section, our focus is on the convergence analysis of the proposed numerical method.

Note that Eq (16) is a nonlinear system in $\mathit{\boldsymbol{x}}^{q}$ that will be solved by a Newton's iterative method. Thus, we need to show that the Jacobian matrix of this nonlinear system is invertible for any feasible $q$. For $i = 1, 2, \ldots, n$, and $q = 1, 2, \ldots, N,$ let

We then have the following theorem.

Theorem 3.1. Let $\mathit{\boldsymbol{F}}(\mathit{\boldsymbol{x}}^{q}) = (F_{1}(\mathit{\boldsymbol{x}}^{q}), F_{2}(\mathit{\boldsymbol{x}}^{q}), \ldots, F_{n}(\mathit{\boldsymbol{x}}^{q}))^{\top}$ for $q\in \{1, 2, \ldots, N\}$. Then, the Jacobian matrix of $\mathit{\boldsymbol{F}}(\mathit{\boldsymbol{x}}^{q})$ is invertible when $h_{q} > 0$ is sufficiently small.

Proof. From direct computation, we see that the Jacobian matrix of $\mathit{\boldsymbol{F}}(\mathit{\boldsymbol{x}}^{q})$, denoted as $\mathit{\boldsymbol{F}}^{\prime}(\mathit{\boldsymbol{x}}^{q})$, can be rewritten as

for $i = 1, 2, \ldots, n, l = 1, 2, \ldots, n$, and $q = 1, 2, \ldots, N$, where $\mathit{\boldsymbol{I}}_{n}$ is the $n \times n$ identity matrix; and $[b_{il}^{q}]_{n\times n}$ is an $n \times n$ matrix with

Since $\mathit{\boldsymbol{f}}$ is thrice continuously differentiable, $\frac{ {\partial{f}_{i}}}{ {\partial{x}_{l}}}$ is bounded on $[t_{0}, t_{f}]$ for $i = 1, 2, \ldots, n,$ and $l = 1, 2, \ldots, n$. Let

where $|\cdot|$ is the Euclidean norm.

Combining Eqs (7), (8), and (17), we obtain

for $i = 1, 2, \ldots, n,$ and $l = 1, 2, \ldots, n$. Furthermore, combining Eqs (12)–(14) and (17), we obtain

for $i = 1, 2, \ldots, n,$ $l = 1, 2, \ldots, n,$ and $q = 2, 3, \ldots, N$. Note that a key step in the derivation of Eq (20) is given in Appendix.

For given constants $\sigma_{i} \in (0, 1), i = 1, 2, \ldots, n$, define

Select $h_{1} \leq \hat{h}$ such that for $i = 1, 2, \ldots, n$,

Choose $h_{q}\leq \hat{h}$ such that for $i = 1, 2, \ldots, n$, and $q = 2, 3, \ldots, N$,

For $i = 1, 2, \ldots, n$, and $q = 1, 2, \ldots, N$, combining Eqs (22) and (23), we obtain

It follows from Eq (24) that $\mathit{\boldsymbol{F}}^{\prime}(\mathit{\boldsymbol{x}}^{q})$ is strictly diagonally dominant. According to the Levy-Desplanques theorem in ^[15], $\mathit{\boldsymbol{F}}^{\prime}(\mathit{\boldsymbol{x}}^{q})$ is non-singular. Therefore, $\mathit{\boldsymbol{F}}^{\prime}(\mathit{\boldsymbol{x}}^{q})$ is invertible, which completes the proof. □

Theorem 3.2. Let $\mathit{\boldsymbol{x}}(t_{q})$ be the exact solution of Eq (3), and let $\mathit{\boldsymbol{x}}^{q}$ be the solution of Eq (16) for $q = 1, 2, \ldots, N$. Then, there exists a positive constant $C$, independent of $\mathit{\boldsymbol{\alpha}}$, such that when $h_{1}$ is chosen to satisfy $h_{1}\leq h^{3/2}$, it holds that

where $\|\cdot\|_{\infty}$ is the infinity norm; $h = \max\limits_{\ell \in \{2, \ldots, N\}}\{h_{\ell}\}$ satisfies the condition $h\leq\hat{h}$ ($\hat{h}$ is as defined in Eq (21)).

Proof. The proof is carried out in three steps.

Step 1. The truncation error $\mathit{\boldsymbol{R}}^{q}$ is of order $\mathcal{O}(h^{3})$

Let $\mathit{\boldsymbol{R}}^{q} = (R_{1}^{q}, R_{2}^{q}, \ldots, R_{n}^{q})^{\top}$ for $q\in\{1, 2, \ldots, N\}$. From the definition of $R_{i}^{q}$ given by Eq (15) and the thrice continuous differentiability of $\mathit{\boldsymbol{f}}$, we deduce that for $i = 1, 2, \ldots, n$, and $q = 1, 2, \ldots, N$,

because $h_{1}\leq h^{3/2}$, where $d = \max\limits_{i\in\{1, \ldots, n\}\atop j\in\{1, \ldots, q-1\}}\{|c_{i0}|, |d_{ij}|\}$; and $h = \max\limits_{\ell \in \{2, \ldots, N\}}\{h_{\ell}\}$.

Recall that $\alpha_{i}\in (0, 1]$ for $i = 1, 2, \ldots, n$. Thus, $2^{\alpha_{i}-1}\leq \Gamma(\alpha_{i}+1)\leq 1$, as established in Theorem 4.1 of ^[16]. Also, it is obvious that $t_{f}^{\alpha_{i}} < \max \{1, t_{f}\}$. Therefore, choosing $\hat{C} = 2d\max\{1, t_{f}\}$, we obtain

for $q = 1, 2, \ldots, N$.

Step 2. The initial step of the mathematical induction process involving Eq (25)

We first show that Eq (25) holds for the case $q = 1$.

Subtracting Eq (16) from Eq (15) and applying the Lagrange mean value theorem, we obtain

for $i = 1, 2, \ldots, n$, and $l = 1, 2, \ldots, n$, where $u_{il}^{1}$ is defined as

Here, for $l = 1, 2, \ldots, n$, $\xi^{1}_{l} = x_{l}(t_{1})+\lambda_{l}^{1}(x_{l}^{1}-x_{l}(t_{1}))$ with $\lambda_{l}^{1}\in (0, 1)$.

Based on Eq (27), we obtain

for $i = 1, 2, \ldots, n$, and $l = 1, 2, \ldots, n$.

Similar to the derivation of Eq (19), it follows that $\big\|[u_{il}^{1}]_{1\times n}\big\|_{\infty}\leq \frac{ {h_{1}^{\alpha_{i}}M}}{ {\Gamma(\alpha_{i}+2)}}$, where $M$ is as defined in Eq (18). Thus, Eq (29) can be rewritten as

for all feasible $i$, which implies that

According to the derivation of Eq (22), if we select $h_{1}\leq\hat{h}$ (recall that $\hat{h}$ is as defined in Eq (21)), then $1-\frac{ {n h_{1}^{\alpha_{i}}M}}{ {\Gamma(\alpha_{i}+2)}} > \min\limits_{1\leq \iota \leq n}\{\sigma_{\iota}\} > 0$ is satisfied for all feasible $i$. Thus, $1-\max\limits_{1\leq i\leq n}\bigg\{\frac{ {n h_{1}^{\alpha_{i}}M}}{ {\Gamma(\alpha_{i}+2)}}\bigg\} > \min\limits_{1\leq \iota \leq n}\{\sigma_{\iota}\} > 0$ holds. From Eqs (26) and (30), we obtain

Step 3. Inductive step for Eq (25)

We now consider the case $q\geq 2$. In Step 2, we have shown that Eq (31), i.e., the case $q = 1$, is valid. To apply the mathematical induction, we assume that

for $j = 1, 2, \ldots, q-1$. Then, we need to show that $\|\mathit{\boldsymbol{x}}(t_{q})-\mathit{\boldsymbol{x}}^{q}\|_{\infty} = \mathcal{O}(h^{3})$.

Subtracting Eq (16) from Eq (15) and applying the Lagrange mean value theorem, we obtain

for $i = 1, 2, \ldots, n$, and $l = 1, 2, \ldots, n$, where $u_{il}^{1}$ is as defined in Eq (28), while $v_{il}^{j-1}$, $w_{il}^{j}$, and $z_{il}^{j+1}$ are defined as

Here, for $l = 1, 2, \ldots, n$, and $k = 0, 1, \ldots, q$, $\xi^{k}_{l} = x_{l}(t_{k})+\lambda_{l}^{k}(x_{l}^{k}-x_{l}(t_{k}))$ with $\lambda_{l}^{k}\in (0, 1)$.

Based on Eq (33), we can deduce that

for $i = 1, 2, \ldots, n$, and $l = 1, 2, \ldots, n$.

Similar to the derivation of Eq (20), we can show that $\big\|[z_{il}^{q}]_{1\times n}\big\|_{\infty} < \frac{ {4h_{q}^{\alpha_{i}}M}}{ {\Gamma(\alpha_{i}+1)}}$. Thus, Eq (34) can be expressed as

for all feasible $i$ and $l$, which implies that

According to the derivation of Eq (23), if we select $h_{q}\leq\hat{h}$, then $1-\frac{ {n4h_{q}^{\alpha_{i}}M}}{ {\Gamma(\alpha_{i}+1)}}\geq \min\limits_{1\leq \iota \leq n}\{\sigma_{\iota}\} > 0$ is satisfied for all feasible $i$. Thus, $1-\max\limits_{1\leq i\leq n}\bigg\{\frac{ {n4h_{q}^{\alpha_{i}}M}}{ {\Gamma(\alpha_{i}+1)}}\bigg\}\geq \min\limits_{1\leq \iota \leq n}\{\sigma_{\iota}\} > 0$ holds. From Eqs (26), (32), and (35) and the boundedness of $u_{il}^{1}$, $v_{il}^{j-1}$, $w_{il}^{j}$ and $z_{il}^{j+1}$ for $i = 1, 2, \ldots, n$, $l = 1, 2, \ldots, n$, and $j = 1, 2, \ldots, q-1$, we obtain

This completes the proof. □

Remark 1. The choice of $h_{1}\leq h^{3/2}$ aligns with the Rannacher time-stepping technique as discussed in ^[17,18]. This ensures that the overall accuracy of the proposed numerical method is not affected by the first time step.

Theorem 3.3. Let $\hat{\mathit{\boldsymbol{x}}}^{q}$ be the numerical solution of Eq (16) obtained by Newton's method for $q = 1, 2, \ldots, N$. Then, there exists a positive constant $\tilde{C}$, independent of both $h$ (as defined in Theorem 3.2) and $\mathit{\boldsymbol{\alpha}}$, such that

Proof. From Theorem 2.1 in ^[19], we note that there exists a well-defined number of iterations such that

where $h_{1}\leq h^{3/2}$ is as defined in Theorem 3.2. Then, Eq (37) can be expressed as

where $h$ is as defined in Theorem 3.2.

Combining Eqs (25) and (38), we obtain

Therefore, we can infer that Eq (36) holds. The proof is complete. □

Remark 2. The error estimate (36) shows that for $q = 1, 2, \ldots, N$, the numerical solution $\hat{\mathit{\boldsymbol{x}}}^{q}$ obtained by Newton's method converges to the exact solution $\mathit{\boldsymbol{x}}(t_{q})$ with third-order accuracy. Note that this convergence rate is independent of $\mathit{\boldsymbol{\alpha}}$, which is a significant distinction from the methods reported in ^[7,8,9,14].

4.   Numerical examples

In this section, we will solve four numerical examples to verify the effectiveness and convergence rate of the numerical method proposed in Section 3. Here, the stopping criterion for Newton's method is set to be $10^{-12}$, and all computations are conducted in the MATLAB 2022b environment on a PC equipped with a 2.80 GHz Intel Core i7-1165G7 CPU and 16.0 GB RAM.

For all test examples, the following graded mesh is used:

where $N$ is a positive integer. This graded mesh satisfies the condition in Theorem 3.2. The computed maximum error and convergence order are defined as

respectively, where $h$ is as defined in Theorem 3.2 under the partition number $N$; and $h^{\prime}$ is the maximum step size under the partition number $N^{\prime}$.

Remark 3. Note that other strategies can also be used in our calculations, as long as they comply with Remark 1.

4.1.   Example 1

Consider the following nonlinear FODE ^[14]:

The exact solution of this nonlinear FODE is $x(t) = t^{8}-3t^{4+\frac{\alpha}{2}}$. Our proposed numerical method is applied to solve this example. For comparison, we also solve this example by the one-step method in ^[13]. The computed numerical results and those reported in ^[14] are listed in Table 1, from which we see that the maximum errors in all methods decrease as the number of mesh points $N$ increases. Notably, the maximum error obtained by our method is significantly smaller than those obtained in ^[13] for different values of $\alpha$. Compared with ^[14], the convergence order computed by our method is comparable to that of Q-PCP for $\alpha = 0.1$ and significantly higher than those of Q-PCD for $\alpha = 0.8$ and $\alpha = 0.9$. This is because the theoretical convergence order of Q-PCP deteriorates as $\alpha$ increases. It is clear from Table 1 that the convergence order obtained by our method approaches to the theoretical value of 3. More importantly, the convergence order is independent of both $\alpha$ and $h$. Furthermore, we perform 50 tests and take the average time as the computational time for all test examples. From Table 1, we see that the one-step method, utilizing a second-order Taylor expansion and an explicit numerical format, achieves slightly shorter computational time compared to our method at the expense of lower accuracy. Figure 1 illustrates the state trajectories for various $\alpha$ when $N = 640$. From Figure 1, we can see that the numerical solutions obtained by our method closely match the exact solutions, which is consistent with the results in Table 1.

4.2.   Example 2

Consider the following nonlinear FODE ^[20]:

where $J_{0}(\cdot)$ is the zero-order Bessel function of the first kind ^[21].

This problem is solved again by using our proposed numerical scheme. For $\alpha = 0.5$, the exact solution of this example is $x(t) = \sin(4\sqrt{t})+0.01t^2+1$. For comparison, we also solve this example by using the implicit product integration trapezoidal rule (PI 2 Impl.) in ^[22]. The obtained numerical results are shown in Table 2. From Table 2, we see that the maximum error obtained by our method is much smaller than that obtained in ^[22] for various values of $N$. We also observe that the convergence order in ^[22] increases with the increasing value of $N$. In contrast, the convergence order computed by our method remains stable and is close to the theoretical value 3. This reconfirms that the convergence order of the proposed numerical method is independent of both $\alpha$ and $h$. Nevertheless, the computational time of the PI 2 Impl. method is significantly shorter compared to our method, attributed to a lower accuracy of the trapezoidal rule and the efficient treatment of persistent memory of fractional integrals. Figure 2 depicts the state trajectories when $N = 5120$. From Figure 2, it can be seen that the numerical solutions obtained by our method exhibit greater accuracy compared to those obtained using the method reported in ^[22].

4.3.   Example 3

Consider the following nonlinear FODE ^[9]:

The exact solution of this example is $x(t) = t^{3+\alpha}$. In contrast to the method proposed in ^[9] for directly approximating fractional derivatives and the method developed in ^[23], which employs the three-step Newton polynomial to approximate Volterra integrals, we use our proposed method and the MATLAB code provided in ^[23] to solve this problem for $\alpha = 0.3$, $\alpha = 0.6$, and $\alpha = 0.9$. The computed numerical results and those reported in ^[9] are listed in Table 3. From Table 3, we observe that the maximum error obtained by our method is significantly superior to those obtained in ^[9,23]. Furthermore, we note that our method exhibits a consistent convergence order close to 3 for various fractional orders, indicating its stability and efficiency. We observe that the convergence order obtained in ^[23] is also close to 3. However, the computational complexity in ^[23] is higher compared to our method. Specifically, for calculating $\mathit{\boldsymbol{x}}^{q}$, the values of $\mathit{\boldsymbol{x}}^{q-3}, \mathit{\boldsymbol{x}}^{q-2}$, and $\mathit{\boldsymbol{x}}^{q-1}$, $q = 3, 4, \ldots, N$, are required in ^[23], while the one-step Euler method and the two-step Adams-Bashforth method are used to obtain the values of $\mathit{\boldsymbol{x}}^{1}$ and $\mathit{\boldsymbol{x}}^2$, respectively. In contrast, our method only requires the values of $\mathit{\boldsymbol{x}}^{q-2}$ and $\mathit{\boldsymbol{x}}^{q-1}$ when calculating $\mathit{\boldsymbol{x}}^{q}$ for $q = 2, 3, \ldots, N$. The value of $\mathit{\boldsymbol{x}}^{1}$ can be calculated by using the second-order Taylor expansion. Thus, it can be seen from Table 3 that the computational time in ^[23] is higher than that in our method. The state trajectories for different fractional orders with $N = 2048$ are depicted in Figure 3. As it can be seen from Figure 3, our method shows satisfactory accuracy.

4.4.   Example 4

Consider the following nonlinear FODEs ^[24]:

The exact solutions of this example are $x_{1}(t) = e^{\frac{t}{2}}$ and $x_{2}(t) = te^{t}$ for $\alpha_{1} = \alpha_{2} = 1.0$. We solve this example for $({\alpha_{1}, \alpha_{2}})^{\top} = (0.7, 0.9)^{\top}$, $({\alpha_{1}, \alpha_{2}})^{\top} = (0.9, 0.7)^{\top}$, and $({\alpha_{1}, \alpha_{2}})^{\top} = (1.0, 1.0)^{\top}$ by our proposed method. Since the exact solutions for the cases $({\alpha_{1}, \alpha_{2}})^{\top} = (0.7, 0.9)^{\top}$ and $({\alpha_{1}, \alpha_{2}})^{\top} = (0.9, 0.7)^{\top}$ are unknown, we use the numerical solutions for $N = 2560$ as approximations to the exact solutions. The computed maximum error, convergence order, and computational time are presented in Table 4. It is observed from Table 4 that the value of the maximum error decreases as the value of $N$ increases. For different fractional orders, the computed convergence order of our method is roughly 3. To visualize the numerical results, we plot the state trajectories corresponding to different fractional orders in Figure 4.

5.   Conclusions

This paper has developed the third-order numerical method for solving nonlinear FODEs in the sense of Liouville-Caputo fractional derivatives. The fractional orders of these nonlinear FODEs can differ from each other. At each mesh point of a given mesh, we approximate the equivalent Volterra integral equations by using the third-order Taylor expansion (for the first subinterval, the second-order Taylor expansion is used). This approximation yields the implicit nonlinear algebraic equations that can be iteratively solved by the Newton's method. Furthermore, the convergence analysis and error estimate are performed, showing that the convergence rate of the proposed method is third order, independent of the fractional order. Finally, four non-trivial numerical examples are solved to illustrate the effectiveness and the convergence of the proposed method. In our future research, we will develop effective numerical methods for solving various fractional optimal control problems.

Author contributions

Xiaopeng Yi: Conceptualization, Investigation, Methodology, Software, Writing-original draft; Chongyang Liu: Conceptualization, Investigation, Methodology, Supervision, Software, Writing-review and editing; Huey Tyng Cheong: Conceptualization, Investigation, Methodology, Supervision, Writing-original draft; Kok Lay Teo: Conceptualization, Investigation, Methodology, Supervision, Writing-review and editing; Song Wang: Conceptualization, Investigation, Methodology, Writing-review and editing. All authors have read and agreed to the published version of the manuscript.

Use of AI tools declaration

The authors declare they have not used Artificial Intelligence (AI) tools in the creation of this article.

Acknowledgments

This work is supported by the Fundamental Research Grant Scheme of Malaysia (grant number FRGS/1/2021/STG06/SYUC/03/1), the National Natural Science Foundation of China (No. 12271307), and the Shandong Province Natural Science Foundation of China (No. ZR2023MA054).

Conflict of interest

Prof. Kok Lay Teo is an editorial board member for AIMS Mathematics and was not involved in the editorial review and/or the decision to publish this article.

The authors declare that they have no competing interests.

Appendix

In Eq (20), the following inequality holds: