Successive midpoint method for fractional differential equations with nonlocal kernels: Error analysis, stability, and applications

Maroua Amel Boubekeur; Salah Boulaaras; Seda Igret Araz

doi:10.1515/phys-2025-0203

Article Open Access

Successive midpoint method for fractional differential equations with nonlocal kernels: Error analysis, stability, and applications

Maroua Amel Boubekeur , Salah Boulaaras and Seda Igret Araz

Published/Copyright: September 1, 2025

Published by

Become an author with De Gruyter Brill

Submit Manuscript Author Information Explore this Subject

From the journal Open Physics Volume 23 Issue 1

Abstract

This research introduces a numerical approach that employs the midpoint technique successively to solve differential equations that contain Atangana–Baleanu and Caputo fractional derivatives. The stability and consistency of the proposed approach has been demonstrated through the application of the technique suggested for solving the perturbed problem. Furthermore, the error analysis has been carried out for a nonlinear equation with Caputo derivatives. The effectiveness of this method is illustrated by examining various examples incorporating nonlocal kernels. Finally, by the use of an example, the well-known Fitzhugh–Nagumo model with Atangana–Baleanu fractional derivatives is simulated by using the proposed scheme. The method exhibits improved accuracy, consistency, and stability compared to the classical midpoint method, with particularly good performance for the Atangana–Baleanu and Caputo cases due to a specially designed predictor component. The novelty of the work lies in extending the method to fractional operators with non-local and non-singular kernels, offering a more robust approach than existing techniques in the literature.

Keywords: successive midpoint; stability; convergence; Fitzhugh–Nagumo model; fractional derivatives

1 Introduction and preliminaries

Numerical methods are mathematical techniques designed to find approximate solutions to complex problems that cannot be solved analytically. Given that real-world challenges in science, engineering, and finance often involve complicated equations, these methods serve as vital tools for modeling and addressing such problems. The evolution of numerical methods has significantly impacted various fields, enabling scientists and engineers to effectively simulate physical processes, analyze data, and optimize systems. This category of methods includes a diverse array of techniques such as interpolation, numerical integration, differentiation, and the resolution of ordinary and partial differential equations. Improvements in algorithms and computing technology have enhanced the accuracy and efficiency of numerical solutions, making them essential in both theoretical and practical contexts. Some of the methods that are fundamental to numerical analysis and are widely used in various scientific and engineering disciplines to solve practical problems are Euler methods [1], Runge–Kutta methods [2], midpoint methods [3], finite difference methods, interpolation techniques such as and Lagrange [4] and Newton interpolation [5–7], numerical integrations like the trapezoidal rule and Simpson’s rule, spectral methods [8], Kudryashov methods [9], and the collocation method [10].

Recently, a numerical technique called the successive midpoint method has been introduced in the literature [11]. While this method has been developed for classical, Caputo Fabrizio, global, and fractal derivatives, but its derivation for Atangana–Baleanu and Caputo fractional derivatives is not available in the literature [12–14]. This study develops a numerical method that utilizes the midpoint method iteratively to solve differential equations with Atangana–Baleanu and Caputo fractional derivatives. Then, we considered the Fitzhugh–Nagumo equation [15,16] which provides a simplified framework for the study of complex biological systems. This equation is extremely important for medical applications such as understanding the behavior of excitable cells, designing drugs that regulate heart rhythm, or studying neurological disorders such as epilepsy. By solving this equation, researchers can gain insight into emergent phenomena observed in biological networks. As a practical implementation of the proposed technique, we have established the scheme for the numerical solution of the Fitzhugh–Nagumo equation with Atangana–Baleanu fractional derivatives and carried out the numerical simulations for such model.

We shall present the definitions of the fractional derivatives that will be used here before proceeding with the derivation of the numerical method.

The definition of the Riemann–Liouville fractional derivative is represented by

(1) D t t 0 R L u ( t ) = 1 Γ ( 1 − α ) d d t ∫ t 0 t u ( τ ) ( t − τ ) − α d τ ,

where the function u ( ⋅ ) is continuous and 0 < α < 1 . The associated integral is given as

(2) I t t 0 R L u ( t ) = 1 Γ ( α ) ∫ t 0 t u ( τ ) ( t − τ ) α − 1 d τ .

The definition of the Caputo fractional derivative is defined by

(3) D t t 0 C u ( t ) = 1 Γ ( 1 − α ) ∫ t 0 t u ′ ( τ ) ( t − τ ) − α d τ ,

where the function u ( ⋅ ) is differentiable and 0 < α < 1 . The definition of the Atangana–Baleanu fractional derivative is given by the following:

(4) D t t 0 A B u ( t ) = 1 1 − α ∫ t 0 t u ′ ( τ ) E α − α 1 − α ( t − τ ) α d τ ,

where the function u ( ⋅ ) is differentiable and 0 < α < 1 . The associated integral is stated by

(5) I t t 0 A B u ( t ) = ( 1 − α ) u ( t ) + α Γ ( α ) ∫ t 0 t u ( τ ) ( t − τ ) α − 1 d τ .

The study is organized as follows: In Section 2, we present the derivation of a successive midpoint formula for a general Cauchy problem with the Atangana–Baleanu fractional derivative. In Section 3, we derive the successive midpoint formula for the Caputo case. Also, in these sections, the error of the exact and approximate solutions of the relevant problem is investigated, and the results are compared with the midpoint and L1 methods to demonstrate the accuracy of the proposed scheme. In Section 4, we demonstrate the stability and consistency of the associated method and present error analysis for a general initial value problem with classical and Caputo fractional derivatives. In Section 5, applications of the suggested method to the Fitzhugh–Nagumo model with Atangana–Baleanu fractional derivatives are presented and the numerical simulations are performed for such model.

2 Successive midpoint formula for a general Cauchy problem with Atangana–Baleanu fractional derivatives

In this section, we apply the successive midpoint formula introduced in the literature [11] to a general Cauchy problem with Atangana–Baleanu fractional derivatives

(6) D t α 0 A B u ( t ) = γ ( t , u ( t ) ) , u ( 0 ) = u 0 .

Applying the Atangana–Baleanu integral yields

(7) u ( t ) − u ( 0 ) = ( 1 − α ) γ ( t , u ( t ) ) + α Γ ( α ) ∫ 0 t γ ( τ , u ( τ ) ) ( t − τ ) α − 1 d τ .

At the point t = t n + 1 , we write

(8) u ( t n + 1 ) = u ( 0 ) + ( 1 − α ) γ ( t n + 1 , u ( t n + 1 ) ) + α Γ ( α ) ∫ 0 t n + 1 γ ( τ , u ( τ ) ) ( t n + 1 − τ ) α − 1 d τ .

Employing the midpoint method consecutively [11], we obtain

(9) u ( t n + 1 ) = u ( 0 ) + ( 1 − α ) γ ( t n + 1 , u ( t n + 1 ) ) + α Γ ( α ) ∑ j = 0 n γ t j + h 2 , u t j + h 2 × ∫ t j t j + 1 ( t n + 1 − τ ) α − 1 d τ = u ( 0 ) + ( 1 − α ) γ ( t n + 1 , u ( t n + 1 ) ) + h α Γ ( α ) ∑ j = 0 n γ t j + h 2 , u t j + h 2 × ( ( n − j + 1 ) α − ( n − j ) α ) .

Dividing the interval second time, we have

(10) u ( t n + 1 ) = u ( 0 ) + α Γ ( α ) ∑ j = 0 n ∫ t j t j + h 2 γ ( τ , u ( τ ) ) ( t n + 1 − τ ) α − 1 d τ + ( 1 − α ) γ ( t n + 1 , u ( t n + 1 ) ) + α Γ ( α ) ∑ j = 0 n ∫ t j + h 2 t j + 1 γ ( τ , u ( τ ) ) ( t n + 1 − τ ) α − 1 d τ = u ( 0 ) + ( 1 − α ) γ ( t n + 1 , u ( t n + 1 ) ) + h α Γ ( α ) ∑ j = 0 n γ t j + h 4 , u t j + h 4 × ( n − j + 1 ) α − n − j + 1 − h 2 α + h α Γ ( α ) ∑ j = 0 n γ t j + 3 h 4 , u t j + h 4 n − j + 1 − h 2 α − ( n − j ) α .

Dividing the third time, we obtain

(11) u ( t n + 1 ) = u ( 0 ) + ( 1 − α ) γ ( t n + 1 , u ( t n + 1 ) ) + α Γ ( α ) ∑ j = 0 n ∫ t j t j + h 4 γ ( τ , u ( τ ) ) ( t n + 1 − τ ) α − 1 d τ + α Γ ( α ) ∑ j = 0 n ∫ t j + h 4 t j + h 2 γ ( τ , u ( τ ) ) ( t n + 1 − τ ) α − 1 d τ + α Γ ( α ) ∑ j = 0 n ∫ t j + h 2 t j + 3 h 4 γ ( τ , u ( τ ) ) ( t n + 1 − τ ) α − 1 d τ + α Γ ( α ) ∑ j = 0 n ∫ t j + 3 h 4 t j + 1 γ ( τ , u ( τ ) ) ( t n + 1 − τ ) α − 1 d τ

= u ( 0 ) + ( 1 − α ) γ ( t n + 1 , u ( t n + 1 ) ) + h α Γ ( α ) ∑ j = 0 n γ t j + h 8 , u t j + h 8 × ( n − j + 1 ) α − n − j + 1 − h 4 α + h α Γ ( α ) ∑ j = 0 n γ t j + 3 h 8 , u t j + h 8 n − j + 1 − h 4 α − n − j + 1 − h 2 α + h α Γ ( α ) ∑ j = 0 n γ t j + 5 h 8 , u t j + 5 h 8 n − j + 1 − h 2 α − n − j + 1 − 3 h 4 α + h α Γ ( α ) ∑ j = 0 n γ t j + 7 h 8 , u t j + 7 h 8 n − j + 1 − 3 h 4 α − ( n − j ) α .

Using the idea of the midpoint consecutively and adding the predictor term, we obtain the following numerical scheme:

(12) u n + 1 = u ( 0 ) + ( 1 − α ) γ ( t n + 1 , u ˜ n + 1 ) + h α Γ ( α ) ∑ j = 0 n ∑ m = 0 2 k − 1 − 1 γ t j + ( 2 m + 1 ) h 2 k , u t j + ( 2 m + 1 ) h 2 k × n − j + 1 − m 2 k − 1 α − n − j + 1 − m + 1 2 k − 1 α .

It is worth noting that we need to calculate the predictor term. The predictor term is evaluated as

(13) u ˜ n + 1 = u ( 0 ) + ( 1 − α ) γ ( t n , u n ) + h α Γ ( α ) ∑ j = 0 n γ ( ( t j , u j ) ) ( ( n − j + 1 ) α − ( n − j ) α ) .

Example 1

We consider the following initial value problem:

(14) D t α 0 A B u ( t ) = t 2 u ( 0 ) = 0 ,

where the exact solution is

(15) u ( t ) = ( 1 − α ) t 2 + α Γ ( 3 ) Γ ( α + 3 ) t α + 2 .

Table 1 presents the numerical errors of the function u ( t ) obtained using the suggested method for various values of the fractional and classical orders. The results are compared with those of the midpoint (MM) [3] and L1 methods, where in the latter, the function γ ( t , u ( t ) ) is approximated using the first-order Lagrange interpolation [4]. Also, the CPU time (for an Intel Core 9, 2.7 GHz) for each simulation was recorded and presented in Table 1 to evaluate the computational efficiency of the proposed method.

Table 1

Error of the function u ( t ) for the suggested method

h	MM for α = 1 , k = 1	L1 method for α = 1	Suggested method for α = 1 , k = 10	CPU time (s)
0.1	8.3333 × 1 0 − 4	0.0078	3.1789 × 1 0 − 9	0.517
0.05	2.0833 × 1 0 − 4	0.0020	7.9473 × 1 0 − 10	0.956
0.01	8.3333 × 1 0 − 6	8.2833 × 1 0 − 5	3.1790 × 1 0 − 11	8.601

h	MM for α = 0.9	L1 method for α = 0.9	Suggested method for α = 0.9 , k = 11	CPU time (s)
0.1	0.0016	0.0265	3.8039 × 1 0 − 9	4.326
0.05	4.4283 × 1 0 − 4	0.0117	1.0326 × 1 0 − 9	16.841
0.01	2.2208 × 1 0 − 5	0.0021	4.9815 × 1 0 − 11	494.945

h	MM for α = 0.8	L1 method for α = 0.8	Suggested method for α = 0.8 , k = 12	CPU time (s)
0.1	0.0026	0.0450	3.4677 × 1 0 − 9	6.083
0.05	7.7191 × 1 0 − 4	0.0213	1.0022 × 1 0 − 9	22.375
0.01	4.5448 × 1 0 − 5	0.0041	5.5959 × 1 0 − 11	536.00

As seen in Table 1, the CPU time required by the suggested method increases significantly as the fractional order α decreases and the parameter k increases. This trend is particularly evident for smaller time step sizes h , where the CPU time grows rapidly – reaching over 536 s when α = 0.8 , k = 12 , and h = 0.01 . The reason for this behavior lies in the nonlocal nature of fractional derivatives: as α decreases, the memory effect becomes more dominant, requiring more historical data to be processed at each time step. Additionally, increasing k introduces further subdivisions and intermediate calculations, contributing to higher computational complexity. In contrast, the midpoint [3] and L1 methods [4] show relatively stable and much lower CPU times, even at smaller h values. However, this comes at the cost of reduced accuracy, as reflected in their higher error values compared to the suggested method. In Figure 1, the numerical simulation is provided for the above problem with the Atangana–Baleanu fractional derivative.

$Figure 1 Numerical simulation for the Cauchy problem with the Atangana–Baleanu fractional derivative for α = 0.8 , k = 11 \alpha =0.8,k=11 .$

Figure 1

Numerical simulation for the Cauchy problem with the Atangana–Baleanu fractional derivative for α = 0.8 , k = 11 .

3 Successive midpoint method for a general Cauchy problem with Caputo fractional derivatives

In this section, we apply the successive midpoint formula introduced in the literature [11] to a general Cauchy problem with Caputo fractional derivative

(16) D t α 0 C u ( t ) = γ ( t , u ( t ) ) , u ( 0 ) = u 0 .

Applying the associated integral yields

(17) u ( t ) − u ( 0 ) = 1 Γ ( α ) ∫ 0 t γ ( τ , u ( τ ) ) ( t − τ ) α − 1 d τ .

At the point t = t n + 1 , we write

(18) u ( t n + 1 ) = u 0 + 1 Γ ( α ) ∫ 0 t n + 1 γ ( τ , u ( τ ) ) ( t n + 1 − τ ) α − 1 d τ .

The above is converted to

(19) u ( t n + 1 ) = u 0 + 1 Γ ( α ) ∑ j = 0 n ∫ t j t j + 1 γ ( τ , u ( τ ) ) ( t n + 1 − τ ) α − 1 d τ .

Using the idea of the successive midpoint, we obtain the following numerical scheme:

(20) u n + 1 = u 0 + h α Γ ( α + 1 ) ∑ j = 0 n ∑ m = 0 2 k − 1 − 1 γ t j + ( 2 m + 1 ) h 2 k , u t j + ( 2 m + 1 ) h 2 k × n − j + 1 − m 2 k − 1 α − n − j + 1 − m + 1 2 k − 1 α .

Example 2

We consider the following initial value problem:

(21) D t α 0 C u ( t ) = t 3 + t 2 , u ( 0 ) = 0 ,

where the exact solution is

(22) u ( t ) = 6 Γ ( α + 4 ) t α + 3 + 2 Γ ( α + 3 ) t α + 2 .

In Table 2, we present the error of the function u ( t ) by employing the suggested method for different values of fractal order and classical order. The total CPU times of the methods were computed and reported in Table 2 in order to compare the computational performance.

Table 2

Error of the function u ( t ) for the suggested method

h	MM for α = 1 , k = 1	L1 method for α = 1	Suggested method for α = 1 , k = 8	CPU time (s)
0.1	0.0021	0.0189	1.2716 × 1 0 − 7	0.376
0.05	5.2083 × 1 0 − 4	0.0050	3.1789 × 1 0 − 8	0.452
0.01	2.0833 × 1 0 − 5	2.0634 × 1 0 − 4	1.2716 × 1 0 − 9	2.712

h	MM for α = 0.9	L1 method for α = 0.9	Suggested method for α = 0.9 , k = 9	CPU time (s)
0.1	0.0044	0.0207	1.4185 × 1 0 − 7	2.337
0.05	0.0012	0.0054	3.8621 × 1 0 − 8	7.374
0.01	6.0988 × 1 0 − 5	2.2211 × 1 0 − 4	1.8742 × 1 0 − 9	250.453

h	MM for α = 0.8	L1 method for α = 0.8	Suggested method for α = 0.8 , k = 10	CPU time (s)
0.1	0.0079	0.0228	7.9473 × 1 0 − 9	4.525
0.05	0.0024	0.0059	1.9868 × 1 0 − 9	14.288
0.01	1.4047 × 1 0 − 4	2.3858 × 1 0 − 4	7.9474 × 1 0 − 11	692.56

As observed from Table 2, the CPU time required for the proposed method increases significantly as the fractional order α decreases and the associated parameter k increases. This is particularly evident for smaller step sizes h where the computational burden becomes more pronounced. The increased in CPU time is due to the enhanced memory effect as the fractional order α decreases and the increased number of divisions resulting from larger values of k . In contrast, classical methods (midpoint and L1) remain relatively stable in terms of CPU time, but exhibit inferior accuracy compared to the suggested approach. In Figure 2, the numerical simulation is provided for the above problem with the Caputo fractional derivative.

$Figure 2 Numerical simulation for the Cauchy problem with the Caputo fractional derivative for α = 0.8 , k = 10 \alpha =0.8,k=10 .$

Figure 2

Numerical simulation for the Cauchy problem with the Caputo fractional derivative for α = 0.8 , k = 10 .

4 Consistency and stability analysis

Consistency of a method refers to its ability to produce results that converge to the exact solution of the ordinary differential equations as the step size (or time increment) approaches zero [5]. In other words, a consistent numerical method ensures that as you refine the discretization (by decreasing the step size), the numerical solution approaches the exact solution of the differential equation [15–24]. In this section, we present the stability and consistency of the suggested numerical scheme. We start with stability, we choose u ˜ n + 1 and u ˜ j + ( 2 m + 1 ) h 2 k the perturbed terms of u n + 1 and u j + ( 2 m + 1 ) h 2 k , respectively. We recall that the approximate solution is given as

(23) u n + 1 = u 0 + h α Γ ( α + 1 ) ∑ j = 0 n ∑ m = 0 2 k − 1 − 1 γ t j + ( 2 m + 1 ) h 2 k , u j + ( 2 m + 1 ) h 2 k × n − j + 1 − m 2 k − 1 α − n − j + 1 − m + 1 2 k − 1 α .

The perturbed equation is

(24) u n + 1 + u ˜ n + 1 = u 0 + u ˜ 0 + h α Γ ( α + 1 ) × ∑ j = 0 n ∑ m = 0 2 k − 1 − 1 γ t j + ( 2 m + 1 ) h 2 k , u j + ( 2 m + 1 ) h 2 k + u ˜ j + ( 2 m + 1 ) h 2 k × n − j + 1 − m 2 k − 1 α − n − j + 1 − m + 1 2 k − 1 α .

Removing the approximate solution from the perturbed equations yields

(25) u ˜ n + 1 = u ˜ 0 + h α Γ ( α + 1 ) ∑ j = 0 n ∑ m = 0 2 k − 1 − 1 γ t j + ( 2 m + 1 ) h 2 k , u j + ( 2 m + 1 ) h 2 k + u ˜ j + ( 2 m + 1 ) h 2 k − γ t j + ( 2 m + 1 ) h 2 k , u j + ( 2 m + 1 ) h 2 k × n − j + 1 − m 2 k − 1 α − n − j + 1 − m + 1 2 k − 1 α .

Taking the absolute value, we obtain

(26) ∣ u ˜ n + 1 ∣ ≤ ∣ u ˜ 0 ∣ + h α Γ ( α + 1 ) ∑ j = 0 n ∑ m = 0 2 k − 1 − 1 γ t j + ( 2 m + 1 ) h 2 k , u j + ( 2 m + 1 ) h 2 k + u ˜ j + ( 2 m + 1 ) h 2 k − γ t j + ( 2 m + 1 ) h 2 k , u j + ( 2 m + 1 ) h 2 k × n − j + 1 − m 2 k − 1 α − n − j + 1 − m + 1 2 k − 1 α .

The Lipschitz condition of the function γ with respect to u leads to

(27) ∣ u ˜ n + 1 ∣ ≤ ∣ u ˜ 0 ∣ + h α Γ ( α + 1 ) ∑ j = 0 n ∑ m = 0 2 k − 1 − 1 u ˜ j + ( 2 m + 1 ) h 2 k n − j + 1 − m 2 k − 1 α − n − j + 1 − m + 1 2 k − 1 α .

To continue our proof, we need to evaluate

(28) n − j + 1 − m 2 k − 1 α − n − j + 1 − m + 1 2 k − 1 α .

Employing the mean value theorem, we have

(29) n − j + 1 − m 2 k − 1 α − n − j + 1 − m + 1 2 k − 1 α = α ( α − 1 ) ξ α − 1 1 2 k − 1 .

By taking z ˜ = max 0 ≤ j ≤ n ∑ m = 0 2 k − 1 − 1 u ˜ j + ( 2 m + 1 ) h 2 k then

(30) ∣ u ˜ n + 1 ∣ ≤ ∣ u ˜ 0 ∣ + ( α − 1 ) h α Γ ( α ) z ˜ ξ α − 1 n 2 k − 1 ≤ c z ˜ ( n + 1 ) ,

where

(31) c = ( α − 1 ) h α Γ ( α ) 2 k − 1 ξ α − 1 .

This shows that the scheme is stable. We next want to show that

(32) lim Δ t → 0 ∣ u ( t n + 1 ) − u n + 1 ∣ = 0 .

Noting that

(33) u ( t n + 1 ) = u 0 + 1 Γ ( α ) ∫ 0 t n + 1 γ ( τ , u ( τ ) ) ( t n + 1 − τ ) α − 1 d τ

and

(34) u n + 1 = u 0 + h α Γ ( α + 1 ) ∑ j = 0 n ∑ m = 0 2 k − 1 − 1 γ t j + ( 2 m + 1 ) h 2 k , u j + ( 2 m + 1 ) h 2 k × n − j + 1 − m 2 k − 1 α − n − j + 1 − m + 1 2 k − 1 α .

Taking difference the above equations yields

(35) ∣ u ( t n + 1 ) − u n + 1 ∣ = 1 Γ ( α ) ∫ 0 t n + 1 γ ( τ , u ( τ ) ) ( t n + 1 − τ ) α − 1 d τ − 1 Γ ( α ) ∑ j = 0 n ∫ t j t j + 1 γ ( t j , u j ) ( t n + 1 − τ ) α − 1 d τ

and

(36) ∣ u ( t n + 1 ) − u n + 1 ∣ = 1 Γ ( α ) ∫ 0 t n + 1 γ ( τ , u ( τ ) ) ( t n + 1 − τ ) α − 1 d τ − ∑ j = 0 i ∫ t j t j + 1 γ ( t j , u j ) ( t n + 1 − τ ) α − 1 d τ ≤ ∫ 0 t n + 1 ∣ γ ( τ , u ( τ ) ) − γ ( t j , u j ) ∣ ( t n + 1 − τ ) α − 1 d τ .

Using the mean value theorem yields

(37) ∣ u ( t n + 1 ) − u n + 1 ∣ ≤ ∣ γ ′ ( c 1 , u ( c 1 ) ) ∣ ∫ 0 t n + 1 ( τ − t j ) ( t n + 1 − τ ) α − 1 d τ .

We calculate the above integral as follows:

(38) ∫ 0 t n + 1 ( τ − t j ) ( t n + 1 − τ ) α − 1 d τ = ( Δ t ) α + 1 Γ ( α + 2 ) ( ( n + 1 ) α ( n + 1 − j ( α + 1 ) ) ) .

Taking the limit of ∣ u ( t n + 1 ) − u n + 1 ∣ , we obtain

(39) lim Δ t → 0 ∣ u ( t n + 1 ) − u n + 1 ∣ = 0 .

4.1 Error analysis with power-law kernel

In this section, we present error analysis for the considered numerical scheme. For error analysis, we consider a non-linear fractional differential equation which is stated by

(40) D t α 0 C u ( t ) = γ ( t , u ( t ) ) , u ( 0 ) = u 0 .

To calculate the error, first, we suppose that the function γ ( t , u ( t ) ) ∈ C 2 ( [ a , b ] ) . Denote the partition points as t j = a + j h , t j + 1 = a + ( j + 1 ) h , and c j = a + 2 j + 1 2 k h the midpoints such that j = 0,1,2 , … , n − 1 . Before proceeding with Caputo case, we first start with classical case

(41) u ( t n + 1 ) − u ( t n ) = ∫ t n t n + 1 γ ( τ , u ( τ ) ) d τ = ∫ t n t n + 1 γ ( τ , u ( τ ) ) d τ = ∫ t n t n + 1 P n ( τ ) + ( τ − t n ) 2 2 ! ∂ 2 ∂ τ 2 γ ( τ , u ( τ ) ) d τ = ∫ t n t n + 1 P n ( τ ) + ( τ − t n ) 2 2 ! ∂ 2 ∂ τ 2 γ ( τ , u ( τ ) ) d τ + E n α ,

where

(42) E n α = ∫ t n t n + 1 ( τ − t n ) 2 2 ! ∂ 2 ∂ τ 2 γ ( τ , u ( τ ) ) d τ .

Then, we write

(43) E n α = 1 2 ∑ j = 1 n ∫ t j t j + 1 ( τ − c j ) 2 γ ″ ( c j , u ( c j ) ) d τ = 1 2 ∑ j = 1 n γ ″ ( c j , u ( c j ) ) ∑ m = 0 2 k − 1 ∫ t j + m h 2 k t j + ( m + 1 ) h 2 k ( τ − c j ) 2 d τ .

We evaluate

∑ m = 0 2 k − 1 ∫ t j + m h 2 k t j + ( m + 1 ) h 2 k ( τ − c j ) 2 d τ = ∑ m = 0 2 k − 1 h 3 3 × 2 3 k .

Then, for classical case, we obtain error as follows:

(44) E n α = 1 2 ∑ j = 1 n γ ″ ( ξ j , u ( ξ j ) ) 2 k h 3 3 × 2 3 k .

Considering h = b − a n , then we have the following:

(45) E n α = 1 6 ∑ j = 1 n γ ″ ( ξ j , u ( ξ j ) ) h 3 2 2 k ≃ γ ″ ( ξ , u ( ξ ) ) ( b − a ) h 2 6 × 2 2 k ,

where ξ ∈ ( a , b ) . When k = 1 which means we use midpoint, then we recover error for midpoint

(46) E n α ≃ ( b − a ) h 2 24 γ ″ ( ξ j , u ( ξ j ) ) .

Now, we continue with the Caputo case

(47) u ( t n + 1 ) − u ( 0 ) = 1 Γ ( α ) ∫ 0 t n + 1 γ ( τ , u ( τ ) ) ( t n + 1 − τ ) α − 1 d τ = 1 Γ ( α ) ∑ j = 0 n ∫ t j t j + 1 γ ( τ , u ( τ ) ) ( t n + 1 − τ ) α − 1 d τ = 1 Γ ( α ) ∑ j = 0 n ∫ t j t j + 1 P n ( τ ) + ( τ − t m ) 2 2 ! ∂ 2 ∂ τ 2 [ γ ( τ , u ( τ ) ) ] τ = γ m × ( t n + 1 − τ ) α − 1 d τ

= h α Γ ( α + 1 ) ∑ j = 0 n ∑ m = 0 2 n − 1 − 1 γ t j + ( 2 m + 1 ) h 2 n , u t j + ( 2 m + 1 ) h 2 n × n − j + 1 − m 2 n − 1 α − n − j + 1 − m + 1 2 n − 1 α + E n α ,

where

(48) E n α = 1 Γ ( α ) ∑ j = 0 n ∫ t j t j + 1 ( τ − c j ) 2 2 ∂ 2 ∂ τ 2 [ γ ( τ , u ( τ ) ) ] τ = γ m × ( t n + 1 − τ ) α − 1 d τ .

We evaluate

(49) E n α = 1 2 Γ ( α ) ∑ j = 0 n γ ″ ( c j , u ( c j ) ) ∫ t j t j + 1 ( τ − c j ) 2 ( t n + 1 − τ ) α − 1 d τ = 1 2 Γ ( α ) ∑ j = 0 n γ ′ ′ ( c j , u ( c j ) ) × ∑ m = 0 2 k − 1 ∫ t j + m h 2 k t j + ( m + 1 ) h 2 k ( τ − c j ) 2 ( t n + 1 − τ ) α − 1 d τ .

If we take norm on both sides of the above equation, we write the following inequality:

(50) ∣ E n α ∣ ≤ ( Δ t ) α + 2 2 Γ ( α + 3 ) sup τ ∈ [ 0 , t n + 1 ] ∂ 2 ∂ τ 2 γ ( τ , u ( τ ) ) ∑ m = 0 2 k − 1 ∫ t j + m h 2 k t j + ( m + 1 ) h 2 k ( τ − c j ) 2 ( t n + 1 − τ ) α − 1 d τ ,

where

(51) ∑ m = 0 2 k − 1 ∫ t j + m h 2 k t j + ( m + 1 ) h 2 k ( τ − c j ) 2 ( t n + 1 − τ ) α − 1 d τ = ( n + 1 ) α ( 6 n 3 + 18 n 2 + 18 n + 6 ) − n α 6 n 3 + 18 n 2 + 3 α 2 + 18 n + 9 α + 3 n α 2 + 6 n 2 α + 15 n α + 6 + 2 − k n α 6 n 3 α + 12 + 12 n 3 + 6 α 2 + 18 α + 36 n 2 + 36 n + 6 n 2 α 2 + 30 n 2 α + 12 n α 2 + 42 n α − ( n + 1 ) α 6 n 3 α + 12 + 12 n 3 + 6 α + 36 n 2 + 36 n + 18 n 2 α + 18 n α + 4 − k ( ( n + 1 ) α − n α ) ( 4 n 3 + 12 n 2 + 11 n + 3 ) × ( α 2 + 3 α + 2 ) .

Thus, the error can be obtained as

(52) ∣ E n α ∣ ≤ ( b − a ) h α + 1 2 Γ ( α + 3 ) γ ″ ( ξ j , u ( ξ j ) ) ( A + 2 − k B + 4 − k C ) ,

where

(53) A = ( n + 1 ) α ( 6 n 3 + 18 n 2 + 18 n + 6 ) − n α 6 n 3 + 18 n 2 + 3 α 2 + 18 n + 9 α + 3 n α 2 + 6 n 2 α + 15 n α + 6 , B = n α 6 n 3 ( α + 2 ) + n 2 ( 6 α 2 + 30 α + 36 ) + n ( 36 + 12 α 2 + 42 α ) + 12 + 6 α 2 + 18 α − 6 ( n + 1 ) α ( n + 1 ) 3 ( α + 2 ) , C = ( ( n + 1 ) α − n α ) ( 4 n 3 + 12 n 2 + 11 n + 3 ) × ( α 2 + 3 α + 2 ) .

5 Application of the proposed scheme to the Fitzhugh–Nagumo model

Simplifying the understanding of excitable media dynamics, the Fitzhugh–Nagumo equation is utilized to depict the behavior of action potentials in neurons and certain chemical reactions. They were proposed independently by Fitzhugh in 1961 [15] and Nagumo et al. in 1962 [16]. The solution of the Fitzhugh–Nagumo equation is of significant importance, as it allows for the acquisition of insights into the fundamental mechanisms underlying excitability in these cells. Therefore, in this section, we aim to solve this equation numerically. The associated model consists of two ordinary differential equations:

The first equation represents the membrane potential V of the excitable cell:
(54) V ′ ( t ) = V − V 3 3 − W + I .
The second equation represents the recovery variable W , which models the inactivation of ion channels:
(55) W ′ ( t ) = ε ( V + a − b W ) .

Here, the constant I represents an external stimulus and a , b , and ε are parameters. These equations capture the basic behavior of excitable systems, including the generation of action potentials and the refractory period. To achieve our aim, we first consider the Fitzhugh–Nagumo equation in the case of the Atangana–Baleanu fractional derivative

(56) D t α 0 A B V ( t ) = ε V − V 3 3 − W + I , D t α 0 A B W ( t ) = V + a − b W ,

where the initial conditions are V ( 0 ) = − 0.5 , W ( 0 ) = − 0.1 and the parameters are chosen as a = 0.001 , b = 0.002 , ε = 9 , and I = 0.32 . Applying the associated integral on the above equation and evaluation at the point t = t n + 1 , we have the following:

(57) V ( t n + 1 ) = V ( 0 ) + ( 1 − α ) γ 1 ( t n + 1 , V ( t n + 1 ) , W ( t n + 1 ) ) + α Γ ( α ) ∫ 0 t n + 1 γ 1 ( τ , V ( τ ) , W ( τ ) ) ( t n + 1 − τ ) α − 1 d τ W ( t n + 1 ) = W ( 0 ) + ( 1 − α ) γ 2 ( t n + 1 , V ( t n + 1 ) , W ( t n + 1 ) ) + α Γ ( α ) ∫ 0 t n + 1 γ 2 ( τ , V ( τ ) , W ( τ ) ) ( t n + 1 − τ ) α − 1 d τ ,

where

(58) γ 1 ( t , V , W ) = V − V 3 3 − W + I , γ 2 ( t , V , W ) = ε ( V + a − b W ) .

Employing the suggested method yields

(59) V n + 1 = V ( 0 ) + ( 1 − α ) γ 1 ( t n + 1 , V ˜ n + 1 , W ˜ n + 1 ) + h α Γ ( α ) ∑ j = 0 n ∑ m = 0 2 k − 1 − 1 γ 1 t j + ( 2 m + 1 ) h 2 k , V t j + ( 2 m + 1 ) h 2 k , W t j + ( 2 m + 1 ) h 2 k × n − j + 1 − m 2 k − 1 α − n − j + 1 − m + 1 2 k − 1 α

W n + 1 = W ( 0 ) + ( 1 − α ) γ 2 ( t n + 1 , V ˜ n + 1 , W ˜ n + 1 ) + h α Γ ( α ) ∑ j = 0 n ∑ m = 0 2 k − 1 − 1 γ 2 t j + ( 2 m + 1 ) h 2 k , V t j + ( 2 m + 1 ) h 2 k , W t j + ( 2 m + 1 ) h 2 k × n − j + 1 − m 2 k − 1 α − n − j + 1 − m + 1 2 k − 1 α .

The predictor components are calculated as

(60) V ˜ n + 1 = V 0 + ( 1 − α ) γ 1 ( t n , V n , W n ) + h α Γ ( α ) ∑ j = 0 n γ 1 ( t j , V j , W j ) × ( ( n − j + 1 ) α − ( n − j ) α ) , W ˜ n + 1 = W 0 + ( 1 − α ) γ 2 ( t n , V n , W n ) + h α Γ ( α ) ∑ j = 0 n γ 2 ( t j , V j , W j ) × ( ( n − j + 1 ) α − ( n − j ) α ) .

In Figures 3, 4, 5, we present the numerical simulation of the Fitzhugh–Nagumo equation with the Atangana–Baleanu fractional derivative.

$Figure 3 Numerical simulation of the Fitzhugh–Nagumo equation with the Atangana–Baleanu fractional derivative for α = 0.9 , k = 4 \alpha =0.9,k=4 .$

Figure 3

Numerical simulation of the Fitzhugh–Nagumo equation with the Atangana–Baleanu fractional derivative for α = 0.9 , k = 4 .

$Figure 4 Numerical simulation of the Fitzhugh–Nagumo equation with the Atangana–Baleanu fractional derivative for α = 0.9 , k = 4 \alpha =0.9,k=4 .$

Figure 4

Numerical simulation of the Fitzhugh–Nagumo equation with the Atangana–Baleanu fractional derivative for α = 0.9 , k = 4 .

$Figure 5 Numerical simulation of the Fitzhugh–Nagumo equation with the Atangana–Baleanu fractional derivative for α = 0.9 , k = 4 \alpha =0.9,k=4 .$

Figure 5

Numerical simulation of the Fitzhugh–Nagumo equation with the Atangana–Baleanu fractional derivative for α = 0.9 , k = 4 .

6 Conclusion

This study presents the derivation and analysis of the successive midpoint method for solving fractional differential equations involving the Atangana–Baleanu and Caputo derivatives. While the method has previously been applied to classical and Caputo–Fabrizio fractional models, this work extends its applicability to a broader class of fractional operators. A comprehensive error analysis was conducted, along with investigations into the method’s consistency and stability for both classical and Caputo-type equations.

Furthermore, numerical experiments were performed to evaluate the effectiveness of the successive midpoint method in comparison to the well-known midpoint method and L1 method. The results demonstrate that the successive midpoint method consistently outperforms the classical midpoint scheme and L1 method in terms of accuracy. Notably, its application to equations involving the Atangana–Baleanu derivative proved particularly advantageous, owing to the integration of a dedicated predictor mechanism tailored to the non-local and non-singular nature of this derivative. Overall, this work contributes to the growing body of the literature on numerical methods for fractional differential equations by offering a robust, accurate, and practically implementable technique. Given the complex and memory-dependent structure of fractional models, the successive midpoint method stands out as a promising tool for researchers and practitioners working with advanced fractional dynamics.

The successive midpoint method offers several advantages, including improved accuracy, consistency, and stability over the classical midpoint method. Its performance is especially notable for equations involving the Atangana–Baleanu derivative, as it incorporates a tailored predictor component that effectively handles the non-local and non-singular nature of this operator. This makes it a robust tool for modeling complex fractional systems with memory effects. However, the method also has some limitations. It is computationally more demanding due to its predictor–corrector structure and may require careful tuning of parameters, such as step size, to maintain efficiency.

Acknowledgments

The researchers would like to thank the Deanship of Graduate Studies and Scientific Research at the Qassim University for financial support (QUAPC2025).

Funding information: The financial support by the Deanship of Graduate Studies and Scientific Research at the Qassim University (QUAPC2025).
Author contributions: Maroua Amel Boubekeur: conceptualization, investigation; visualization; roles/writing – original draft; writing – review, and editing. Salah Boulaaras: conceptualization, investigation, methodology, software, supervision, visualization, writing – review, and editing. Seda İğret Araz: conceptualization, investigation, methodology, software, supervision, visualization, writing – review, and editing. All authors have accepted responsibility for the entire content of this manuscript and approved its submission.
Conflict of interest: The authors state no conflict of interests.
Data availability statement: All data generated or analysed during this study are included in this published article.

References

[1] Endre S, David M. An introduction to numerical analysis. Cambridge, United Kingdom: Cambridge University Press, ISBN 0-521-00794-1, 2003. Search in Google Scholar

[2] Ralston A. Runge–Kutta methods with minimum error bounds. Math Comput. 1962;92:339. 10.2307/2003133Search in Google Scholar

[3] Argyros IK, Chen D. The midpoint method for solving nonlinear operator equations in Banach space. Appl Math Lett. 1992;5:7–9. 10.1016/0893-9659(92)90076-LSearch in Google Scholar

[4] Toufik M, Atangana A. New numerical approximation of fractional derivative with non-local and non-singular kernel: application to chaotic models. Europ Phys J Plus. 2017;132:1–16. 10.1140/epjp/i2017-11717-0Search in Google Scholar

[5] Atangana A, Igret Araz S. New numerical scheme with newton polynomial: theory, methods, and applications. London, England: Elsevier, Academic Press; 2021. Search in Google Scholar

[6] Arık IA, Sarı HK, Igret Araz S. Numerical simulation of Covid-19 model with integer and non-integer order: The effect of environment and social distancing. Results Phys. 2023;51:106725. 10.1016/j.rinp.2023.106725Search in Google Scholar

[7] Iğret Araz S, Çetin MA. Fractal-fractional modeling of the Covid-19 spread with deterministic and stochastic approaches. Int J Appl Comput Math. 2025;11:4. 10.1007/s40819-024-01815-4Search in Google Scholar

[8] Yazdani AB, Kiasari MM, Jafari H. A novel approach for solving time fractional nonlinear Fisher’s equation by using Chebyshev spectral collocation method. Progress Fract Differ Appl. 2021;7(2):97–102. 10.18576/pfda/070204Search in Google Scholar

[9] Ali KK, Karakoç SBG, Rezazadeh H. Optical soliton solutions of the fractional perturbed nonlinear Schrodinger equation, TWMS J Appl Eng Math. 2020;10(4):930–9. Search in Google Scholar

[10] Karakoç SBG, Bhowmik SK, Omrani K, Yıldırım Sucu D. A numerical treatment for the Gilson-Pickering equation using collocation method with error estimation. Sigma J Eng Nat Sci. 2023;41(4):875–84. 10.14744/sigma.2023.00094Search in Google Scholar

[11] Atangana A, Igret Araz S. A successive midpoint method for nonlinear differential equations with classical and Caputo-Fabrizio derivatives. AIMS Math. 2023;8(11):27309–27. 10.3934/math.20231397Search in Google Scholar

[12] Podlubny I. Fractional Differential Equations, vol. 198 of Mathematics in Science and Engineering. San Diego: Academic Press; 1999. ISBN: 9780125588409. Search in Google Scholar

[13] Caputo M, Linear model of dissipation whoseQ is almost frequency independent. II. Geophys J Int. 1967;13(5):529–39. 10.1111/j.1365-246X.1967.tb02303.xSearch in Google Scholar

[14] Atangana A, Baleanu D New fractional derivatives with non-local and non-singular kernel. Theory Appl Heat Transfer Model Thermal Sci. 2016;20(2):763–9. 10.2298/TSCI160111018ASearch in Google Scholar

[15] FitzHugh R. Impulses and physiological states in theoretical models of nerve membrane. Biophys J. 1961;1(6):445–66. 10.1016/S0006-3495(61)86902-6Search in Google Scholar

[16] Nagumo J, Arimoto S, Yoshizawa S. An active pulse transmission line simulating nerve axon. Proc IRE. 1962;50(10):2061–70. 10.1109/JRPROC.1962.288235Search in Google Scholar

[17] Atangana A, Igret Araz S. Theory and methods of piecewise defined fractional operators. Elsevier. 2025. Search in Google Scholar

[18] Arnold DN. Stability, consistency, and convergence of numerical discretizations, encyclopedia of applied and computational mathematics. Berlin, Heidelberg: Springer; 2015. p. 1358–64. 10.1007/978-3-540-70529-1_407Search in Google Scholar

[19] Eslahchi MR, Esmaili S. The convergence and stability analysis of a numerical method for solving a mathematical model of language competition. Appl Numer Math. 2020;151:119–40. 10.1016/j.apnum.2019.12.015Search in Google Scholar

[20] Li M, Shi D, Pei L. Convergence and superconvergence analysis of finite element methods for the time fractional diffusion equation. Appl Numer Math. 2020;151:141–60. 10.1016/j.apnum.2019.12.023Search in Google Scholar

[21] Costa VAF. On the convergence of the iterative methods. Numer Heat Transfer Part B Fund. 2022;81:1–6, 1–13. 10.1080/10407790.2021.1997029Search in Google Scholar

[22] Wang J, Yang JZ, Zhang J. Stability and convergence analysis of high-order numerical schemes with DtN-type absorbing boundary conditions for nonlocal wave equations. IMA J Numer Anal. 2024;44(1):604–32. 10.1093/imanum/drad016Search in Google Scholar

[23] Franyelit SC, Luis RR. Convergence and stability criteria for numerical solutions of partial differential equations in science and engineering. Univ Ciencia Y Tech. 2024;28(122):7–16. 10.47460/uct.v28i122.761Search in Google Scholar

[24] Nasresfahani F, Eslahchi MR. Numerical convergence and stability analysis for a nonlinear mathematical model of prostate cancer. Numer Methods Partial Differ Equ. 2023;39(4):3064–88. 10.1002/num.23000Search in Google Scholar

Received: 2025-05-19

Revised: 2025-07-04

Accepted: 2025-07-28

Published Online: 2025-09-01

This work is licensed under the Creative Commons Attribution 4.0 International License.

Articles in the same Issue

https://doi.org/10.1515/phys-2025-0203

Keywords for this article

successive midpoint; stability; convergence; Fitzhugh–Nagumo model; fractional derivatives

Creative Commons

BY 4.0